Neimark-Sacker bifurcations in a host-parasitoid system with a host refuge

In this work, the effect of a host refuge on a host-parasitoid inter- 
action is investigated. The model is built upon a modi ed Nicholson-Bailey 
system by assuming that in each generation a constant proportion of the host 
is free from parasitism. We derive a sucient condition based on the model 
parameters for both populations to coexist. We prove that it is possible for 
the system to undergo a supercritical and then a subcritical Neimark-Sacker 
bifurcation or for the system only to exhibit a supercritical Neimark-Sacker 
bifurcation. It is illustrated numerically that a constant proportion of host 
refuge can stabilize the host-parasitoid interaction.


1.
Introduction. There are many mathematical models of predator-prey interactions in the literature investigating either predator's response to prey such as using different functional and numerical responses or predator's response to other predators such as considering mutual interference. The effects of prey responds to predators on the other hand have not received much attention. One of the mechanisms by which prey responses to predation is that of using spatial refuges. In many communities, prey populations often have access to areas where the prey are safe from their predators [8,22,23]. It is widely believed that prey refuges play important roles in natural communities and the spatial arrangement of these resources can have critical outcomes for individual and population fitness.
Traditionally, there are two common ways of incorporating spatial refuges into mathematical models [7,21]. Namely, either a fixed number of the prey is invulnerable to predation or a fixed fraction of the prey is free from predators. The existence of refuges clearly have important effects on the coexistence of predators and prey. In [16], Maynard Smith considered a simple Lotka-Volterra predator-prey model with these two types of refuges. It is well known that in the absence of refuges the model has neutrally stable cycles. Maynard Smith showed that a constant fraction of refuge does not change this qualitative behavior. However, a constant number of refuges does make the coexisting equilibrium to become asymptotically stable. Since then, many continuous-time predator-prey models with prey refuges have appeared in the literature. See [3,5,10,13,14,15,17,18,19,20] and references therein for a partial list of the continuous-time mathematical models with prey refuges.
Although there are many populations in nature that do not reproduce continuously, discrete-time predator-prey models with spatial refuges, on the other hand, have not been extensively investigated. The earliest research dates back to Hassell [7], where stability is performed on the Nicholson-Bailey host-parasitoid model with either a constant number of prey refuge or a constant proportionality of prey refuge. Hassell concluded that proportionate refuges can lead to stability.
In this work, we modify the model of Hassell [7] by incorporating densitydependence into the host population and assume that in each generation a constant proportion of the host population is free from parasitism. Specifically, the growth of the host population in the absence of the parasitoids follows the Beverton-Holt equation. We prove that the system can either undergo a supercritical and then a subcritical Neimark-Sacker bifurcation or only a supercritical Neimark-Sacker bifurcation appears. We also provide bifurcation diagrams using the constant proportion as the bifurcation parameter. It is concluded that the constant proportion of prey refuge can stabilize the host-parasitoid interaction.
In the following section we present the model and its analysis. We derive conditions for the coexistence of both populations. Stability of the unique interior steady state is also studied. Section 3 provides criteria for Neimark-Sacker bifurcations. Numerical examples using Matlab are presented to illustrate theoretical findings. The final section provides a brief summary and conclusions. 2. The model and analysis. Let x(t) and y(t) be the host and parasitoid populations at generation t = 0, 1, · · · , respectively. In [7], Hassell proposes the following model to investigate the impact of constant proportionate refuges upon the hostparasitoid interaction: where γ, 0 < γ < 1, is the constant fraction of the hosts that are available to the parasitoids in each generation. The parameter λ > 0 is the intrinsic growth rate of the hosts and a > 0 denotes searching efficiency of the parasitoids. When there is no spatial refuge, γ = 1, (1) becomes the classical Nicholson-Bailey model and it is well known that the coexisting steady state is always unstable [1].
It can be easily seen that if (1 − γ)λ > 1, then (1) and thus lim t→∞ x(t) = ∞. Therefore, solutions of (1) are not bounded if (1 − γ)λ > 1. On the other hand, (1) always has the extinction steady state (0, 0) which is globally asymptotically stable if λ < 1. Let λ > 1. A simple calculation shows that (1) has a unique interior steady stateÊ = (x,ŷ) if and only ) and x =ŷ γ(1 − e −aŷ ) . The Jacobian matrix of (1) atÊ is J(Ê) = 1 −λγax It is clear that |trJ(Ê)| < 1 + detJ(Ê) holds trivially and thereforeÊ will lose its stability only if detJ(Ê) > 1. We expect that (1) undergoes a Neimark-Sacker bifurcation when detJ(Ê) = 1. In [7, page 67], a stability region in terms of the parameters λ and 1 − γ is graphed and it is concluded that the host refuge can stabilize the host-parasitoid interaction. The author also indicates that there is still a large parameter space with which there is either no coexisting steady state or the steady state is unstable. The Beverton-Holt equation is one of the simplest discrete-time models. It is regarded as the discrete analogue of the well-known continuous-time logistic equation. It models contest competition between individuals within the same population when competition involves resources such as food or mates that are stable [2]. Since majority of the continuous-time systems use logistic equations to model prey's dynamics in the absence of the predator [5,10,13,15], it is natural to propose a parallel discrete-time system to study the impact of host refuge. Therefore, in this study we will focus on the effects of refuge when the host population is governed by the Beverton-Holt equation. Since in the absence of the parasitoids the host population is always stabilized, any destabilization will be due to the interaction with the parasitoids.
The proposed host-parasitoid model is given by the following system where λ, k, a and β are positive constants and 0 < γ < 1. The parameters λ, γ and a have the same biological interpretations as those in (1). The parameter β is the number of parasitoids produced by each parasitized host and k is related to the carrying capacity of the hosts. Observe that the prey population of both types, namely those refuged and those that are under predation pressure, has the same growth rate. This simple assumption is also assumed in many of the continuoustime models [3,5,7,10,13,14,15]. Furthermore, the host's carrying capacity in (2) is λ − 1 k . If k = 0, then resources of the host population become unlimited and we recover model (1). There are five parameters in model (2). Lettingŷ = ay,x = kx,β = βa/k, and ignoring the hats, (2) is converted into the following system with only three parameters: where 0 < γ < 1, λ > 0 and β > 0. Our goal is to study system (3). First, let us make some observations for the limiting cases of model (3). If γ = 0, then all of the hosts are safe within a refuge and the parasitoids go extinct since there are no hosts available for parasitism and the parasitoid is specialized to this particular host population. As a result, the host population grows according to the Beverton-Holt model and the population is stabilized. On the other extreme, γ = 1, then there is no spatial refuge for the hosts and the host-parasitoid interaction reduces to the following model Dynamics of (4) were discussed and compared with the corresponding model of Allee effects in [12]. In particular, both populations in (4) can coexist under some restrictions of the parameters, which contrast to the model of Allee effects. Using β as the bifurcation parameter, it is clear that system (4) undergoes a Neimark-Sacker bifurcation when the unique interior steady state loses its stability.
To study model (3), first observe that solutions of (3) remain nonnegative and are bounded. If there are no parasitoids initially, y(0) = 0, then y(t) = 0 for t > 0 and (3) reduces to the Beverton-Holt equation , and the host population will eventually be stabilized. For the full system (3), there always exists a trivial steady state E 0 = (0, 0) and its stability depends on the Jacobian matrix . It follows that E 0 is asymptotically stable if λ < 1 and is a saddle point if λ > 1. It can be easily shown that E 0 is globally asymptotically stable in R 2 + if λ < 1. In such case, the host population goes extinct prior to parasitism since its intrinsic growth rate is too small.
Let λ > 1. The y-component of an interior steady state (x, y) satisfies and the x-component is given by , y > 0 where lim y→∞ h(y) = ∞, lim y→0 + h(y) = 1 βγ , and Though strictly decreasing, g(x) may not be defined for all x ≥ 0. To discuss the existence of an interior steady state, we consider the following two cases separately, namely In case (a), g(x) is strictly decreasing for x ≥ 0 with Therefore, (5) and (6)  is defined and strictly decreasing for Thus, model (3) has an interior steady state if and only if βγx > 1 for case (b) as well, and the interior steady state is unique whenever it exists. We conclude that if λ > 1, then (3) has a unique interior steady state if and only if βγx > 1.
In the following approach, we will fix λ and let E * = (x * , y * ) move on the curve y = g(x) given in (5). Hence, x * = x * (β) and y * = y * (β). We first study the magnitude of E * with respect to β. Observe that By (6), x * satisfies g(x * ) = βγx * (1 − e −g(x * ) ). Differentiating both sides of the above equation with respect to β and simplifying, we get Similarly, e y * = γλ . Differentiating both sides with respect to β, we obtain from (10) that Consequently, x * is strictly decreasing and y * is strictly increasing in β.
It is expected that both host and parasitoid populations can coexist if βγx > 1. A useful mathematical tool to prove coexistence of the interacting populations is the concept of uniform persistence. System (3) is said to be uniformly persistent if there exists η > 0 such that lim inf t→∞ x(t) ≥ η and lim inf t→∞ y(t) ≥ η for all solutions with x(0) > 0, y(0) > 0. There are two methods to claim uniform persistence. One method is by applying the average Liapunov function [11] and the other is by studying boundary dynamics of the system [4,9]. We apply Theorem 5.1 of [4] to show persistence of (3) by verifying that the map restricted on the boundary is isolated and acyclic. Proof. Let Y be the boundary of R 2 + . Then Y is a closed subset of R 2 + and Y is invariant for (3). Clearly the map F induced by (3) is dissipative and the maximal We may assume x(0), y(0) > 0. Then (x(t), y(t)) ∈ K 0 so that 0 < x(t) < δ and 0 < y(t) < δ for t ≥ 0. It follows from the first equation of (3) that This inequality then implies lim t→∞ x(t) = ∞ and we obtain a contradiction.
Although both host and parasitoid populations can coexist indefinitely when βγx > 1, the population interactions may be complicated. By Lemma 2.2 we shall study detJ(E * ) given in (13) in more detail to determine whether the model changes its stability as the parameter β varies.
Since y * is a strictly increasing function of β, we have β = β(y * ) and may express detJ(E * ) as a function of y * . Using equilibrium equations (5) and (6), we can first write x * and then β as functions of y * . By omitting the * s, detJ(E * ) in (13) can be rewritten as .
The results given in terms of y in Theorem 2.5 can be carried over to the pa- , where dy dβ > 0 by (11), implies that ddetJ(E * ) dβ and dD(y) dy have the same signs. For the interior steady state E * = (x * , y * ), we know y * = y * (β) is strictly increasing in β. In particular, such β is unique. Let y * i = y * (β i ) for i = 1, 2. Hence, Then we have the following results.
It is difficult to solve E(z) = 0 defined in (15) analytically. But given λ and γ, it can be done numerically without difficulty. When detJ(E * ) crosses the line detJ(E * ) = 1 twice as in Theorem 2.5, where D(z) > 1 and (1 − γ)λ ≥ 1, it is expected that the interior steady state E * is destabilized first and then stabilized again as β is increased. If there exists a unique β 1 > 0 with detJ(E * ) = 1, then it is suspected that E * only undergoes a supercritical Neimark-Sacker bifurcation and E * cannot be stabilized again once it undergoes a stability change from stable to unstable.
To determine a supercritical Neimark-Sacker bifurcation, we first move the unique interior steady state (x * , y * ) to the origin by letting u = x − x * and v = y − y * : The Taylor expansion of system (30) about (0, 0) is where a 11 a 12 a 21 a 22 is the Jacobian matrix of system (3) evaluated at the interior steady state E * , and In order to put the linear part of (31) into the Jordan canonical form, we let µ ± iω be the eigenvalues of J(E * ) with µ 2 + ω 2 = 1 when β = β 1 , and let T = a 12 0 µ − a 11 ω .
Notice (34) is of the same form as the expression given in (15.24) of [6, page 475]. The existence of a supercritical Neimark-Sacker bifurcation is determined by where and all the partial derivatives are evaluated at (X, Y ) = (0, 0) with β = β 1 [6]. The assertion then follows directly from Theorem 15.31 of [6].
If β 2 derived in Theorem 2.6 is considered, then since ddetJ(E * ) dβ | β=β2 < 0, the conclusions in Theorem 3.1 are reversed. That is, if θ < 0, then (3) has an asymptotically stable invariant circle for β < β 2 and near β 2 . It is not easy to determine the value of θ given in (36) analytically. We use numerical examples to illustrate Theorem 3.1. We adopt the same parameter values as those in Figure 2(c)-(d), i.e., λ = 7.0 and γ = 0.85. The first β value such that detJ(E * ) = 1 is about 1.4575 and the θ value given in (36) is roughly θ = 8.825 > 0. Therefore, the system has an attracting invariant curve when β is larger than 1.4575 and close to 1.4575. As β is increased, detJ(E * ) crosses detJ = 1 again when β is around 4.351 with θ approximately equal to −44.5115 < 0. Therefore the system has an attracting invariant loop when β is less than 4.351 but close to it. See Figure  3  To study the impact of constant proportional refuge, we first observe that components of the interior steady state of (3) with respect to γ satisfy respectively and Therefore, the x component of the interior steady state is an increasing function of the fraction 1 − γ. On the other hand, since 1 + x * − λ = x * −x < 0 and dx * dγ < 0, the effect of γ on y * is not monotone as illustrated in (39). To investigate the stabilization/destabilization effects of the spatial refuge, we first let λ = 3.5 and β = 3.0 and use 1 − γ as the bifurcation parameter. The diagrams are shown in Figure 3(c)-(d). As the proportion 1 − γ of refuge increases, the interior steady state is stabilized and the magnitude of the host population increases. However, the magnitude of the parasitoid population at the coexisting steady state is not monotone with respect to 1 − γ. The same conclusion is obtained if λ = 7.0 and β = 3.0 as shown in Figure 3(e)-(f).

4.
Conclusions. There is a large number of theoretical studies in an attempt to understand the many different biological mechanisms that affect predator-prey interactions. One of the important aspects in this direction is that of prey's response to predators by using refuges. Although significant efforts have been carried out for populations that reproduce continuously, little attention has been made for populations with non-overlapping generations. One conclusion drawn from these studies of continuous-time models is that prey refuges can stabilize predator-prey interactions and prevent prey populations from extinction.
Insect parasitoids make up about 14% of the one million or so known insect species [7]. In addition, insects are subject to attack by parasitoids and are thus considered as hosts for the parasitoids. These insect species usually reproduce at certain time of their life stages and continuous-time models are not adequate to model such populations. It is the goal of this study to investigate the impact of host spatial refuge upon the host-parasitoid interaction using discrete-time systems.
Our model is based on a modified Nicholson-Bailey system and we assume that a constant proportion of the host is free from parasitism in each generation. We prove that both hosts and parasitoids can coexist if βγ(λ − 1) > 1, or equivalently, > 0 (cf. Theorem 2.3). If the inequality is reversed, then the parasitoid population goes extinct and the host population is stabilized as presented in Theorem 2.1. It is clear that the unique coexisting steady state is asymptotically stable when it just appears, i.e, when the quantity βγ(λ − 1) exceeds 1 but near 1. Using β as the bifurcation parameter, we show that the model can undergo a supercritical and then a subcritical Neimark-Sacker bifurcation in Theorems 2.6 and 3.1. It is also possible for the model to exhibit only the supercritical Neimark-Sacker bifurcation, depending on the parameter regimes. In the former situation, the parameter β can act as both a destabilizing and a stabilizing mechanism. For the latter case, increasing β is always destabilizing to the population interaction. It is clear that for each fixed β and λ, there exists a critical threshold γ 0 = 1 β(λ − 1) , such that the parasitoids go extinct if the constant proportionate refuges exceeds 1−γ 0 . To further investigate the effects of spatial host refuge, we prove that the host population at the coexisting steady state is always an increasing function of the proportionate refuges 1 − γ. However, the effect on the parasitoids is not monotone. Using 1 − γ as the bifurcation parameter, it is illustrated in Figure  3(c)-(f) that proportionate refuges can stabilize the host-parasitoid interaction.
To compare models (1) and (2)/(3), first notice that the host population in model (1) can grow unboundedly large if the product of the host growth rate and the proportionate host refuge exceeds one, λ(1 − γ) > 1. This is due to the exponential growth given by the host equation. Since a Beverton-Holt equation is adopted for the host dynamics in (2), the host population is always bounded and stabilized. For both models (1) and (2) the host population goes extinct and so are the parasitoids if the host's growth rate is smaller than one. There is no other non-extinction boundary steady state in system (1) while there exists another boundary steady state (x, 0) in model (3) for which the parasitoid population is extinct. The boundary steady state (x, 0) is globally asymptotically stable under the condition given in Theorem 2.1. In particular, the parasitoid population will become extinct if the constant proportional host refuge 1 − γ is large. This biological phenomenon seems reasonable but it cannot be captured in model (1). Furthermore, system (3) is able to undergo both a supercritical and a subcritical Neimark-Sacker bifurcation as one of the parameters is varied. From the graph given in [7, page 67], model (1) can however only exhibit a supercritical Neimark-Sacker bifurcation.
We conclude from this investigation that the constant proportionate refuges can drive the parasitoids to extinction if the proportion is large. If both populations are in the coexisting steady state, then an increase of the proportionality can increase the host population level. Moreover, proportionate refuges can stabilize the hostparasitoid interaction by damping irregular oscillations.