Refuge-mediated predator–prey dynamics and biomass pyramids

Accepted Manuscript

Refuge-Mediated Predator-Prey Dynamics and Biomass Pyramids Hao Wang, Silogini Thanarajah, Philippe Gaudreau PII: DOI: Reference:

S0025-5564(17)30004-4 10.1016/j.mbs.2017.12.007 MBS 8014

To appear in:

Mathematical Biosciences

Received date: Revised date: Accepted date:

4 January 2017 9 December 2017 26 December 2017

Please cite this article as: Hao Wang, Silogini Thanarajah, Philippe Gaudreau, Refuge-Mediated Predator-Prey Dynamics and Biomass Pyramids, Mathematical Biosciences (2017), doi: 10.1016/j.mbs.2017.12.007

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights • Develop a mathematical model describing the predator-prey interactions when prey have access to a refuge. • Investigate the effects of the prey refuge on the predator-prey dynamics and the occurrence of an inverted biomass pyramid. • Study the impacts of fishing efforts and predator migrations on the predator-prey dynamics and the existence of an inverted biomass pyramid.

CR IP T

• Stability and bifurcation analysis

• Numerical test and sensitivity analysis for the predator-prey cycle in the refuge case.

AC

CE

PT

ED

M

AN US

• Observe that high fishing pressure or strong predator migration will change an inverted biomass pyramid into a regular biomass pyramid.

1

ACCEPTED MANUSCRIPT

Refuge-Mediated Predator-Prey Dynamics and Biomass Pyramids

CR IP T

Hao Wang∗, Silogini Thanarajah, Philippe Gaudreau Department of Mathematical and Statistical Sciences, University of Alberta Abstract.

M

AN US

Refuge can greatly influence predator-prey dynamics by movements between the interior and the exterior of a refuge. The presence of refuge for prey decreases predation risk and can have important impacts on the sustainability of a predator-prey system. The principal purpose of this paper is to formulate and analyze a refuge-mediated predator-prey model when the refuge is available to protect a portion of prey from predation. We study the effect of the refuge size on the biomass ratio and extend our refuge model to incorporate fishing and predator migration separately. Our study suggests that decreasing the refuge size, increasing the predator fishing, and increasing the predator emigration stabilizes the system. Here, we investigate the dependence of Hopf bifurcation on refuge size in the presence of fishing or predator migration. Moreover, we discuss their effects on the biomass pyramid and establish a condition for the emergence of an inverted biomass pyramid. We perform numerical test and sensitivity analysis to check the robustness of our results and the relative importance of all parameters. We find that high fishing pressure may destroy the inverted biomass pyramid and thus decrease the resilience of reef ecosystems. In addition, increasing the emigration rate or decreasing the immigration rate decreases the predator-prey biomass ratio. An inverted biomass pyramid can occur in the presence of a stable limit cycle.

Introduction

PT

1

ED

Keywords. predator-prey model with refuge, inverted biomass pyramid, Hopf bifurcation, effects of fishing, predator migration.

AC

CE

The study of effects of refuge use by prey on the dynamics of predator-prey interactions is identified as a vital issue in mathematical ecology [12, 14, 24]. Refuges have many possible influences on predatorprey population dynamics in which prey use refuges [5, 24]. One way for prey to diminish the risk of predation is to seek shelter. The existence of prey refuge can have significant effects on the coexistence of predators and their prey by preventing prey extinction [26]. These types of ecosystems are relatively newly discovered and have generated plenty of new research and open questions in ecology. In this paper, our focus and parametrization is the study of coral reefs, although our theoretical framework can be applied to other refuge cases such as burrows [4], trees [6], rock talus [15], and cliff faces [2]. Within the area, a refuge may fulfill as source habitats while the refuge’s exterior may fulfill as a sink habitats. Thus, part of the refuge in predator-prey interactions may strengthen source-sink dynamics in both prey and predators. In most ecosystems, the biomass pyramid is widest at the base. This appears reasonable due to energy loss across tropic cascades. At first glance, inverted biomass pyramids seem quite paradoxical. However, such pyramids have been found in pristine coral reefs [29]. Previous studies have examined the effect of refuge through theoretical models and laboratory experiments and suggested that refuges are vital for prey persistence [13, 16]. Even though a number of theoretical studies have designated the effects of refuge on predator-prey coexistence, assumptions of these models are commonly based on mathematical perspective rather than on biological phenomena. Previous studies [32] assumed that a constant number or proportion of prey is found in refuge. In ∗ Corresponding

author: [email protected]

2

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

reality, the constant number of prey in refuge often varies with predation risk and resource level. They concluded that the presence of constant number of prey protected by refuge has a stronger stabilizing effect on population dynamics than that of protecting a constant proportion of prey [11, 24, 32]. Smith et al. [30] showed that the presence of a constant proportion of prey protected by refuge does not change the stability of the Lotka-Volterra system. Hassel et al. [12] presented that adding a large refuge to a model excludes the stable limit cycle. These experimental and mathematical results indicate that refuge has a stabilizing effect on predator-prey dynamics. In spite of these results, other work question its validity [8]. These conflicting conclusions as well as the interactions between populations inside and outside the refuge in ecological systems need further exploration. In order to model the emergence of inverted biomass pyramids, Singh et al. [29] derived a refuge-modulated predator-prey model. That model was an extended Lotka-Volterra model with a new predation function. In reality, the habitat should have two separate regions: interior refuge and exterior refuge, which have been included in some previous models [3, 7]. Our model also incorporates prey inside and outside refuge. The main feature of our model is the explicit incorporation of the refuge size, which allows us to explore the dependence of predator-prey dynamics and the biomass ratio on the refuge size. The prey carrying capacity inside refuge and the prey immigration into refuge are naturally assumed to be dependent on the refuge size. In this paper, we introduce a new mathematical model that modifies previous models to better describe reality. We extend our model to include fishing and predator migration to explore the effect of refuge. First, we study the effects of refuge on predator-prey dynamics and on biomass pyramids. Secondly, we investigate the impacts of fishing and predator migration on refuge dependences. Furthermore, we apply the qualitative methods of differential equations and use the Hopf bifurcation theorem to study refuge-mediated predator-prey dynamics. We demonstrate the existence of a Hopf bifurcation point and the periodic solution of the Hopf type. We study the existence of equilibrium and explain the emergence of inverted biomass pyramids and their dependence on the refuge size. We focus more on the internal equilibrium, of the model, that can generate the limit cycle via Hopf bifurcation. Also, we examine the stability conditions of all equilibria and the existence of a limit cycle. We explain the emergence of inverted biomass pyramids from our theoretical results. Our model supports previous research which has found that a larger refuge size yields higher biomass ratios. Further, we observe that the introduction of fishing and predator migration decreases the biomass ratio, and sufficiently strong fishing and predator emigration changes an inverted biomass pyramid into a regular biomass pyramid. Finally, the presence of refuge for part of the prey population replaces the stable internal equilibrium with a stable limit cycle. Many factors can influence predator-prey dynamics in coral reefs. One of these important factors is fishing. Another factor which affects dynamical properties is predator migration. Effects of fishing on refuge-mediated prey-predator models have been considered by a few researchers but effects of predator migration on refuge-mediated prey-predator models have not been considered in any research yet. Kar et al. [20] showed that the presence of independent fishing on predators and prey break the cyclic behavior of the system. Ji et al. [19] considered a predator-prey model with constant number of prey protected by refuge and a constant rate of prey fishing and showed that the system is controlled by constant number of prey protected by refuge or constant rate of prey fishing. Here, we consider the combined effect of fishing and prey refuge and the combined effect of predator migration and prey refuge on predator-prey dynamics. The main questions addressed in this paper are the following: What is the dependence of predatorprey dynamics on the refuge size? How does fishing or predator-migration regulate its dependence? How does the biomass pyramid depend on the prey refuge? The explicit incorporation of the prey refuge in our model allows us to answer these questions directly. The paper is organized as follows. In section 2, we study the dynamics of the refuge-mediated predator-prey model, the occurrence of an inverted biomass pyramid, periodic solutions, and sensitivity analysis. In section 3, we incorporate fishing into our model and study the refuge effect on the system. In section 4, we incorporate predator emigration and immigration into our model and study the refuge effect on the system. We conclude and discuss our results in section 5.

2

A Refuge-Mediated Predator-Prey Model

Coral reefs provide homes to many fish species. The diversity of fish found on reefs ranges from large predators that prey on other fish to small prey that eat plankton. We consider fish biomass directly in 3

ACCEPTED MANUSCRIPT

CR IP T

our model rather than through a proxy of fish population. It is assumed that the habitat is divided into two areas, namely, inside the refuge and outside the refuge. It is assumed that predator species are not allowed to enter the refuge. The presence of refuge manages prey biomass dynamics in two different ways. In the following section, we introduce the Refuge-Predator-Prey model, its variables and the assumptions that are made for its development. Singh et al. [29] constructed a mathematical model with an explicit refuge to represent a new biologically reasonable mechanism that can explain stable inverted biomass pyramids in pristine coral reefs. They assumed that the prey growth rate depends on the refuge size, i.e., a(r). Wang et al. [33] originally proposed three refuge-mediated predation functional responses with their underlying biological mechanisms. We construct a more realistic model with the assumption that the prey carrying capacity in the refuge and the prey migration rate into the refuge explicitly depend on the refuge size, i.e., K(r),β(r) respectively. It is assumed that the prey species in the refuge grow logistically but never grow outside the refuge because most resources are stored inside the refuge. For the development of this model, we make the following key assumptions: • Prey can only find food and grow in the refuge, and the carrying capacity is an increasing function of the refuge size. • The prey growth inside the refuge is logistic.

AN US

• In the absence of prey, the predator decays exponentially. • Predators can only eat prey outside the refuge.

• Prey can migrate in or out of the refuge, and only the immigration rate depends on the refuge size.

PT

ED

M

The predator-prey model in coral reef is provided by   xr dxr = axr 1 − − αxr + β(r)x, dt K(r) dx (1) = αxr − β(r)x − f (x)y, dt dy = cf (x)y − dy, dt where K(r), β(r) are increasing nonnegative continuous functions. We also assume that the parameters a, d, α, c > 0. There are many forms of the predator functional response in the literature. We assume that f (x) and f 0 (x) are both positive for x > 0 and that f (0) = 0. For instance, we have the Holling γx type I (f (x) = γx) and the Holling type II (f (x) = ∆+x ), where γ, ∆ are positive constant. The variables and functions in the model (1) are all positive and defined bellow:

CE

• xr - prey density in the refuge.

• x - prey density outside the refuge.

AC

• y - predator density.

• K(r) - carrying capacity of prey in the refuge.

• a - per-capita growth rate of prey inside refuge. • α - emigration rate from the refuge.

• β(r) - immigration rate into the refuge.

• f (x) - predator response function.

• c ∈ (0, 1) - conversion efficiency.

• d - density independent per-capita death rate of predator. • r - refuge size for prey (kg/m2 ).

• γ - per capita predation rate (in type I) or the maximum predation rate (in type II). 4

ACCEPTED MANUSCRIPT

2.1

Mathematical Analysis

In this section, we study the qualitative behavior of our refuge-dependent predator-prey model, including the positivity of solutions, the existence and local stability of all steady states, and the existence of periodic solutions of Hopf type. First, we establish positivity for this model. Otherwise stated, given positive initial conditions, the solutions will remain positive for all time. Theorem 2.1. All solutions of the system (1) which start in R3+ = [0, ∞)3 will remain nonnegative. dxr dt xr =0 dx dt x=0 dy dt y=0 Hence, all solutions will remain positive for all

=

CR IP T

Proof. It is easy to verify that:

β(r)x ≥

0

=

αxr

≥

0

=

0

≥

0

times.

(2)

AN US

We now study the existence of the steady states. Model (1) has only three nonnegative equilibria, namely E0 (0, 0, 0), E1 (˜ xr , x ˜, 0) and E∗ (¯ xr , x ¯, y¯). E0 (0, 0, 0)

=

E1 (˜ xr , x ˜, 0)

=

(0, 0, 0),   αK(r) ,0 , K(r), β(r)

y¯

ED

M

as well as one internal steady state denoted by E∗ (¯ xr , x ¯, y¯)) where s "  #   K(r) d β(r) −1 2 x ¯r = (a − α) + (a − α) + 4a K(r) f , 2a c
x ¯ = f −1 dc , =

c α¯ xr − β(r)f −1 d

d c




(3) (4)

(5)

.

(6)

CE

PT

The equilibria E0 and E1 always exist, and we show the existence of E∗ as follows: E∗ are feasible if and only if     d α¯ xr −1 α¯ xr − β(r)f > 0 ⇔ d < cf . c β(r)

AC

The steady state in equation (3) corresponds to the total extinction equilibrium where both prey and predators go extinct. The steady state in equation (4) correspond to the predator free steady state. The last steady state in equation (5) corresponds to the coexistence equilibrium. The Jacobian of this system is given by     2a (a − α) − x β(r) 0 r   K(r) . J(xr , x, y) =  (7) 0  α −β(r) − f (x)y −f (x)  0 cf 0 (x)y cf (x) − d We perform a local stability analysis of all three steady states shown in equations (3), (4) and (5). However, to perform stability analysis, we require the Routh-Hurwitz criterion. Theorem 2.2 (Routh-Hurwitz criterion). Given a second degree polynomial of the form P2 (λ) = λ2 + a1 λ + a2 , all the roots of the polynomial P2 (λ) are negative or have negative real parts if and only if 5

(8)

ACCEPTED MANUSCRIPT

• a1 > 0, • a2 > 0. Given a third degree polynomial of the form P3 (λ) = λ3 + a1 λ2 + a2 λ + a3 ,

(9)

all the roots of the polynomial P3 (λ) are negative or have negative real parts if and only if

CR IP T

• a1 > 0, • a3 > 0, • a1 a2 > a3 .

Using Theorem 2.2, we have all the tools necessary to establish stability properties of our three steady states (3), (4) and (5). We begin by studying the total extinction equilibrium.

AN US

Theorem 2.3. The extinction steady state E0 = (0, 0, 0) is a saddle with one-dimensional unstable manifold. Proof. Evaluating the Jacobian at the point E0 = (0, 0, 0), we  a − α β(r) −β(r) J(0, 0, 0) =  α 0 0 The characteristic polynomial of this system is given by

obtain  0 0 . −d

(10)


P (λ) = (λ + d) λ2 + (β (r) − a + α) λ − aβ (r) .

M

(11)

AC

CE

PT

ED

Hence, we deduce that the eigenvalues satisfy: λ1 = −d < 0 and λ2 λ3 = −aβ(r) < 0. If λ2 < 0, then λ3 > 0, which implies that E0 is a saddle with one-dimensional unstable manifold.   αK(r) , 0 is locally asymptotically stable Theorem 2.4. The predator free steady state E1 = K(r), β(r) if   αK (r) cf < d. (12) β (r)   Proof. Evaluating the Jacobian at the point (xr , x, y) = K(r), αK(r) , 0 , we obtain β(r)  a−α   αK(r) α ,0 =  J K(r),  β(r) 0 

β(r)

−β(r) 0

0





 −f αK(r) . β(r)    αK(r) cf β(r) − d

(13)

The characteristic polynomial of this system is given by    
αK (r) − d) λ2 + (β (r) + a + α) λ + β (r) a . (14) P (λ) = λ − (cf β (r)   −d is an eigenvalue. This eigenvalue is negative From equation (14), we can deduce that λ = cf αK(r) β(r) if and only if   αK (r) cf < d. (15) β (r)

6

ACCEPTED MANUSCRIPT

Using the Routh-Hurwitz stability criterion shown in Theorem 2.2, we know that the second degree polynomial in equation (14) has roots with negative real parts if and only if β (r) + a + α

> 0,

(16)

β(r)a

> 0.

(17)

when cf

αK(r) β(r)

> d.

CR IP T

These conditions are always will have negative real parts.   satisfied, hence both roots of this polynomial   αK(r) Hence K(r), < d and unstable (saddle) , 0 is locally asymptotically stable when cf αK(r) β(r) β(r)  

Theorem 2.5. The internal steady state E∗ = (¯ xr , x ¯, y¯) exists if and only if x ¯> locally asymptotically stable if

β(r) −1 α f

  d . E∗ is c

(d − v1 )v22 + (β(r)(α + d) − v12 )v2 + αβ(r)v1 > 0, where =

v2

=



   β(r) d f −1 , K(r) c    d 0 −1 y¯. −β(r) − f f c

AN US

v1

s

− (a − α)2 + 4a

ED

M

Proof. Evaluating the Jacobian at the point E∗ = (¯ xr , x ¯, y¯), we obtain r   
2 + 4a β(r) f −1 d − (a − α) β(r) 0 K(r) c  
J(¯ xr , x ¯, y¯) =  d  α −β(r) − f 0 (f −1 dc )¯ y − c
0 cf 0 (f −1 dc )¯ y 0 The characteristic polynomial of this Jacobian is given by

CE

where

PT

P (λ) = λ3 + A1 λ2 + A2 λ + A3 ,

A1 A2 A3

(18)

(19) (20)



  . 

(21)

(22)

= −(v1 + v2 ),

(23)

= dv1 (β(r) + v2 ),

(25)

= −(β(r)(d + α) + v2 (d − v1 )),

(24)

AC

It is important to note that both v1 < 0 and v2 < 0. Using the Routh-Hurwitz stability criterion shown in Theorem 2.2, we know that equation (22) will have roots with negative real part if and only if −(v1 + v2 ) > 0,

dv1 (β(r) + v2 ) > 0, (−(v1 + v2 ))(−(β(r)(d + α) + v2 (d − v1 )))

>

dv1 (β(r) + v2 ).

(26) (27) (28)

From the definition of v1 and v2 in equation (19) and (20), we can deduce that conditions (26) and (27) will always be satisfied. Moving the right hand side term to the left hand side in (28) and simplifying we obtain the condition (d − v1 )v22 + (β(r)(α + d) − v12 )v2 + αβ(r)v1 > 0.

7

(29)

ACCEPTED MANUSCRIPT

Biologically, if the value of prey refuge size increases, then the equilibrium density of predator population decreases due to lack of prey resources.

a α β γ ∆ d K c

= 0.007 = 0.035 = 0.0119 = 0.0112 = 1 = 0.00007 = 0.8 = 0.04.

CR IP T

For the suitable values of parameters, we verify the existence and stability of the internal equilibrium for the system as well as the existence of a stable limit cycle. In numerical simulations, we assume K(r) = Kr and choose the following parameter values:

(30)

AN US

The initial conditions are xr (0) = 1/10, x(0) = 1/20 and y(0) = 1/6. We vary the refuge size (r) from 0.01 to 1 to study the transition of dynamics in Figures 1 - 3. When the refuge size is very small, all solutions tend to the predator free equilibrium E1 (see Figure 1). If we slightly increase the refuge size, solutions tend to a stable internal node (see Figure 2). If we further increase the refuge size, solutions tend to a stable internal spiral (see Figure 3). If the refuge size is large enough, the system (1) has an unstable internal equilibrium point and an asymptotically stable limit cycle (see Figure 4).

M

0.015

0.005

PT

0 1

ED

y

0.01

0.5

0

0.2

0.1

0

0.3

0.4

xr

CE

x

E1

AC

Figure 1: The phase portrait of the system (1) with r = 0.01. All solutions directly converge to the predator free equilibrium E1 .

8

ACCEPTED MANUSCRIPT

0.2 0.15 E*

y

0.1

0 1 0.5 0

x

0.1

0

CR IP T

0.05

0.2

0.3

0.4

xr

AN US

Figure 2: The phase portrait of the system (1) with r = 0.3. All solutions directly converge to the internal equilibrium E∗ .

0.4

M

y

0.3 0.2

0 1.5

E*

ED

0.1

1

PT

0.5

x

0

0.2

0.1

0

0.3

0.4

xr

AC

CE

Figure 3: The phase portrait of the system (1) with r = 0.5. All solutions spirally converge to the internal equilibrium E∗ .

9

ACCEPTED MANUSCRIPT

0.5 0.4

y

0.3 0.2 0.1

CR IP T

E*

0 3 2 1 0

x

0.4

0.2

0

xr

0.6

0.8

AN US

Figure 4: The solution curves of the system (1) with r = 1. There is a stable periodic solution around the unstable internal spiral E∗ . The Existence and Properties of the Periodic Solution of Hopf Type We consider the existence and properties of the periodic solutions generated from Hopf bifurcation when f (x) = γx. In this case, we have s

v1

= −

v2

= − (β(r) + γ y¯) .

M

(a −

α)2

+ 4a



β(r) K(r)



 d , cγ

ED

By substituting v2 into (23), (24), (25), we obtain A1 , A2 , and A3 as follows: A1

PT

A2 A3

= −(v1 − β(r) − γ y¯),

(31)

= −dv1 γ y¯,

(33)

= −β(r)(d + α) + (β(r) + γ y¯)(d − v1 ),

(32)

CE

Theorem 2.6. Suppose that the coefficients of P (λ) = 0 at (¯ xr , x ¯, y¯) are functions of r and that Aj (r)(j = 1, 2, 3) at r = r¯ satisfy the following conditions: Aj (¯ r) > 0,

j = 1, 2, 3

AC

d(A1 A2 − A3 ) |r=¯r 6= 0. dr Then the Hopf bifurcation from (¯ xr , x ¯, y¯) occurs at r = r¯. As the value of r passes through r¯, the periodic solution of the Hopf type of system (1) appears in the neighborhood of (¯ xr , x ¯, y¯). A1 A2 − A3 |r=¯r = 0,

10

ACCEPTED MANUSCRIPT

Y 1

0.8

0.4

0.2

0 0.2

0.4

0.6

0.8

1 r

1.2

1.4

1.6

AN US

0

CR IP T

0.6

1.8

2

Figure 5: The supercritical Hopf bifurcation diagram for predator vs refuge size. Note:



 αK(r) There is no periodic solution of the Hopf type bifurcating from λ0 = (0, 0, 0) and λ1 = K(r), ,0 β(r) for any parameter because P (λ)|λ0 = 0 or P (λ)|λ1 = 0 cannot have a pair of pure imaginary roots.

AC

CE

PT

ED

M

The Jacobian matrix J(¯ xr , x ¯, y¯) has a single pair of complex conjugate eigenvalues η(r) = Re(η) + ¯ = Re(η) − Im(η). Suppose that this pair of eigenvalues has the largest real part of all Im(η), η(r) eigenvalues and is such that in a small neighborhood of a bifurcation value r¯, (i) Re(η) < 0 if r < r¯, (ii) Re(η) = 0, Im(η) 6= 0 if r = r¯, (iii) Re(η) > 0 if r > r¯. Then, in a small neighborhood of r¯ , r > r¯ , the steady state is unstable by growing oscillations and at least a small amplitude limit cycle periodic solution exists around the equilibrium point. The appearance of a periodic solution from an equilibrium is called a Hopf bifurcation. When the parameter r is continuously varied near the criticality r¯, a periodic solution emerges for r < r¯ (this case is referred to as a supercritical bifurcation) or for r > r¯ (which is referred to as a subcritical bifurcation). Figure 5 shows the Hopf bifurcation of the system (1) with refuge size (r) as the bifurcation parameter. The theorems for the existence and stability of steady states provide information to possible predatorprey dynamics in comparison to the classical predator-prey dynamics. The Hopf bifurcation theorem guarantees the appearance of limit cycles. These mathematical results give insights into the predatorprey dynamics, while the main scientific insights come from the biomass ratio in the next few sections.

2.2

Predator-Prey Biomass Ratio

In this section, we provide sufficient conditions for the emergence of inverted biomass pyramids of the system (1). Our model has three equilibrium points. The unstable equilibrium point (0, 0, 0) corresponds to a reef with no fish. The equilibrium point (K(r), αK(r) β(r) , 0) corresponds to the absence of predators γx and is rarely seen in reefs. The inner equilibrium for a Holling type II functional response, f (x) = ∆+x

11

ACCEPTED MANUSCRIPT

is given by x ¯r

=

r     K(r) β(r) d∆ (a − α) + (a − α)2 + 4a K(r) cγ−d , 2a

x ¯

=

d∆ cγ−d ,

(34)

i ch d∆ α¯ xr − β(r) cγ−d . d The predator-prey biomass ratio at the inner equilibrium point is defined by =

d∆ c xr − β(r) cγ−d ] d [α¯ . d∆ ¯r cγ−d + x

3.5 Type I Type II Type III TypeIV

2.5 2 1.5

0 1

M

1 0.5

(35)

AN US

3

1.2

ED

Predator−prey Biomass Ratio

y¯ = x ¯+x ¯r

CR IP T

y¯

1.4

1.6

Refuge Size

1.8

2

PT

Figure 6: The biomass pyramid is inverted and the predator-prey biomass ratio at an inner equilibrium point is an increasing function of refuge size. 2 √ βr Type I: β(r) = βr, Type II: β(r) = ξ+r , Type III: β(r) = ξ2βr +r 2 , Type IV: β(r) = β r.

AC

CE

Since there exist no empirical data or evidence for the function β(r), we test different possible types in Figure 6, which illustrates the dependance of the predator-prey biomass ratio on the refuge size for type I, type II, type III and type IV functional forms for β(r). The predator-prey biomass ratio is an increasing function of refuge size for all four functional forms, and the result on biomass ratio is robust with respect to the format of β(r). For simplicity, we use type I functional response for β(r). We will investigate the sufficient conditions for the existence of an inverted biomass pyramid. Theorem 2.7. Given the model (1) with the proper conditions satisfied for the stability of the internal steady state shown in (5), there will be an inverted biomass pyramid if the following conditions are satisfied:
!   2af −1 dc cβ(r) + d − (a − α) K(r) cα − d d r < 1. (36) <α and  
c β(r) d 2 −1 (a − α) + 4a K(r) f c

Proof. For the existence of an inverted biomass pyramid, we require that   y¯ c α¯ xr − β(r)¯ x 1< = x ¯+x ¯r d x ¯+x ¯r 12

(37)

ACCEPTED MANUSCRIPT

Rearranging equation (37), we obtain: x ¯r − β(r) x ¯ x ¯r +1 x ¯ α + β(r) α− x ¯r +1 x ¯ α + β(r) x ¯r +1 x ¯ 1 x ¯r +1 x ¯ α

<

⇒

d c

<

⇒

cα − d c

>

⇒

cα − d c(α + β(r))

>

(38)

CR IP T

d c

> 0

AN US

Equation (38) establishes a condition for the existence of an inverted biomass pyramid. For the emergence of an inverted biomass pyramid, we require that cα − d > 0. Assuming this to be true and isolating for x ¯r /¯ x, we obtain: cβ(r) + d x ¯r (39) > x ¯ cα − d Hence, r    
β(r) K(r) (a − α) + (a − α)2 + 4a K(r) f −1 dc cβ(r) + d (40)
> cα − d 2af −1 dc Simplifying, we obtain

d c




>

M

r   β(r) (a − α)2 + 4a K(r) f −1

ED

r   β(r) f −1 Dividing both sides by (a − α)2 + 4a K(r)

2af −1

K(r)

d c




d c


!

cβ(r) + d cα − d



− (a − α)

(41)

, we obtain our desired result.

2.3

PT

This theorem provides more realistic conditions on the refuge size for the appearance of an inverted biomass pyramid in comparison to a few cases given in Wang et al. [33]. The first condition dc < α is similar to the classical condition in Result 1 of Wang et al. [33], and the second condition shows a very nonlinear and complicated dependence on the refuge size.

Numerical Test

AC

CE

We will now perform some numerical test on a specific case of the model presented in (1). In this section, we focus on the following model:   dxr xr = axr 1 − − αxr + β(r)x dt K(r) dx γ xy (42) = αxr − β(r)x − dt ∆+x dy cγ xy = − dy. dt ∆+x

As we can see, in the model (42), we make the following assumptions regarding the function f (x). For the predator response function f (x), we assume a Holling type 2 function. Utilizing all the theorems presented in section 2.1, we can arrive at the following conclusions. The model (42) has three steady states. (xr , x, y) (xr , x, y)

=

(0, 0, 0)   αK(r) = K(r), ,0 β(r) 13

(43) (44)

ACCEPTED MANUSCRIPT

as well as one internal steady state denoted by (¯ xr , x ¯, y¯) where: s " #   d∆   K(r) β(r) x ¯r = (a − α) + (a − α)2 + 4a K(r) 2a γc − d x ¯

=

y¯

=

d∆ γc − d    c d∆ α¯ xr − β(r) , d γc − d

(45)

0<

d∆ αK(r) < γc − d β(r)

CR IP T

which is feasible only if

d γαK(r) < . c ∆β(r) + αK(r)

=⇒

(46)

The steady state (43) is always unstable. The steady state (44) is stable only if d γαK(r) < . ∆β + αK(r) c

AN US

Moreover, the internal steady state is stable if

(d − v1 )v22 + (β(r)(α + d) − v12 )v2 + αβ(r)v1 > 0, where

v2



β(r) = − (a − + 4a K(r)   (γc − d)2 = −β(r) − y¯. ∆γc2 α)2



(48)

 d∆ , γc − d

(49) (50)

M

v1

s

(47)

PT

ED

Additionally, the required conditions for the emergence of an inverted biomass pyramid are given by:     2a d∆ cβ(r) + d − (a − α) d K(r) γc − d cα − d s α> and < 1. (51)    d∆  c β(r) 2 (a − α) + 4a K(r) γc − d

CE

Since the scale of these parameters cab vary significantly, the visualizations of these solutions becomes troublesome. For this reason, we propose the following nondimensionalization of the model (42).

AC

xr (t)

= K(r)Xr (t)

t

=

e3

=

τ a γ a

x(t)

= K(r)X(t)

e1

=

e4

=

α a ∆ K(r)

y(t)

= K(r)Y (t)

e2

=

e5

=

β(r) a d a

(52)

Although not unique, we choose this specific change of variable in order for the biomass ratio to remain invariant up to a scaling in time i.e. Biomass Ratio(t)

y(t) xr (t) + x(t) y(t)/K(r) = xr (t)/K(r) + x(t)/K(r) Y ( τa ) = Xr ( τa ) + X( τa ) = Normalized Biomass Ratio( τa ) =

14

(53)

ACCEPTED MANUSCRIPT

Applying the changes of variables in system: dXr dτ dX dτ dY dτ

(52) to the model in (42), we arrive at the nondimensionalized = Xr [1 − Xr ] − e1 Xr + e2 X   XY = e1 Xr − e2 X − e3 e4 + X   XY = ce3 − e5 Y. e4 + X

(54)

(Xr , X, Y )

=

(Xr , X, Y )

=

CR IP T

With this in mind, we find the steady states for this model to be given by (0, 0, 0)   e1 1, , 0 e2

(55) (56)

¯ r , X, ¯ Y¯ ) where: as well as one internal steady state denoted by (X =

¯ X

=

Y¯

=

i p 1h ¯ (1 − e1 ) + (1 − e1 )2 + 4e2 X 2

AN US

¯r X

e4 e5 ce3 − e5 ¯ r − e2 X)(e ¯ 4 + X) ¯ (e1 X , ¯ e3 X

(57)

e4 e5 e1 < . ce3 − e5 e2

(58)

which is feasible only if

M

0<

CE

PT

ED

The most interesting visual aspect of this model arise when plotting the solutions of the nondimensionalize model (54) when the internal steady state (57) is unstable. In such instances, we obtain behaviour similar to that shown in Figures 7 and 8. The values used to create Figure 7 are the following: e1 e2 e3 e4 e5 c

= = = = = =

5 1.7 1.6 1.25 0.01 0.04.

(59)

3 0.9 4 1.1 0.07 0.1.

(60)

AC

The values used to create figure 8 are the following: e1 e2 e3 e4 e5 c

= = = = = =

Both cases had the initial conditions, Xr (0) = 1/10, X(0) = 1/20 and Y (0) = 1/6. The case in Figure 7 has the conditions for the emergence of an inverted biomass pyramid satisfied. The case in Figure 8 does not have the conditions for the emergence of an inverted biomass pyramid satisfied. Although not proven in this paper, it would appear as though there is an emergence of a stable limit cycle as seen in both Figures 7 (a) and 8 (a). In Figure 7 (b) and 8 (b), the normalized biomass ratio in equation (53) is plotted over time in terms of units of the per-capita growth rate of prey inside refuge. When this ratio is above the line, this ratio is greater than one which implies the presence of an inverted biomass pyramids. Oppositely, when this ratio is bellow the line, the biomass pyramid is regular. As we can see 15

ACCEPTED MANUSCRIPT

Predator−Prey interaction over time

1

10

0

10

−1

10

−2

10

−3

10

0

200

400 600 time ( units of a )

(a)

800

1000

AN US

Predator−Prey Ratio over time

CR IP T

Normalized population density

Normalized prey population in refuge Normalized prey population outside refuge Normalized predator population

2

10

1

M

10

0

10

ED

Predator−Prey Ratio

Predator−Prey Ratio over time Biomass Ratio index

−1

0

PT

10

200

400 600 time ( units of a )

800

1000

(b)

AC

CE

Figure 7: Panel (a) displays the normalized population densities for of the nondimensionalize system (54) over time in terms of units of the per-capita growth rate of prey inside refuge. Panel (b) displays the ratio of the normalised population densities over time in terms of units of the per-capita growth rate of prey inside refuge. In this case the conditions for the emergence of an inverted biomass pyramid are satisfied. from both Figures 7 and 8, there appears to be an emergence of a inverted biomass pyramid in both cases. This would imply that an inverted biomass pyramid can occur in the presence of a stable limit cycle however, they are not permanently sustainable. This observation leads to many interesting questions. We are interested in determining if this limit cycle is spent mostly as an inverted biomass pyramid or a regular biomass pyramid. To answer this question, we propose the following measure.   Z y(s) 1 t+T ds, (61) MeanRatio = lim t→∞ T t x(s) + xr (s) where T is the period of the limit cycle. The quantity M eanRatio gives the average value 16

y(s) x(s) + xr (s)

ACCEPTED MANUSCRIPT

Predator−Prey interaction over time

0.9

Normalized prey population in refuge Normalized prey population outside refuge Normalized predator population

0.8

0.6 0.5 0.4 0.3 0.2 0.1 0

0

200

400

600

800 1000 1200 time ( units of a )

1400

(a) Predator−Prey Ratio over time 3.5

3

1800

2000

2

1.5

M

Predator−Prey Ratio

2.5

1

ED

0.5

0

200

400

600

800 1000 1200 time ( units of a )

1400

1600

1800

2000

(b)

PT

0

1600

AN US

Predator−Prey Ratio over time Biomass Ratio index

CR IP T

Normalized population density

0.7

CE

Figure 8: Panel (a) displays the normalized population densities for of the nondimensionalize ODE system (54) over time in terms of units of the per-capita growth rate of prey inside refuge. Panel (b) displays the ratio of the normalised population densities over time in terms of units of the per-capita growth rate of prey inside refuge. In this case the conditions for the emergence of an inverted biomass pyramid are not satisfied.

AC

over one period of the limit cycle. If MeanRatio > 1 then the majority of the limit cycle period is spent as an inverted biomass pyramid. Conversely, if MeanRatio < 1, then the majority of the period of one limit cycle is spent as a regular biomass pyramid. Other quantities of interest are the maximum and minimum values of the steady limit cycle. These quantities can be expressed as follows: y(t) x(t) + xr (t)

(62)

y(t) . x(t) + xr (t)

(63)

MaxLimitCycle = lim sup t→∞

MinLimitCycle = lim inf t→∞

17

ACCEPTED MANUSCRIPT

Their normalized counterparts may be written as: NormalizedMeanRatio

=

NormalizedMaxLimitCycle

=

NormalizedMinLimitCycle

=

1 ˜ τ →∞ T

Z

τ +T˜



Y (s) X(s) + Xr (s) τ Y (τ ) lim sup , τ →∞ X(t) + Xr (τ ) Y (τ ) lim inf . τ →∞ X(τ ) + Xr (τ ) lim



ds,

(64) (65) (66)

T˜ NormalizedMeanRatio NormalizedMaxLimitCycle NormalizedMinLimitCycle

= = = =

CR IP T

where T˜ = Ta is the normalized period of the limit cycle. For numerical purposes, these quantities can be estimated by using a large value of time t. For the parameters displayed in (59), we solved the ODE for 0 ≤ τ ≤ 1000 and obtained the following values: 166.6 1.204 26.99 0.1842.

(67)

AN US

For the parameters displayed in (60), we solved the ODE for 0 ≤ τ ≤ 2000 and obtained the following values: T˜ = 59.21 NormalizedMeanRatio = 0.6585 (68) NormalizedMaxLimitCycle = 3.003 NormalizedMinLimitCycle = 0.0830.

2.4

ED

M

As we can see, if the condition for the emergence of a inverted biomass pyramid are satisfied, the mean value of time spent in one period of the steady limit cycle is spent as an inverted biomass pyramid. On the other, if the condition for the emergence of a inverted biomass pyramid are not satisfied, the mean value of time spent in one period of the steady limit cycle is spent as a regular biomass pyramid. In the following section, we will investigate the sensitivity of the input variables on the model.

Sensitivity Analysis

PT

In this section we will perform some local sensitivity analysis on important parameters discussed before. Suppose we are given a general ODE model of the form: dy = F(y, p, t), dt

y(0) = y0 ,

(69)

AC

CE

where p are parameters of the model. We wish to measure how much the solutions yi change when we perturb the parameters p. In order to do so, we will investigate the relative sensitivity matrix defined by:    ∂yi pj ˜ S = {˜ si,j } = . (70) ∂pj yi

When we cannot compute the derivatives explicitly, we can approximate them using a central difference like so: ∂yi yi (pj + ∆pj ) − yi (pj − ∆pj ) = + O((∆pj )2 ). (71) ∂pj 2∆pj In this paper, we will analyse the sensibility of the non-dimensionalized parameters {ei }5i=1 and c on the ¯ internal steady states defined in equation (57), the biomass ratio X¯ rY+X¯ , the period of the stable limit cycle T , the measure value A defined in equation (61) and the maximum and minimum values of the biomass ratio defined in equations (62) and (63), respectively. We will be performing this local sensitivity analysis for the values presented in equation (59). When computing the relative sensitivity matrix we use the exact derivative for the internal steady states and the biomass ratio since exact expression of

18

ACCEPTED MANUSCRIPT

these quantities are available. For the remaining with the step sizes: ∆e1 ∆e2 ∆e3 ∆e4 ∆e5 ∆c

quantities, we approximate the derivatives using (71) = = = = = =

0.1 0.1 0.1 0.1 0.001 0.01.

(72)

 Xr?    X?    Y?    Y?  ? +X ?  Xr    T˜    NormalizedMeanRatio     NormalizedMaxLimitCycle  NormalizedMinLimitCycle

e1 −1.1938

e2 0.97781

e3 −1.1589

e4 0.97781

e5 1.1589

0

0

−1.1852

1

1.1852

−1.0754

0.87693

−1.0392

0.87694

0.039344

−0.72548

0.59032

0.13815

−0.11657

−1.1381

1.7159

−1.3433

0.2179

−0.9906

−0.7336

11.4066

−8.4258

3.4773

−9.9574

−4.5730

17.7083 −6.2537

AN US



CR IP T

These step sizes were chosen in order for the condition (48) on the stability of the internal steady to remain unsatisfied. The relative sensitivity matrix we obtain is presented bellow in the matrix (73).

−15.1882

5.4408

−18.3120

−7.4648

4.8763

−1.3664

4.1360

0.3211

c −1.1589



   −1.1852    −0.039350      1.1381     0.2087    6.0462     9.3164  −0.4609

(73)

From the relative sensitivity matrix in (73), we can deduce this local sensitivity analysis for these parameters. most positively correlated e5 e5 e2 , e4 c e1 e1 e1 e2

M

Xr? X? Y?

Y? Xr? +X ?

(74)

PT

ED

T˜ NormalizedMeanRatio NormalizedMaxLimitCycle NormalizedMinLimitCycle

most negatively correlated e1 e3 , c e1 e5 e2 e4 e4 e1

3

CE

Sensitivity indices are positive or negative which indicates the nature of the relationship, and it is the magnitude that ranks the strength of the relationship as compared to the other parameters.

Impact of Fishing

AC

Overfishing depletes various segments of the marine food chain and is the biggest threat to coral reefs. We add fishing to our model and show that fishing pressure will change the inverted biomass pyramid. Also, we observe that only prey fishing will destroy the inverted biomass pyramid. In this model, we assume that the predator fishing rate is proportional to the predator biomass and prey fishing is similar to predator hunting. We test different forms of prey fishing and show that our results are qualitatively robust. The model equations incorporating fishing are provided by dxr dt dx dt dy dt

  xr = axr 1 − − αxr + β(r)x, K(r) = αxr − β(r)x − f (x)y − mf (x), = cf (x)y − dy − ly.

19

(75)

ACCEPTED MANUSCRIPT

where m is a prey fishing effort (per day) and l is the predator fishing effort (per day). We assume that the functional responses of prey to the predator or the human fishing are the same.

Y 1 0.8

CR IP T

0.6 0.4 0.2

0.6

0.8

1

1.2 r

1.4

1.6

1.8

2

AN US

0 0.4

Figure 9: The Hopf bifurcation diagram for predator vs refuge size.

X

1.4 1.2

M

1

0.6 0.4

ED

0

0.8

0.2 0

5e-05 0.0001 0.00015 0.0002 0.00025 l

(a) Predator vs predator fishing effect.

0.8

AC

0.6 0.4 0.2 0

0

5e-05

0.0001 l

0.00015

0.0002

1.4 Population density

CE

X_R 1

0

(b) Outside refuge Prey vs predator fishing effect.

PT

Y 0.5 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05 0

Prey inside refuge prey outside refuge Predator

1.2 1 0.8 0.6 0.4 0.2

5e-05 0.0001 0.00015 0.0002 0.00025 l

0 0

(c) Inside refuge prey vs predator fishing effect.

1

2

Time

3

4

5

5

x 10

(d) Population density vs time

Figure 10: The bifurcation diagram and the population density for the system (75) at r = 1. Parameters: a = 0.007, K = 0.8, α = 0.035, β = 0.0119, γ = 0.0112, ∆ = 1, d = 0.00007, c = 0.04, l = 0.00005, m = 0.0005.

20

ACCEPTED MANUSCRIPT

Figure 10 illustrates the bifurcation diagram for different predator fishing efforts and the total density of predator prey population over time. We see that increasing fishing efforts and decreasing refuge size stabilizes the system (see Figures 9, 10). Therefore, the predator-prey fishing efforts and the refuge size determine the stability of the system. From Figures 5 and 9, we conclude that we need larger refuge size to reach the limit cycle in the presence of fishing. Hence, increasing fishing effort eliminates the limit cycle. In this case, more predators and the external prey will be eliminated from the corel reef which leads to extiction of predators.

3.1

Mathematical Analysis

CR IP T

In this section, we study the positivity of solutions, the existence and local stability of all steady states. For simplicity, we assume f (x) = γx. The following theorem shows the positivity. Theorem 3.1. All solutions of the system (75) which start in R3+ = [0, ∞)3 will remain nonnegative. Proof. Similar to the proof of Theorem 2.1.

AN US

Model (75) has only three nonnegative equilibria, namely E0 (0, 0, 0), E1 (˜ xr , x ˜, 0) and E∗ (¯ xr , x ¯, y¯). The equilibrium E0 always exists and we show the existence of E1 and E∗ as follows: The steady states for this model are given by

E0 (0, 0, 0) = (0, 0, 0), E1 (˜ xr , x ˜, 0) = (

αK(r) K(r) ((a − α)mγ + aβ(r)), ((a − α)mγ + aβ(r)) , 0), a(β(r) + mγ) a(β(r) + mγ)2

(76)

M

as well as one internal steady state denoted by E∗ (¯ xr (r, l), x ¯(r, l), y¯(r, l)) where: r     K(r) β(r) (d+l) 2 x ¯r (r, l) = (a − α) + (a − α) + 4a K(r) , cγ 2a (77)

=

d+l cγ ,

y¯(r, l)

=

ca¯ xr (β(r) + mγ) − . (d + l) γ

ED

x ¯(r, l)

PT

Note that the equilibrium E1 is feasible if and only if (a − α)mγ + aβ(r) > 0,

AC

CE

and the equilibrium E∗ is feasible if and only if x ¯r >

(d + l)(β(r) + mγ) . caγ

Let the parameter l be fixed. We can observe that x ¯ is independent of the parameter m, whereas the value of y¯ depends on m. Therefore m has no effect on the interior equilibrium level on prey outside refuge. To check the effect of m on predator, differentiating y¯ with respect to m, we get d¯ y < 0. dm

This shows that y¯ is a strictly decreasing function of m. Therefore, as the prey fishing effort m increases, the prey species remain unchanged but the predator population decreases. This happens due to the loss of food for the predator species. Moreover, the predator population goes to extinction when m is large. Now let m be fixed, we can show that d¯ x > 0 and dl 21

d¯ xr > 0. dl

ACCEPTED MANUSCRIPT

Therefore x ¯r and x ¯ are strictly increasing functions of l. Also, s " # ca β(r)K(r) d¯ y (d + l) = −x ¯r . dl (d + l)2 acγ d¯ x dl

> 0, then d¯ y > 0 if d + l > x ¯r dl d¯ y < 0 if d + l < x ¯r dl

r r

acγ , β(r)K(r) acγ . β(r)K(r)

CR IP T

Since

AN US

Thus as the predator fishing effort l increases, the predator population increases when l is smaller than the threshold value. But if the fishing effort l gradually increases above the threshold value, that is, as the fishing effort l becomes large enough, then the increasing amount of the fishing effort l can decrease predator population. The Jacobian of the system (75) is given by     2a x β(r) 0 (a − α) − r   K(r) . (78) J(xr , x, y) =    α −β(r) − γy − mγ −γx 0 cγy cγx − d − l Theorem 3.2. The extinction steady state E0 is an unstable node if a − α > 0.

(79)

M

Proof. Evaluating the Jacobian at the point E0 , we obtain   a−α β(r) 0 . −(β(r) + mγ) 0 J(0, 0, 0) =  α 0 0 −(d + l) The eigenvalues are −(d + l), and the roots of the polynomial

ED

p(λ) = λ2 + (β (r) + mγ − (a + α)) λ − (aβ (r) + (a − α)mγ).

(80)

PT

The roots of a polynomial of order two have negative real parts if and only if its coefficients are positive. (81)

β(r) + mγ − (a − α) > 0,

(82)

−(aβ (r) + (a − α)mγ) > 0,

CE

Since the second equation can never be satisfied, there exists an eigenvalue of equation (75) with positive real part which implies E0 is an unstable node.

AC

Theorem 3.3. The predator free steady state E1 = (˜ xr , x ˜, 0) is locally asymptotically stable for αK(r) 2˜ xr d+l − 1 < 0 and > x ˜ , where x ˜ = ((a − α)mγ + aβ(r)). K(r) cγ a(β(r)+mγ)2 Proof. Evaluating the Jacobian at the point E1 = (˜ xr , x ˜, 0), we obtain   2a −(a − α) − K(r) x ˜r β(r) 0 . J (˜ xr , x ˜, 0) =  α −(β(r) + mγ) −γ x ˜ 0 0 cγ x ˜ − (d + l)

The eigenvalues are −(d + l) + cγ x ˜ < 0, and the roots of the characteristic polynomial      2˜ xr 2˜ xr 2 p(λ) = λ + β(r) + mγ + α + a − 1 λ + a(β(r) + mγ) − 1 + mαγ. K(r) K(r)

(83)

The roots of a polynomial of order two have negative real parts if and only if its coefficients are positive. 2˜ xr In our case, both coefficients are positive if and only if ( K(r) − 1) > 0. Therefore, the predator-free

2˜ xr equilibrium is locally asymptotically stable for ( K(r) − 1) > 0 and

22

d+l cγ

>x ˜.

ACCEPTED MANUSCRIPT

To perform stability analysis for the inner equilibrium, we require Li and Muldowney’s approach [22] and McCluskey and van den Driessche’s lemma [1, 23] for second and third degree polynomials. Lemma 3.4. Let A by an 3 × 3 matrix with real entries. If tr(A), detA and detA[2] are all negative, then all of the eigenvalues of A have negative real part. Theorem 3.5. The internal equilibrium steady state (xr , x, y) = (¯ xr , x ¯, y¯) exists if and only if x ¯r > (d + l)(β(r) + mγ) . E∗ is locally asymptotically stable if AB > αβ(r), and (d + l) − cγ x ¯>0 ca 2a¯ xr where A = K(r) − a + α > 0, B = β(r) + γ y¯ + mγ > 0.

implies



and

2a¯ xr K(r)

 −(A + B) −γ x ¯ 0 , cγ y¯ −(A + P ) β(r) J [2] (¯ xr , x ¯, y¯) =  0 α −(P + B)

(84)

(85)

(86)

− a + α > 0, B = β(r) + γ y¯ + mγ > 0, and P = (d + l) − cγ x ¯ > 0.

M

where A =



 β(r) 0 −B −γ x ¯ , cγ y¯ −P

AN US

−A J(¯ xr , x ¯, y¯) =  α 0

CR IP T

Proof. Consider the internal equilibrium E∗ = (¯ xr , x ¯, y¯).   2a¯ xr − a + α] β(r) 0 −[ K(r) , J(¯ xr , x ¯, y¯) =  α −(β(r) + γ y¯ + mγ) −γ x ¯ 0 cγ y¯ −[(d + l) − cγ x ¯]

ED

We can calculate the trace and determinant of J and the determinant of J [2] to determine the stability of E∗ : tr(J(E∗ )) = −(A + B + P ) < 0,

with

PT

  det(J(E∗ )) = − Acγ 2 x ¯y¯ + P (AB − β(r)α) ,

and

 det(J [2] (E∗ )) = − (A + B)[(A + P )(B + P ) − αβ(r)] + cγ 2 x ¯y¯(B + P ) .

Predator-Prey Biomass Ratio

AC

3.2

CE

Here, tr(J), detJ and detJ [2] are all negative if AB − αβ(r) > 0, then all of the eigenvalues of J have negative real parts.

In this section, we provide sufficient conditions for the emergence of inverted biomass pyramids of the system (75). Also, we study the influence of fishing efforts on an inverted biomass pyramid. The predator-prey biomass ratio at the inner equilibrium point is given by y¯(r, l) cα(cγ − d − l)¯ xr (r, l) − (d + l)cβ∆ − m(d + l)(cγ − d − l) = , x ¯(r, l) + x¯r (r, l) (d + l)[(cγ − d − l)¯ xr (r, l) + (d + l)∆]

23

(87)

ACCEPTED MANUSCRIPT

where =

x ¯r (x, l)

=

x ¯(x, l)

=

y¯(x, l)

=

γx , ∆+x s " #   K(r) β(r) ∆β(r)(d + l) 2 (a − α) + (a − α) + 4a , 2a K(r) cγ − d − l ∆(d + l) , cγ − d − l cα c∆β(r) x ¯r (x, l) − − m, d+l cγ − d − l

with cγ > d + l.

CR IP T

f (x)

AN US

0

ED

M

10

0

1

Time

2

3

5

x 10

Figure 11: The ratio of the population densities over time.

AC

CE

PT

Predator−prey Biomass Ratio

Figure 11 shows the predator-prey biomass ratio over the time. Additionally, In Figure 12, we display the predator-prey biomass ratio for various refuge sizes and fishing rates. We study the effect of fishing on the predator-prey biomass ratio by comparing equations (87) and (35) and the Figure 12. We observe that the predator-prey biomass ratio is a decreasing function of the fishing effort and the biomass ratio is less than one in the presence of fishing pressure. In the presence of predator fishing, the effect of prey fishing is not very strong.

24

2 1.5

m=0 m=0.07 m=0.08

1 0.5 1

1.2

1.4 1.6 Refuge Size

1.8

0.16

0.12

m=0 m=0.020 m=0.026

0.1 0.08 0.06 1

1.2

AN US

Predator−prey Biomass Ratio

(a) l=0

0.14

2

CR IP T

Predator−prey Biomass Ratio

ACCEPTED MANUSCRIPT

1.4 1.6 Refuge Size

1.8

2

(b) l=0.0002

M

Figure 12: Predator-prey biomass ratio at an internal equilibrium point as a function of refuge size with different prey fishing rates (m) and two predator rates (l).

PT

ED

Our results are independent of the form of prey fishing. let p(x) be the general prey fishing rate. The modified equations are   dxr xr − αxr + β(r)x, = axr 1 − dt K(r) dx γx (88) = αxr − β(r)x − ∆+x y − p(x), dt dy γx = c ∆+x y − dy − ly. dt

CE

The predator-prey biomass ratio at the internal equilibrium point is

AC

  y¯(r, l) c α(cγ − d − l)x¯r (r, l) − β(r)∆(d + l) − (cγ − d − l)p(¯ x(r, l)) . = x ¯(r, l) + x¯r (r, l) d+l (cγ − d − l)x¯r (r, l) + ∆(d + l)

y¯(r,l) The biomass ratio at the fished reef ( x¯(r,l)+ x¯r (r,l) ) is less than the biomass ratio at a reef without y¯ fishing ( x¯+¯ ) xr   y¯(r, l) c α¯ xr − β(r)∆d y¯ ≤ = . x ¯(r, l) + x¯r (r, l) d x ¯r + ∆d x ¯+x ¯r

As a result of fishing, the predator-prey biomass ratio is less than the biomass ratio at reefs without fishing. This result is robust under different forms of prey fishing. Moreover, high fishing pressure will destroy the inverted biomass pyramid. We now study the sufficient conditions for the existence of an inverted biomass pyramid.

Theorem 3.6. Given the model (75) with the proper conditions satisfied for the stability of the internal steady state shown in (77), there will be an inverted biomass pyramid if the following conditions are

25

ACCEPTED MANUSCRIPT

satisfied: d+l <α c

and

q

where p(x) = vx.

2a (c(β−v)+(d+l))∆(d+l) K(r) (cα−d−l)(cγ−d−l) β(r) (d+l)∆ (a − α)2 + 4a K(r) (cγ−d−l)

< 1,

(89)

Rearranging and obtain

(r,l) − β(r) − v α x¯x¯r(r,l) d+l < , x ¯ (r,l) r c +1 x ¯(r,l)

α + β(r) − v d+l < α − x¯ (r,l) , r c +1 x ¯(r,l)

x ¯r (r,l) x ¯(r,l)

<α−

d+l cα − d − l = , c c

AN US

α + β(r) − v

cα − d − l > c(α − v + β(r)) Assuming cα − d − l > 0, we obtain

1

x ¯r (r,l) x ¯(r,l)

.

x ¯r (r, l) cβ(r) + (d + l) − cv > . x ¯(r, l) cα − d − l

M

Simplifying, we obtain s

CR IP T

Proof. For the existence of an inverted biomass pyramid, we require that   y¯(r, l) c α¯ xr (r, l) − β(r)¯ x(r, l) − v¯ x(r, l) 1< = . x ¯(r, l) + x¯r (r, l) d+l x ¯(r, l) + x ¯r (r, l)

β(r) (d + l)∆ 2a (c(β − v) + (d + l))∆(d + l) > , K(r) (cγ − d − l) K(r) (cα − d − l)(cγ − d − l)

PT

ED

(a − α)2 + 4a

q

2a (c(β−v)+(d+l))∆(d+l) K(r) (cα−d−l)(cγ−d−l)

β(r) (d+l)∆ (a − α)2 + 4a K(r) (cγ−d−l)

< 1.

AC

CE

This theorem shows how the appearance of an inverted biomass pyramid depends on fishing. The first condition shows that the fishing effect plays the same role as the predator death rate, and hence fishing can destroy an inverted biomass pyramid. The second condition shows a very nonlinear and complicated dependence of an inverted biomass pyramid on the fishing effect.

4

Impact of Predator Migration

Migration is a common animal behavior and occurs in all kinds of habitats including coral reef. It may be defined as a synchronized fish movement between two habitats. In this section, we discuss the impact of predator migratory behavior on the population dynamics and the biomass ratio of a predator-prey system. The model with migrations of the predator is given by   dxr xr − αxr + β(r)x, = axr 1 − dt K(r) dx (90) = αxr − β(r)x − f (x)y, dt dy = cf (x)y − dy + I(x) − Ey, dt 26

ACCEPTED MANUSCRIPT

CR IP T

where f (x) = γx, the immigration rate I(x) = ζx is an increasing function of x, while E is the assumed constant per capita emigration rate (when the predator is full, they leave the studied coral reef region). Figure 13 illustrates the predator-prey population density for different predator migration rates and Figure 14 shows the bifurcation diagram for the system (90) with refuge size (r) as a bifurcation parameter. Also, we plot the bifurcation diagram for the system (90) with emigration rate (E) as a bifurcation parameter (see Figure 15). We observe that increasing emigration rates and decreasing refuge size stabilize the system. From Figures 5 and 14, we conclude that we need larger refuge size to reach the limit cycle in the presence of predator emigration. Hence, increasing predator emigration eliminates the limit cycle. Our results are qualitatively the same as the effect of fishing but quantitatively different. Prey inside refuge prey outside refuge Predator

1.2 1 0.8 0.6 0.4 0.2 0 0

1

AN US

Population density

1.4

2

Time

3

4

5

5

x 10

(a) E=0.00005, ζ = 0.00001

1 0.8 0.6

M

1.2

ED

Population density

1.4

0.4 0.2

Population density

CE

PT

0 0

AC

Prey inside refuge prey outside refuge Predator

2

4

Time

6

8

10 4 x 10

(b) E=0.0005, ζ = 0.0001 1.4 1.2 1

Prey inside refuge prey outside refuge Predator

0.8 0.6 0.4 0.2 0 0

2

4

Time

6

8

10 4 x 10

(c) E=0.005, ζ = 0.001

Figure 13: Density of predator-prey population for predator migration. Parameters: a = 0.0007, K = 0.8, α = 0.035, β(r) = 0.0119, γ = 0.0112, ∆ = 1, d = 0.00007, r = 1, c = 0.04.

27

ACCEPTED MANUSCRIPT

Y 1 0.8 0.6

CR IP T

0.4 0.2 0 0.4

0.6

0.8

1

1.2 r

1.4

1.6

1.8

2

AN US

Figure 14: The Hopf bifurcation diagram for predator vs refuge size.

X

Y

0.5

0.5

0.45

0.45

0.4

0.4

0.35

0.35

0.3

0.3

0.25 0.2

M

0.25 0.2 0.15 0.1

0.15 0.1

0.05 0 0

1e-05

2e-05

3e-05

ED

0.05

4e-05

5e-05 r

6e-05

7e-05

0 0

8e-05

9e-05

PT

(a) Predator vs predator emigration rate. 0.5 0.45

0.35 0.3

AC

0.25

3e-05

4e-05

5e-05 r

6e-05

7e-05

8e-05

9e-05

0.0001

1.4

CE

0.4

2e-05

(b) Outside refuge prey vs predator emigration rate.

Population density

X_R

1e-05

0.0001

0.2

0.15

0.1

Prey inside refuge prey outside refuge Predator

1.2 1 0.8 0.6 0.4 0.2

0.05

0 0

0 0

1e-05

2e-05

3e-05

4e-05

5e-05 r

6e-05

7e-05

8e-05

9e-05

0.0001

(c) Inside refuge prey vs predator emigration rate.

1

2

Time

3

4

5

5

x 10

(d) Population density vs time.

Figure 15: The bifurcation diagram and the population density of the system (90) with r = 1 and ζ = 0.00001. Parameters: a = 0.007, K = 0.8, α = 0.035, β = 0.0119, γ = 0.0112, ∆ = 1, d = 0.00007, r = 0.8, c = 0.04, E = 0.00005.

28

ACCEPTED MANUSCRIPT

4.1

Mathematical Analysis

In this section, we study the positivity of solutions, the existence and local stability of all steady states. Theorem 4.1. All solutions of the system (90) which start in R3+ = [0, ∞)3 will remain nonnegative. Proof. Similar to the proof of Theorem 2.1.

E0 (0, 0, 0) = (0, 0, 0), as well as one internal steady state denoted by E∗ (x¯r , x ¯, y¯) = (

CR IP T

It can be checked that model (90) only has two nonnegative equilibria, namely E0 (0, 0, 0) and E∗ (¯ xr , x ¯, y¯). The equilibrium E0 exists and we shall show the existence of E∗ as follows: The steady states for this model are given by

γζ ζ β(r) + x ¯2 , x ¯, x ¯), α α(d + E − cγ x ¯) d + E − cγ x ¯

AN US

and x ¯ satifies the equation

f (¯ x) = A¯ x3 + B x ¯2 + C x ¯ + D = 0,

(91)

(92)

where A = aγ 2 (ζ − cβ(r))2 > 0, B = −2a(d + E)γ(c − ζβ(r)) − K(r)α2 c2 γ 2 β(r)2 , C = aβ(r)2 (d + E)2 + aK(r)α2 (d + E)cγβ(r)2 − K(r)αγ(a − α),

M

and D = K(r)αβ(r)(d + E)[α(1 − (d + E)β(r) − a)].

Note that equation (92) has a unique positive solution x = x ¯ if the following inequalities hold:

ED

c > ζβ(r)

PT

β(r)2 (d + E)2 + aK(r)α2 (d + E)cγβ(r)2 < K(r)αγ(a − α) α(1 − (d + E)β(r) < a

Knowing the value of x ¯, the value of y¯ can be computed. Note that for y¯ to be positive, we must have

CE

d + E > cγ x ¯

AC

The Jacobian of the system (90) is given by    2xr − 1) + α − a(  K(r) J(xr , x, y) =   α 0

β(r)

0

−(β(r) + γy) −γx cγy + ζ −(d + E − cγx)



 . 

(93)

Theorem 4.2. The total extinction equilibrium steady state E0 = (0, 0, 0) is an unstable node if 1 − c > 0. Proof. Evaluating the Jacobian at the point (0, 0, 0), we obtain   a − α β(r) 0 . −β(r) 0 J(0, 0, 0) =  α 0 ζ −(d + E)

(94)

The characteristic polynomial of this system is given by

p(λ) = (λ + d + E) λ2 + (cβ (r) − a + α) λ − (cβ (r) (a − α) + β(r)α). 29




(95)

ACCEPTED MANUSCRIPT

The eigenvalues are −(d + E), and the roots of the polynomial λ2 + (cβ(r) − a + α) λ − (cβ(r)(a − α) + β(r)α)

(96)

will have roots with negative real parts if and only if (97)

cβ(r) − (a − α) > 0,

(98)

−(c + (1 − c)α) > 0.

CR IP T

Since the second equation can never be satisfied because (1 − c) > 0. There exist an eigenvalue of equation (90) with positive real part which implies that E0 is an unstable node. Theorem 4.3. The internal steady state E∗ = (¯ xr , x ¯, y¯) is positive if and only if c > ζβ(r), β(r)2 (d + 2 2 2 2 E) + aK(r)α (d + E) cγβ(r) < K(r)αγ(a − α), α(1 ¯. E∗ is locally  − d − E)β(r) < a and d + E > cγ x 2¯ xr asymptotically stable if J1 J3 > αβ(r), where J1 = a( − 1) + α) > 0, J2 = (β(r) + γ y¯) > 0, K(r) J3 = (d + E − cγ x ¯) > 0 and J4 = cγ y¯ + ζ > 0.

AN US

Proof. Consider the internal equilibrium E∗ = (¯ xr , x ¯, y¯).    2¯ xr − a( − 1) + α) β(r) 0  K(r) J(¯ xr , x ¯, 0) =   α −(β(r) + γ y¯) −γ x ¯ 0 cγ y¯ + ζ −(d + E − cγ x ¯) that is,

and

  

(99)



a(

2¯ xr K(r)

 −J1 β(r) 0 −J2 −γ x ¯ , J(¯ xr , x ¯, 0) =  α (100) 0 J4 −J3  − 1) + α) > 0, J2 = (β(r) + γ y¯), J3 = (d + E − cγ x ¯) > 0 and J4 = cγ y¯ + ζ,

M



 −(J1 + J2 ) −γ x ¯ 0 . ζ −(J1 + J3 ) β(r) J [2] (E∗ ) =  0 α −(J2 + J3 ) 

ED

where J1 =



(101)

CE

PT

We can compute the trace and determinant of J and the determinant of J [2] evaluated at E∗ to determine the stability of E∗ : with

AC

and

tr(J(E∗ )) = −2(J1 + J2 + J3 ),

det(J(E∗ )) = − [J1 J2 γ x ¯ + J2 (J1 J3 − αβ(r))] ,

det(J [2] (E∗ )) = − [(J1 + J3 ) ((J1 + J2 )(J2 + J3 ) − αβ(r)) + γJ4 x ¯(J2 + J3 )] .

Here, tr(J(E∗ )), det(J(E∗ )) and det(J [2] (E∗ )) are all negative if J1 J3 > αβ(r). In this case, the eigenvalues of J have negative real parts.

4.2

Predator-Prey Biomass Ratio

In this section, we provide sufficient conditions for the emergence of inverted biomass pyramids of the system (90). Also, we study the influence of predator migration on an inverted biomass pyramid. The predator-prey biomass ratio at the internal equilibrium point is given by y¯ ζα = . x ¯+x ¯r (d + E − cγ x ¯)(α + β(r)) + γζ x ¯ 30

(102)

ACCEPTED MANUSCRIPT

0

10

0

1

2

Time

0

10

3

5

0

1

x 10

(a) E=0.00005, ζ = 0.00001. 1

0

10

0

PT

Time

2

3

5

x 10

(b) E=0.00005, ζ = 0.000001.

M

10

ED

Predator−prey Biomass Ratio

1

10

CR IP T

Predator−prey Biomass Ratio

1

10

AN US

Predator−prey Biomass Ratio

We plot the predator-prey biomass ratio over time for different predator migration rates in Figure 16. We study the effect of predator migration on the predator-prey biomass ratio by comparing expression (102) and Figure 16. We observe that the predator-prey biomass ratio is a decreasing function of predator emigation and increasing function of predator immigration. When we slightly increase the emigration rate we see a large decrease in the biomass ratio (see 16(a), (c)) but when we decrease the immigration rate sufficiently, we see a slight decrease in biomass ratio (see 16(a), (b)). Thus, we conclude that the predator emigration rate, E, plays an important role in the predator-prey biomass ratio.

1

Time

2

3

5

x 10

(c) E=0.00007, ζ = 0.00001.

CE

Figure 16: Predator-prey biomass ratio vs emigration rate and immigration rate. Parameters: a = 0.0007, K = 0.8, α = 0.035, β(r) = 0.0119, γ = 0.0112, ∆ = 1, d = 0.00007, r = 1, c = 0.04. In the following, we obtain the sufficient conditions for the existence of an inverted biomass pyramid.

AC

Theorem 4.4. Given the model (90) with the proper conditions satisfied for the stability of the positive internal steady state shown in (91), an inverted biomass pyramid occurs if the following conditions are satisfied: (d + E)(α + β(r)) − ζα E+d >x ¯ and x ¯> . (103) cγ γ[(α + β(r))c − ζ] Proof. For the existence of biomass pyramid, we require that 1< Rearranging and obtain

y¯ ζα = . x ¯r + x ¯ (d + E − cγ x ¯)(α + β(r)) + γζ x ¯ 

 (d + E − cγ x ¯)(α + β(r)) + γζ x ¯ < 1, ζα 31

(104)

ACCEPTED MANUSCRIPT

x ¯>

(d + E)(α + β(r)) − ζα . γ[(α + β(r))c − ζ]

From the expression (104) in the proof, we can observe that increasing the predator emigration rate always decreases the biomass ratio but decreasing the predator immigration can either increase or decrease the biomass ratio.

Discussion

CR IP T

5

AC

CE

PT

ED

M

AN US

In this paper, we develop a model describing the predator-prey interactions when prey have access to a refuge. We investigate the effects of the prey refuge on the predator-prey population dynamics and the biomass ratio. Furthermore, we study the predator-prey model with fishing efforts and predator migrations. The model can be applied to many examples such as prey fish hidden in coral reefs, crabs hidden under rocks, and lemmings hidden underground. Of course, for a specific example, our model may need to be modified slightly. Arguably the most interesting aspect of this model occurs when the internal equilibrium state is unstable. We derive the conditions for the existence and stability of steady states together with some critical values of the parameters at which the system undergoes a Hopf bifurcation around the internal equilibrium. Specifically, this periodic solution can only bifurcate from the internal equilibrium point. The existence of a Hopf bifurcation determines the case of the oscillatory coexistence of predator and prey in the presence of a refuge. In addition, increasing predator emigration or predator fishing eliminates the limit cycle. Previous studies showed that the refuge used by prey have a stabilizing effect on simple models. The influence of refuge on the predator functional response and on prey productivity was studied in [29]. The results of this paper suggest that the relationship between the refuge size and population dynamics may be quite complex where the amount of protected prey biomass is a nonlinear function of the total prey biomass and the available refuge size. We show that the use of a refuge for prey has a strong influence on the predator-prey population dynamics. We analyze the dynamics of the biomass ratio. We obtain the sufficient conditions for the existence of an inverted biomass pyramid. We observe that fairly high fishing pressure or strong predator migration will change the inverted biomass pyramid into a regular biomass pyramid. This is consistent with previous studies [27]. As seen in Figure 4, there exists a stable limit cycle. Populations of predator and prey can cycle between an inverted biomass pyramid and a regular biomass pyramid throughout the limit cycle. For the majority of ecological systems, there will be less biomass at each tropic level resulting in a regular biomass pyramid. Our model suggests that inverted biomass pyramids are possible in the presence of a stable limit cycle. In the presence of a stable limit cycle, we calculate the mean value, the maximum and minimum values of the biomass ratio over one period of the limit cycle. We perform sensitivity analysis on a normalized version of our model to investigate the effect of our normalized parameters on the studied quantities. Also, we show that this biomass ratio decreases when the fishing effort takes place. Our analytical and numerical results indicate that the dynamic behaviors of the model not only depend on the prey refuge size r but also on the fishing efforts l and m as well as the predator migration rates E and ζ. Hence, our study indicates that the refuge plays an important role to guarantee the coexistence of the predator and prey. In order to show the global stability, we need to find a Lyapunov function since our system is nonmonotone three-dimensional. This is an open mathematical question. The next step is to incorporate the maturation delay of the predator into our model. Introducing delay will force limit cycles more likely to occur. In addition, we can expect the phase shift in the predator-prey cycle [25].

32

ACCEPTED MANUSCRIPT

Acknowledgements We would like thank the postdoc Chunhua Shan for his help on XppAut. We also would like to thank the editor and the referees for their insightful comments. The first author’s research is partially supported by NSERC.

References

CR IP T

[1] Arino, J., McCluskey, C.C., van den Driessche, P., (2003). Global results for an epidemic model with vaccination that exhibits backward bifurcation, SIAM J. APPL. MATH. 64, 260-276. [2] Berger, J., (1991). Pregnancy incentives, predation constraints and habitat shifts: experimental and field evidence for wild bighorn sheep, Animal Behaviour. 41, 61-77. [3] Berezovskaya, F.S., Song, B., and Castillo-Chavez, C., (2010). Role of prey dispersal and refuges on predator-prey dynamics, SIAM J. APPL. MATH. 70, 1821-1839.

AN US

[4] Clarke, M.F., da Silva, K.B., Lair, H., Pocklington, R., Kramer, D.L., Mclaughlin, R.L., (1993). Site familiarity affects escape behaviour of the eastern chipmunk, Tamius striatus, Oikos. 66, 533-537. [5] Collings, J.B., (1995). Bifurcation and stability analysis of a temperature-dependent mite predatorprey interaction model incorporating a prey refuge, Bulletin of Mathematical Biology. 57, 63 - 76. [6] Dill, L.M., Houtman, R., (1989). The influence of distance to refuge on flight-initiation distance in the grey squirrel (Sciurus carolinensis), Canadian Journal of Zoology. 67, 232-235. [7] Dubey, B., (2007). A Prey-Predator Model with a Reserved Area, Nonlinear Analysis: Modelling and Control. 12, 479-494.

M

[8] Ellner, S.P., McCauley, E., Kendall, B.E., Briggs, C.J., Hosseini, P.R., (2001). Habitat structure and population persistence in an experimental community. 412, 538 - 543.

ED

[9] Lotka, A.J., (1925). Elements of Physical Biology, Williams and Wilkins Company, Baltimore. [10] Freedman, H.I. and Waltman, P., (1984). Persistence in models of three interacting predator-prey population, Math Biosci. 68, 213 - 231.

PT

[11] Gonzalez-Olivars, E. and Ramos-Jiliberto, R., (2003). Dynamics consequences of prey refuges in a simple model system: more prey, few predators and enhanced stability, Ecol. Model. 166, 135.

CE

[12] Hassell, M. P., (1978). The dynamics of arthropod predator-prey systems, Princeton University Press, Princeton, N.J. [13] Hixon, M.A., Beets, J.P., (1993). Predation, Prey Refuges, and the Structure of Coral-Reef Fish Assemblages, Ecological Monographs. 63,77 - 101

AC

[14] Holling, C.S., (1959). Some characteristics of simple types of predation and parasitism. 91, 385.

[15] Holmes, W.G., (1991). Predator risk affects foraging pikas: observational and experimental evidence, Animal Behaviour. 42, 111-119. [16] Huffaker, C. B., (1958). Experimental studies on predation: dispersion factors and predator-prey oscilla- tions, J. Hilgardia. 27, 343 - 383. [17] Aoki, I., (2006). Ecological pyramid of dissipation function and entropy production in aquatic ecosystems, J. Ecological Complexity. 3, 104 - 108. [18] Anagnost, J.J. and Desoer, C.A., (1991). An elementary proof of Routh-Hurwitz stability criterion, J. Circuits, Systems and Signal Processing. 10, 101 - 114.

33

ACCEPTED MANUSCRIPT

[19] Ji, L., and Wu, C., (2010). Qualitative analysis of a predator-prey model with constant-rate prey harvesting incorporating a constant prey refuge, Nonlinear Analysis: Real World Applications. 11, 2285 - 2295. [20] Kar, T.K., (2006). Modelling and analysis of a harvested prey-predator system incorporating a prey refuge, Journal of Computational and Applied Mathematics. 185, 19 - 33. [21] Knowlton, N and Jackson., J, (2008). Shifting Baselines, Local impacts, and global change on coral reefs, PLoS Biol. 6, e54.

CR IP T

[22] Li, M.Y., Muldowney, J.S., (1995). Global stability for the SIER model in epidemiology, Math. Biosci. 125, 155-164. [23] McCluskey, C.C., van den Driessche, P., (2004). Global analysis of two tuberculosis models, Journal of Dynamics and Differential Equations 16, 139-166. [24] McNair, J. N., (1986). The effects of refuges on predator-prey interactions: a reconsideration, Theor. Popul. Biol. 29, 38 - 63.

AN US

[25] Li, M., Lin, X., and Wang, H., (2014) Global Hopf branches and multiple limit cycles in a delayed Lotka-Volterra predator-prey model, Discrete and Continuous Dynamical System Series B. 19, 747760. [26] Oaten, A. and Murdoch W.W., (1975). Functional response and stability in predator-prey systems, Amer. Nar. 109, 289 - 298. [27] Sandin, S., Smith, J., DeMartini, E., Dinsdale, E., Donner, S., Friedlander, A., Konotchick, T., Malay, M., Maragos, J., and Obura, D., (2008). Baselines and degradation of coral reefs in the nothern line islands, PLoS ONE. 3(2), e1548.

M

[28] Sakanoue, S., (2009). Resource-based approach to modelling the dynamics of interacting populations, J. Ecological Modelling. 220, 1383 - 1394.

ED

[29] Singh, A., Wang, H., Morrision, W.M. and Weiss, H., (2012). Modeling Fish Biomass Structure at Near Pristine Coral Reefs and Degradation by Fishing, J. Biological Systems. 20, 20 - 36. [30] Smith, M.J., (1974). Models in ecology, Cambridge:Cambridge University Press.

PT

[31] Trebilco, R., Baum, J.K., Salomon, A.K. and Dulvy, N.K., (2013). Ecosystem ecology: Size-based constraints on the pyramids of life, Trends in Ecology & Evolution. 28, 423 - 431.

CE

[32] Taylor, R.J., (1984). Predation, Chapman and Hall, New York.

AC

[33] Wang, H., Morrison, W., Singh, A. and Weiss, H., (2009). Modeling inverted biomass pyramids and refuges in ecosystems, Ecological Modelling. 220, 1376 - 1382.

34