Mathematical Biosciences 261 (2015) 62–73

Contents lists available at ScienceDirect

Mathematical Biosciences journal homepage: www.elsevier.com/locate/mbs

Effects of additional food in a delayed predator–prey model Banshidhar Sahoo∗, Swarup Poria Department of Applied Mathematics, University of Calcutta, Kolkata, West Bengal, India

a r t i c l e

i n f o

Article history: Received 12 October 2014 Revised 3 December 2014 Accepted 17 December 2014 Available online 27 December 2014 Keywords: Predator–prey Habitat complexity Additional food Time delay Hopf-bifurcation

a b s t r a c t We examine the effects of supplying additional food to predator in a gestation delay induced predator–prey system with habitat complexity. Additional food works in favor of predator growth in our model. Presence of additional food reduces the predatory attack rate to prey in the model. Supplying additional food we can control predator population. Taking time delay as bifurcation parameter the stability of the coexisting equilibrium point is analyzed. Hopf bifurcation analysis is done with respect to time delay in presence of additional food. The direction of Hopf bifurcations and the stability of bifurcated periodic solutions are determined by applying the normal form theory and the center manifold theorem. The qualitative dynamical behavior of the model is simulated using experimental parameter values. It is observed that fluctuations of the population size can be controlled either by supplying additional food suitably or by increasing the degree of habitat complexity. It is pointed out that Hopf bifurcation occurs in the system when the delay crosses some critical value. This critical value of delay strongly depends on quality and quantity of supplied additional food. Therefore, the variation of predator population significantly effects the dynamics of the model. Model results are compared with experimental results and biological implications of the analytical findings are discussed in the conclusion section. © 2015 Elsevier Inc. All rights reserved.

1. Introduction The traditional approach of modeling predator–prey interaction is often based on an organism’s primary function within a food web (e.g., prey and predator). Additional food is an important component of most predators (e.g., coccinellid) diet, and although they receive less attention than prey in the scientific literature, these foods fundamentally shape the life histories of many predator species. The availability of suitable additional food (non-prey food) in a predator–prey system can have significant impact on the dynamics of the system. The consequences of providing additional food to predator and the corresponding effects on the predator prey dynamics and its utility in biological control (such as species conservation and pest management) have been a topic of great attention for many scientists. In recent years, many biologist, experimentalists, and theoreticians investigated the consequences of providing additional food to predators in predator– prey systems [1–8]. Huxel and McCann [2] investigated the impact of additional food on the stability of a simple food web model. Huxel et al. [3] proposed a food web model with variable allochthonous inputs which are either one type available to both consumer and predator or



Corresponding author. Tel.: +91 9933 590 675. E-mail addresses: [email protected] (B. Sahoo), [email protected] (S. Poria).

http://dx.doi.org/10.1016/j.mbs.2014.12.002 0025-5564/© 2015 Elsevier Inc. All rights reserved.

two distinct types, one for consumer and one for predator. Srinivasu et al. [6] examined the qualitative behavior of a predator–prey system in presence of additional food to the predator. They concluded that handling times for the available foods to the predator play the key role in determining the dynamical behavior of the system. Haque and Greenhalgh argued that alternative food source may play an important role in promoting the persistence of predator–prey systems [9]. Guin et al. [10] investigated the significant role of self and crossdiffusion coefficients in a prey-dependent predator–prey model in which predator has alternative source of food. Recently, a chaos control mechanism of a predator–prey system incorporating additional food to predator is reported by Sahoo and Poria [5]. Habitat complexity is the structural complexity of habitats. Habitat complexity can strongly mediate predator–prey interactions, affecting not only total predation rates, but also modifying selectivities for different prey species or size classes [11–16]. Pennings [17] and Grabowski [18] found that habitat complexity reduces encounter rates of predators with prey. For example, aquatic habitat becomes structurally complex in presence of submerged vegetation or aquatic weeds. It is observed that structural complexity of the habitat stabilizes the predator–prey interaction between piscivorous perch (predator) and juvenile perch and roach (prey) by reducing predator foraging efficiency. Luckinbill prolonged the coexistence of Paramecium aurelia (prey) and Didinium nasutum (predator) in laboratory system by increasing strength of habitat complexity using

B. Sahoo, S. Poria / Mathematical Biosciences 261 (2015) 62–73

methyl cellulose in the Cerophyl medium (nutrient) [19]. Therefore, it is important to incorporate the effect of habitat complexity when predator–prey interaction is studied by means of theoretical models. The models with delay are much more realistic, as in reality time delays occur in almost every biological situation [20–27]. It is more realistic to assume that the reproduction of predator after predating the prey will not be instantaneous, but mediated through some time lag which is required for gestation of the predator [28]. After predation, some amounts of energy in the form of biomass of prey assimilate into the predator’s energy also in the form of biomass. But this bio-physiological process is not simple; the conversion of prey energy to predator energy is not instantaneous, and several processes are involved in this mechanism. Caperon [29] examined the Isochrysis galbana growing in a nitrate limited chemostat, over two experiments (one of 56 days and another of 80 days). He observed that there is a time delay between the changes in substrate concentration and the corresponding changes in the bacterial growth rate. So, inclusion of time delay will certainly make the predator–prey model one step closer to real situation. Predators are predominantly valued for their ability to control prey and as a result to keep high levels of biodiversity. Effects of variation of predator population on the model are investigated. Additional food helps to control predator population in the model. Good quality or high quantity of additional food favor rapid growth of predator population. Majority of the predator–prey models has not incorporated the effects of habitat complexity as well as the effects of additional food in biological systems. The main goal of this paper is to investigate the role of additional food as a biological controller in a delayed predator–prey system with habitat complexity. The organization of the paper is as follows: We propose a food chain model incorporating effects of additional food to predator in Section 2. In Section 3 the existence conditions of interior equilibrium point are derived and role of additional food on the existence conditions are reported. Effects of additional food on the stability condition of the interior equilibrium point is derived. Role of additional food on the stability and direction of periodic solutions are investigated by using the normal form theory and the center manifold theorem due to Hassard in Section 4. Numerical simulation results are supplied in Section 5 in support of the theoretical analysis. Finally conclusion is drawn in Section 6. 2. The model We shall now derive the modified form of functional response considering the effects of additional food to predator with habitat complexity. The rate of prey consumption by an average predator is known as the functional response. This can be classified as: (a) prey dependent, when prey density alone determines the response; (b) predator dependent, when both predator and prey populations affect the response; and (c) multi-species dependent, when species other than the predator and its prey species influence the functional response [30]. The most commonly used functional response in a predator–prey model is the Holling Type-II functional response [31]. The Holling type-II functional response is defined as

f (x) =

ax , 1 + ahx

where f (x) is the amount of food consumed, x is the amount of food offered, a is a proportionality constant related to the attack rate, h is the handling time per food item. Since the existence of habitat complexity reduces the probability of capturing a prey by reducing the searching efficiency of predator and habitat complexity affects the attack coefficient [32]. Therefore, the attack coefficient a has to be replaced by a(1 − c), where c (0 < c < 1) is a dimension less parameter that measures the degree or strength of habitat complexity. We assume that the complexity is homogeneous throughout the habitat. Then following Kot [33], the total number of prey caught(V), is given

63

by

V = a(1 − c)Ts x,

where Ts = T − hV.

Here T is the total time, Ts is the available searching time. Solving for V we get the modified Holling type-II response function as

V=

Ta(1 − c)x . 1 + a(1 − c)hx

Since, predator’s functional response is defined as the number of prey caught by a predator at unit time, so the functional response in presence of habitat complexity is given by

f (x) =

a(1 − c)x . 1 + a(1 − c)hx

The term c (0 < c < 1) measures the strength of habitat complexity, which reduces the predation rate. Notice that for c = 0, there is no complexity, we get the original Holling Type II response function. Assuming density-dependent logistic growth of prey with intrinsic growth rate r, the predator–prey model is of the form

  x x dx a(1 − c)xy = rx 1 − = rx 1 − − − f (x)y, dt k 1 + a(1 − c)hx k dy θ a(1 − c)xy = − dy = g(x)y − dy, dt 1 + a(1 − c)hx

(1)

where x denotes the density of prey, y is the density of the predator, k is the carrying capacity of the prey in ecosystem, θ is the conversion efficiency of prey into the predator and d is the mortality rate of predator. One aspect of habitat manipulation is the addition of floral resources to agro-ecosystems to provide additional food to predators, potentially enhancing their fitness and efficacy [34,35]. These techniques could also be used to improve the success of classical biological control attempts [36]. The mirid predator Macrolophus pygmaeus is a natural enemy of major economic importance for the control of white flies and other small arthropod pests in Europe [37–40]. It is observed that provision of a minimum of 40 eggs per individual predator for three days is required for optimal development and reproduction of this mirid predator. Providing the predator with lower quantities of eggs resulted in higher mortality, slower development and lower adult weights [40]. Since provision of eggs to the predator proved expensive, experiments were conducted to find if pollen can be a supplementary food for this predator [40]. It is observed that food consisting of 10 eggs and 15 mg of pollen was needed for optimal development of the predator, which was relatively a cheaper alternative [40,41]. Thus availability of additional food of a fixed quality appears to be vital in the development, conservation and sustainability of the species both ecologically and economically. Both theoretical studies [6] and experimental results [42–44] established that provision of additional food to predators mediates indirect interactions between the species of the ecosystem, ultimately affecting the population dynamics of the predator and prey. The above fact motivate us to incorporate effects of additional food in our model. We now modify the model (1) by supplying “additional food” [6,8] to predator. The predator is provided with a constant additional food in case of extinction of the prey as predator has alternate source of food other than the prey available to it. We make the following assumptions: (a) The predator is provided with additional food of constant biomass A which is assumed to be distributed uniformly in the habitat. The constant biomass assumption is valid for many arthropod predators because they can feed on plant-provided alternative food sources such as pollen or nectar which approximately remains constant [1]. Therefore, the predator is a generalist with a resource other than the prey available to it. (b) The number of encounters per predator with the additional food is proportional to the density of the additional food. The

64

B. Sahoo, S. Poria / Mathematical Biosciences 261 (2015) 62–73

proportionality constant characterizes the ability of the predator to identify the additional food [6]. Let us consider that h1 represents handling time (time spent on processing a food item) of the predator y per prey item in presence of additional food and e1 represents the ability (the rate at which the consumer encounters food items per unit of food density) of the predator to detect the prey respectively. Let the constant h2 represents the handling time of the predator y per unit quantity of additional food and the term e2 represents the ability of the predator y to capture additional food. Therefore, the functional response with additional food may be written in the form

f (x) =

e1 a1 (1 − c)x , 1 + e2 h2 A + e1 h1 a1 (1 − c)hx

y(φ) = ψ2 (φ) ≥ 0, φ ∈ [−τ , 0), ψ1 (0) > 0, ψ2 (0) > 0, where C+ = {(ψ1 (φ), ψ2 (φ)) ∈ C ([−τ , 0], R2+ ) : x(φ) = ψ1 (φ), y(φ) = ψ2 (φ)} is the Banach space of continuous functions form the interval [−τ , 0] into R2+ = {(x, y) : x(t) > 0, y(t) > 0}.

a2 a1 (1 − c)x , b1 + αμA + a1 (1 − c)hx

where a2 = 1/ h1 , b1 = 1/e1 h1 , μ = e2 /e1 and α = h2 / h1 .

2.1. Existence and positive invariance

a2 (a1 /b1 )(1 − c)x . 1 + α(μA/b1 ) + (a1 /b1 )(1 − c)hx

The system (3) can be written as X¯˙ = F (X¯ ) with X¯ (t) = X¯ (t, X¯ (φ)), where X¯ (φ) = (ψ1 (φ), ψ2 (φ)) ∈ C+ , ψi (0) > 0, for i = 1, 2, 3 and X¯ = (x, y)T ∈ R2+ . Now F (X¯ ) is given by

Hence, the functional response is modified into the form

f (x) =

a2 a(1 − c)x , 1 + αξ + a(1 − c)hx

where a =

a1 b1



and ξ = μb A . The term μA represents effectual addi1

tional food for the predators y. However, the term α =

h2 h1

is directly

proportional to the handling time h2 and so it represents the “quality” of additional food. If the relation h2 < h1 holds, then the predator can easily capture additional food than prey species and it implies that the additional food is of high quality. Therefore, for high quality of additional food α is less than 1 [4]. The non-dimensionalized term ξ = μb A is directly proportional to the biomass of the additional food 1

(A) and so it represents the “quantity” of additional food. Hence, the functional response in terms of quality and quantity of additional food is equivalent to

f (x) =

a(1 − c)x . 1 + αξ + a(1 − c)hx

If  represents the efficiency (effectual conversion rate of food to predator) with which the food consumed by predator gets converted into the predator biomass and therefore, the growth rate of predator y is of the form

g(x) =

θ a(1 − c)(x + ξ )y , 1 + αξ + a(1 − c)hx

where, θ =  a2 =  / h1 . Incorporating the above facts, the modified form of the model (1) is the following:

 x a(1 − c)xy dx = rx 1 − − , dt k 1 + αξ + a(1 − c)hx dy θ a(1 − c)(x + ξ )y = − dy. dt 1 + αξ + a(1 − c)hx

⎞  a(1 − c)xy x − rx 1 − ⎜ k 1 + αξ + a(1 − c)hx ⎟ ⎟, F = F (X¯ ) = ⎜ ⎝ θ a(1 − c)(x(t − τ ) + ξ )y(t − τ ) ⎠ − dy 1 + αξ + a(1 − c)hx(t − τ ) where F : C+ → R2 and F ∈ C ∞ (R2 ). It is easy to show that whenever choosing X¯ (φ) ∈ C+ such that Xi = 0, then Fi (X¯ )|X =0,X¯ t ∈C+ ≥ 0, for i = 1, 2, 3. Due to Lemma-2 of Yang i et al. [23] it is observed that any solution of X¯˙ = F (X¯ ) with X¯ (φ) ∈ C , say X¯ (t) = X¯ (t, X¯ (φ)), is such that X¯ (t) ∈ R2+ for all t > 0.

(2)

It is important and easy to show that model (2) can be reduced to model (1) by the simple variable alternation 1+aαξ → a. This indicates that the proposed supply of an additional food is equivalent to the

+

2.2. Boundedness Boundedness is a necessary condition for the system (3) to be biologically realistic. The following propositions ensure the boundedness of the system (3) with non-zero time lag. Proposition 2.1. The prey population is always bounded from above. Proof. From the first equation of the system (3) it follows that

 x dx ≤ rx 1 − , dt k

 a2 a(1 − c)(x + ξ ) . 1 + αξ + a(1 − c)x

That is, g(x) =

(3)

x(φ) = ψ1 (φ) ≥ 0,

Therefore, it can be written as

f (x) =

 x a(1 − c)xy dx = rx 1 − − , dt k 1 + αξ + a(1 − c)hx dy θ a(1 − c)(x(t − τ ) + ξ )y(t − τ ) = − dy. dt 1 + αξ + a(1 − c)hx(t − τ ) The initial conditions for the system (3) are chosen as:

(1/ h1 )a1 (1 − c)x i.e., f (x) = . (1/e1 h1 ) + (h2 / h1 )(e2 /e1 )A + a1 (1 − c)hx f (x) =

change of the attack rate a. It is also clear that increase of additional food reduces the predatory attack rate to prey. As the reproduction of the predator population after predation the prey will not be instantaneous, but mediated by some constant time lag τ > 0 for gestation of predator, we incorporate time delay in the predator equation into the model (2). Therefore the model (2) becomes the following

limsupt→∞ x(t) ≤ k.

Proposition 2.2. If the quantity of additional food ξ satisfies the condid tion 0 < ξ < aθ (1−c ) , then the system (3) is bounded. Proof. We first consider the function χ = θ x(t − τ ) + y(t). Calculating the time derivative along the solution of the system (3), for t > t0 − τ letting M = min{r, d − aξ θ (1 − c)} we have

dχ θ a(1 − c)ξ y(t − τ ) = rθ x(t − τ ) + < −rθ x(t − τ ) dt 1 + αξ + ah(1 − c)x(t − τ ) −{d − θ aξ (1 − c)}y + 2rθ x(t − τ ) < −Mχ + 2rθ k, where we used Proposition 2.1 for the inequality x(t − τ ) < k for all t > t0 − τ . Thus there exists a positive constant , such that χ(t) <

for all large t.

B. Sahoo, S. Poria / Mathematical Biosciences 261 (2015) 62–73

(i) If C 2 < D2 , then the interior equilibrium point E∗ of the system (3) is locally asymptotically stable for τ < τ ∗ and unstable when τ > τ ∗ where

3. Stability of the equilibrium and local Hopf bifurcations In this section, we discuss the stability of the coexisting equilibrium point of the system (3) only. The interior equilibrium point αξ )−θ a(1−c)ξ and y∗ = of the system is E∗ (x∗ , y∗ ), where x∗ = d(1+ a(1−c)(θ −hd)

τ∗ =

r(k−x∗ )[1+αξ +a(1−c)hx∗ ] . The interior equilibrium point (x∗ , y∗ ) will be biak(1−c) )(θ −hd)−d , c < 1 − d and ologically meaningful if 0 < ξ < ka(α1−c d−aθ (1−c) ak(θ −hd) d hd + ak < θ < 1. The existence conditions of interior equilibrium

point depends on quality and quantity of additional food supplied to predator. Let us consider X (t) = x(t) − x∗ , Y (t) = y(t) − y∗ as perturbed variables. Then from the system (3) we obtain the time evolution rule for X (t), Y (t) after linearization as:

˙



X (t) X (t) X (t − τ ) = A1 + A2 , Y (t) Y (t − τ ) Y˙ (t)

(4)

where, ⎛ ⎜ A1 = ⎝

rx∗ [ahk(1 − c) − (1 + αξ ) − 2ah(1 − c)x∗ ] k[1 + αξ + a(1 − c)hx∗ ]



d

θ

+

1 + αξ + a(1 − c)



0

⎜ A2 = ⎝ θ r[1 + αξ − a(1 − c)hξ ](k − x∗ ) k[1 + αξ + a(1 − c)hx∗ ]

hx∗

⎟ ⎠

−d

0

and



ξ

⎞ 0 ⎟ ⎠. d

det(A1 + A2 e−λτ − λI) = 0,

λ − (A + Be−λτ )λ + (C + De−λτ ) = 0, 2

(5)

where,

rx∗ [ahk(1 − c) − (1 + αξ ) − 2ah(1 − c)x∗ ] − d, k[1 + αξ + a(1 − c)hx∗ ] B = d, rdx∗ [akh(1 − c) − (1 + αξ ) − 2ah(1 − c)x∗ ] C=− , k[1 + αξ + a(1 − c)hx∗ ] rdx∗ [akh(1 − c) − (1 + αξ ) − 2ah(1 − c)x∗ ] D= k[1 + αξ + a(1 − c)hx∗ ]

d ξ − − + . θ 1 + αξ + a(1 − c)hx∗

(k − x∗ )rθ (1 + αξ − a(1 − c)hξ ) . × ∗ k(1 + αξ + a(1 − c)hx ) A=

ω∗

cos−1



ω∗2 (D − AB) − CD . D2 + ω∗2 B2

When τ = τ ∗ , a Hopf-bifurcation occurs as τ passes through the critical value τ ∗ [45]. (ii) If C 2 > D2 , then the interior equilibrium point E∗ of the system (3) is locally asymptotically stable for all τ ≥ 0. Proof. See Appendix A. 4. Direction and stability of Hopf bifurcation In this section, we have studied the direction of the Hopf bifurcation and stability of bifurcated periodic solutions arising through Hopf bifurcation by applying the normal form theory and center manifold theorem introduced by Hassard [45]. Throughout this section, we always assume that system (3) undergoes Hopf bifurcation at the positive equilibrium point E∗ (x∗ , y∗ ) for τ = τ ∗ , at which the characteristic Eq. (5) has a pair of purely imaginary roots ±iω∗ . The conditions for direction and stability of Hopf bifurcation are summarized in the following theorem.

(i) The direction of Hopf bifurcation is determined by the sign of μ2 : if μ2 > 0 (μ2 < 0), then the Hopf bifurcation is supercritical (subcritical) and the bifurcating periodic solutions exist for τ > τ ∗ (τ < τ ∗ ). (ii) The stability of bifurcated periodic solutions are determined by β2 : the periodic solutions are stable if β2 < 0 and unstable if β2 > 0. (iii) The period of bifurcated periodic solutions are determined by T2 : the period increases if T2 > 0 and decreases if T2 < 0. Proof. The values of μ2 , β2 and T2 are determined in Appendix B. 5. Simulation results

The characteristic equation of the system (4) without time delay is

λ2 − [A + B]λ + [C + D] = 0,

1

Theorem 4.1.

The characteristic equation of system (4) is

i.e.,

65



(6)

(A+B)± (A+B)2 −4(C+D)

and the corresponding eigenvalues are λ1,2 = . 2 Since, for the stability of the equilibrium point, real parts of the above two roots must be negative. Here (C + D) > 0 with supply of additional food or without additional food, so eigenvalue of (6) will have negative real parts if (A + B) < 0. After some simple calculations, one can show that the non-delayed system will be stable if the conditions of the following theorem hold. Theorem 3.1. The system (3) is locally asymptotically stable without delay at the coexisting equilibrium point E∗ if ka(1−c)(θ −hd)−d α d−aθ (1−c) , θ +hd 1 − ahk(θ −hd) < c < 1 − ak(θd−hd) , ), α > aθ (1−c d d hd + ak < θ < 1.

We verify our theoretical analysis with the help of numerical simulation results taking parameter values from the experimental data (Table 1). Initial values are kept fixed in all cases at (15, 5.83) [19]. The parameters α (quality of additional food) and ξ (quantity of additional food) will vary throughout numerical simulations. For the numerical simulations, we use MATLAB 2009a software. From Theorem 3.1, we see that for the stability of coexisting equilibrium point, the quantity of additional food has to be bounded by a value which depends on quality of additional food. For the above set of parameter values, α should be greater than 0.005. Therefore, prey and predator coexist in an ecosystem if we increase predator population slightly. In absence of additional food (i.e., for α = 0, ξ = 0), the system (3) is asymptotically stable at interior equilibrium E∗ (x∗ , y∗ ) = (253.9056, 97.8867) when τ = 0, since A + B = −0.3398 < 0 and C + D = 1.5807 > 0. When τ = 0, it is observe that the parameter Table 1 Ranges and default values of parameters with their sources. Parameters

Description

Range

Default value

r k c h a d

Intrinsic growth rate of prey Carrying capacity Strength of habitat complexity Handling time Attack rate Predator death rate Conversion efficiency

3.3 [46] 898 [19] 0 τ ∗ in absence of additional food. For simulating the dynamics of the system (3) in presence of additional food, we take α = 2, ξ = 0.1 < 1.2716 satisfying the existence conditions for coexisting equilibrium point (Theorem 3.1). With such additional food the system (3) is asymptotically stable at E∗ = (304.5592, 108.2160) for τ = 0, since A + B = −0.5216 < 0 and C + D = 1.4371 > 0. Therefore, the system (3) is asymptotically stable for this type of additional food without delay. For τ = 0, we plot dynamic behavior of x(t) versus x(t − τ ) in Fig. 2(a) and y(t) versus y(t − τ ) in Fig. 2(b). The time evolution of x(t − τ ) and y(t − τ ) are presented in Fig. 2(c). All these figures show the equilibrium point is asymptotically stable for τ = 0.5 < τ ∗ = 1.192. But, for τ = 1.2 > τ ∗ = 1.192, the equilibrium point is unstable when α = 2, ξ = 0.1 (Fig. 3). Therefore, from Figs. 2 and 3 it is evident that coexisting equilibrium point is stable for τ = 0.5 < τ ∗ = 1.192 and unstable for τ = 1.2 > τ ∗ = 1.192. Thus, the system (3) undergoes a Hopf bifurcation at E∗ when τ = τ ∗ = 1.192 with additional food having α = 2, ξ = 0.1.

Next we take α = 3 and ξ = 0.1 < 0.8470 satisfying the existence conditions for coexisting equilibrium point (Theorem 3.1). From Fig. 4, one can observe that the system (3) is asymptotically stable at the equilibrium point E∗ when τ = 2.5 < τ ∗ = 3.802 and unstable when τ = 4.3 > τ ∗ = 3.802. Therefore, the system (3) undergoes a Hopf bifurcation at E∗ = (329.9498, 112.2195) when τ = τ ∗ = 3.802 for α = 3, ξ = 0.1. We have constructed a bifurcation diagram to observe the dynamics of the system (3) in Fig. 5 with respect to gestation delay τ for different values of α and ξ . Figs. 2–4 are verified again with Fig. 5(a) and (b). We also present the dynamics of the system (3) in Fig. 5(c) for high quality of additional food (α = 0.8 < 1, ξ = 0.2). From Fig. 5(c), it is evident that the system (3) undergoes a Hopf bifurcation at E∗ when τ = τ ∗ = 0.8298. In Fig. 6, the stability region is plotted in the (τ , ξ ) plane for α = 2. Fig. 6 indicates that small quantity of additional food is required for stability of the system with low gestation delay τ . We observe that for stability of the system with higher delay, the quantity of additional food supply ξ should be increased (Table 2). The stability region in the (τ , α) plane is also plotted in Fig. 7 for ξ = 0.1. For small delay high quality of additional food is required, but for high delay low quality of additional food is sufficient for stability of coexisting equilibrium point. Therefore, from Figs. 6 and 7 we conclude that stable steady state can be reached with the supply of additional food for all delay.

B. Sahoo, S. Poria / Mathematical Biosciences 261 (2015) 62–73

67

(a) x(t−τ)

800 600

(b) 150

α=2, ξ=0.1 τ=1.2

y(t−τ)

1000

400

α=2, ξ=0.1 τ=1.2

100

50

200 0 100

200

300

400

500

600

700

800

0 0

900

20

40

60

x(t)

80

100

120

140

y(t)

Populations

(c) 1000

α=2, ξ=0.1 τ=1.2

500

Prey x

Predator y

0 0

100

200

300

400

500

600

700

800

Time(days) Fig. 3. Behavior of the system (3): (a) x(t) versus x(t − τ ), (b) y(t) versus y(t − τ ) (a) τ = 0.2, and (c) depicts that the system (3) is oscillatory when τ = 1.2 > τ ∗ = 1.192.

(a)

(b)

900

900

800

800

α=3, ξ=0.1 τ=2.5

600 500 400

α=3, ξ=0.1

700

Populations

Populations

700

Prey x

300 200

τ=4.3

600 500

Prey x

400 300 200

100

100

Predator y

0 0

200

400

600

Predator y

0 0

800

100

Time(days)

200

300

400

Time(days)

Fig. 4. Behavior of the system (3) for different τ : (a) τ = 2.5, (b) τ = 4.3 for α = 3, ξ = 0.1. This figure shows that the system (3) is steady when τ = 2.5 < τ ∗ = 3.802 and oscillatory when τ = 4.3 > τ ∗ = 3.802.

(a)

(b)

550 500

(c)

450

α=2, ξ=0.1

400

600

α=3, ξ=0.1

550

Prey x

300 250

Prey x

450

Populations

350

400

Populations

Populations

450

350

300 250 200

200 150

400 350

Prey x

300 250 200

150

Predator y

100 50 0

α=0.8, ξ=0.2

500

Predator y

150

100 0.5

1

τ (days)

1.5

2

0

Predator y

100 1

2

3

τ (days)

4

5

50 0

0.5

1

1.5

2

τ (days)

Fig. 5. Bifurcation diagram of the prey and predator populations with respect to the gestation delay τ for different α and ξ . Other parameters are as in Table 1. (a) τ ∗ = 1.192 for α = 2 > 1, ξ = 0.1, (b) τ ∗ = 3.802 for α = 3 > 1, ξ = 0.1 and (c) τ ∗ = 0.8298 for α = 0.8 < 1, ξ = 0.2.

68

B. Sahoo, S. Poria / Mathematical Biosciences 261 (2015) 62–73 3.5

0.16

Steady region

Quality of additional food α

Quantity of additional food ξ

0.18

0.14 0.12 0.1

Oscillatory region 0.08 0.06 0.04 0.02 0

1

2

3

4

5

6

7

8

Steady region

3 2.5 2

Oscillatory region 1.5 1 0.5 0

9

5

10

τ(days)

Table 2 Effects of supplying additional food on existence of predator–prey population.

0.5

1.5

2

Additional food

Predator–prey coexistence

αξ αξ αξ αξ αξ αξ αξ αξ αξ

Both exist Both exist Prey exists, predator extinct Both exist Both exist Prey exists, predator extinct Both exist Both exist Prey exists, predator extinct

= 0.4 = 1.6 = 3.6 = 0.4 = 1.6 = 3.2 = 0.4 = 1.6 = 2.4

We increase the strength of complexity c from 0.45 to 0.63, satisfying the conditions of Theorem 3.2(ii). We observe that the interior equilibrium point E∗ of the system (3) becomes asymptotically stable for all admissible values of τ using any supply of additional food levels, which is shown in Fig. 8. In this case, the additional food has no effect on the delayed system. 6. Conclusions Effect of variation of predator population on the model is an interesting area of research. Predator population can be controlled

30

supplying additional food to predator. Considering the effects of supplying additional food to predator a delay induced predator–prey model with habitat complexity is proposed. We obtain sufficient conditions for asymptotic stability of the coexisting equilibrium point and notice that it is strongly dependent on quality and quantity of additional food. Therefore, variation of predator population by supplying additional food has strong influence on the model. We observe that there exists a critical value of the delay parameter below which the system has steady state and above which the system has oscillatory solution. Therefore, Hopf bifurcation occurs in the system at this critical delay value. The Hopf bifurcation point of the system changes with the quality and quantity of additional food. Biologically predator population has significant effect in the appearance of oscillatory solution of the system. By applying the normal form theory and the center manifold theorem, we determine the stability and direction of the bifurcating periodic solutions. For numerical simulation, we choose the parameter values of the P. aurelia (prey) and D. nasutum (predator) experiment. We predict that in absence of additional food, when the gestation delay τ exceeds the critical value τ ∗ (= 0.3412 days approx.) for a biologically meaningful parameter set, the system (3) bifurcates from stable focus to limit cycle. But in presence of additional food, the critical values τ ∗ depends on quality and quantity of additional food, which is shown in Fig. 5. Also, population size of the system depends 900

800

Populations

α=2, ξ=0.1 τ=1.5

700

Populations

25

Fig. 7. Graph of τ (days) versus quality of additional food α for ξ = 0.1 and other parameters are as in Table 1. The τ α -plane is divided into steady and oscillatory regions. This figure gives the minimum values of α for different τ so that the system (3) becomes steady within a stipulated time.

900

600 Prey x

500 400 300 200

800

α=2, ξ=0.1

700

τ=3

600 Prey x

500 400 300 200

100 0 0

20

τ(days)

Fig. 6. Graph of τ (days) versus quantity of additional food ξ for α = 2 and other parameters are as in Table 1. The τ ξ -plane is divided into steady and oscillatory regions. This figure gives the minimum values of ξ for different τ so that the system (3) becomes steady within a stipulated time.

τ

15

100

Predator y 50

100

Time(days)

150

0 0

Predator y 50

100

150

Time(days)

Fig. 8. Time evolution of the prey and predator populations of the system (3): (a) τ = 1.5, (b) τ = 3 for c = 0.63, α = 2, ξ = 0.1 and other parameters are as in Table 1. This figure shows that for all permissible delay the system (3) converges to the equilibrium E∗ = (452.7852, 120.6868) with initial condition (15, 5.83).

B. Sahoo, S. Poria / Mathematical Biosciences 261 (2015) 62–73

crucially on quality and quantity of additional food if the gestation delay period exceeds the critical value τ ∗ . Fig. 1(b) shows 3–4 oscillations over a period of 15–20 days in absence of additional food. But, in presence of additional food, the number of oscillations will decrease over the same period (as observed from Figs. 3(c) and 4(b)). It is well known that the predator–prey interaction exhibits oscillations in nature. Oscillatory population may be driven to extinction in presence of environmental stochasticity if the population density becomes very low. Therefore, we always prefer steady state behavior of a system for conservation of species. Our theoretical study (Theorem 3.2(ii)) and simulation result (Fig. 8) indicate that one of the reasons for fluctuations in the population size is the gestation delay. The fluctuation of the population size can be completely removed either by supplying additional food suitably or by increasing the strength of habitat complexity. Predator population control is possible through the supply of additional food. Therefore, only through predator control we can stabilize a predator–prey system. Our results are consistent with the observations that species losses within ecosystem heavily depends on predator population [51–55]. Therefore, supply of suitable additional food may be very useful way to maintain biodiversity of an ecological system. Luckinbill [19] was successful in developing an experiment to prolonged coexistence of Didinium (predator) and Paramecium (prey) in a homogeneous environment by the addition of methyl cellulose. Luckinbill reported coexistence of predator–prey in the medium with 0.75 g/l H2 O (half-strength medium) for the methyl cellulose and Cerophyl mixture. In systems with full-strength medium Didinium (predator) became extinct and the prey survived. The increase of methyle cellulose concentration in the medium reduces the predatory attack rate in the experiment. Therefore, coexistence is possible if predatory attack rate is moderately low but very low attack rate results extinction of predator. Again, Veilleux [56] studied the interaction between the predator D. nasutum and its prey P. aurelia in the laboratory. The Paramecium was fed the bacterium Aerobacter aerogenes (additional food). The author pointed out that in methyl cellulose regime, the outcome of the interaction depended upon the concentration of the bacterial nutrient Cerophyl in the medium. The Cerophyl concentration regulates the amount of bacterial food available to the Paramecium and, therefore, can be related to the degree of Paramecium starvation. This study reported that at low Cerophyl concentrations, the predator and prey coexisted at a numerically stable equilibrium, on the other hand, at high Cerophyl concentrations, the species were unable to coexist. This experiment also established that the coexistence of

species is possible if predatory attack rate is moderately low but for very low attack rate predator extinction is noted. In our model, the predatory attack rate decreases with the increase of additional food in the system. Our numerical simulation results show that both species exist for small quantity of additional food but predator extinction is observed at high quantity of additional food (Table 2 and Fig. 9). Therefore, results of our model is consistent with the experimental results of [19] and Veilleux [56]. Results presented in the present paper provided a useful platform to understand the role of additional food as a population controller in a delayed food chain model in presence of habitat complexity. A possible interesting aspect for future development is to study the effects of predator control in a metapopulation model of prey predators with species dispersion and migration. Appendix A. Proof of Theorem 3.2 For the delay-induced system (3), the interior equilibrium E∗ will be asymptotically stable if all the roots of the corresponding characteristic Eq. (5) have negative real parts. We start with the assumption that E∗ is asymptotically stable in case of non-delayed system and then find conditions for which E∗ remains stable for all delays [57]. By Rouche’s theorem [58] and the continuity in τ , the transcendental Eq. (5) has roots with positive real parts if and only if it has purely imaginary roots. Now, we shall be able to find conditions for which eigenvalues have negative real parts.

Let,

λ(τ ) = μ(τ ) + iω(τ ),

where μ and ω are real. As the interior equilibrium E∗ of the nondelayed system is stable, we have μ(0) < 0. By continuity, if τ (> 0) is sufficiently small, we still have μ(τ ) < 0 and E∗ is still stable. The change of stability will occur at some values of τ for which μ(τ ) = 0, ω(τ ) = 0, that is λ will be purely imaginary. Let τ ∗ be such that, μ(τ ∗ ) = 0, ω(τ ∗ ) = ω∗ = 0, so that λ = iω(τ ∗ ) = iω∗ . In this case, the steady state loses stability and eventually becomes unstable when μ(τ ) becomes positive. In other words, if such an ω(τ ∗ ) does not exist, then the steady state E∗ will be stable for all τ . Now, if iω is a root of (5), separating real and imaginary parts, we have

−ω2 + C = −D cos ωτ + Bω sin ωτ ,

(b)

τ=0.5

(c) 900

τ=1.5

800

−−−−− prey x

500

− − − predator y

400 300

700

600

−−−−− prey x

500

− − − predator y

400 300

Populations

600

600

300 200

100

100

100

40

60

Time

80

100

0 0

20

40

60

Time

80

100

− − − predator y

400

200

20

−−−− prey x

500

200

0 0

τ=2

800

700

Populations

700

(8)

Squaring and adding the two Eqs. (7) and (8), we get

900

800

(7)

ωA = −(Bω cos ωτ + D sin ωτ ).

(a) 900

Populations

69

0 0

20

40

60

80

100

Time

Fig. 9. Time evolution of the prey and predator populations of the system (3) using parameter values taken from Table 1: (a) αξ = 0.4 (magenta line, --), αξ = 1.6 (blue line, --), αξ = 3.6 (red line, -) (b) αξ = 0.4 (magenta line, --), αξ = 1.6 (blue line, --), αξ = 3.2 (red line, -) and (c) αξ = 0.4 (magenta line, --), αξ = 1.6 (blue line, --), αξ = 2.4 (red line, -) with initial condition (15, 5.83). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

70

B. Sahoo, S. Poria / Mathematical Biosciences 261 (2015) 62–73

ω4 + (A2 − 2C − B2 )ω2 + (C 2 − D2 ) = 0. Note that A2 − 2C − B2 = [ C2

(9)

rx∗ [ahk(1−c)−(1+αξ )−2ah(1−c)x∗ ]

D2

k[1+αξ +a(1−c)hx∗ ]

]2 is always

positive. If > then all roots of Eq. (9) have negative real parts for all delay and the equilibrium E∗ is locally asymptotically stable. If C 2 < D2 , Eq. (9) has unique positive root. It follows that Eq. (9) will have a positive root ω∗ . This implies that the characteristic Eq. (5) will have a pair of purely imaginary roots ±iw∗ . The corresponding critical value of time delay τj ∗ , for j = 0, 1, 2, 3, . . . at ω = ω∗ is

τj ∗ =



1

ω

cos−1 ∗

dλ dτ

−1

=



2

2λ − A B − − − C λ − λ3 Dλ − Bλ2

(10)

τ . λ

(11)

At τ = τ ∗ , λ = iω∗ and then



−1  d(Reλ(τ )) dλ sign = sign Re , dτ dτ τ =τ ∗ ,λ=iω∗ τ =τ ∗ ,λ=iω∗  2  (A − 2C − B2 )ω∗2 + 2ω∗4 = sign > 0, D2 ω∗2 + ω∗4 B2

Let τ = τ ∗ + μ, μ ∈ R such that Hopf bifurcation occurs at μ = 0. By change of variables X (t) = x(t) − x∗ and Y (t) = y(t) − y∗ and by rescaling the time t → τt , in system (3) the transformed system becomes the following functional differential equation (FDE) in C = C ([−1, 0], R2 ):

u˙ (t) = Lμ (ut ) + f (μ, ut ).

(12)

˙



X (t) X (t) X (t − 1) i.e., = τ M(τ ) + τ N(τ ) + τ f (X, Y ), Y (t) Y (t − 1) Y˙ (t) where, d

θ

+

−d,

0

⎜ N(τ ) = ⎝ θ r[1 + αξ − a(1 − c)hξ ](k − x∗ ) k[1 + αξ + a(1 − c)hx∗ ]

0 d

⎞ ⎟ ⎠,

⎞ a(1 − c)xy rx2 − ⎟ ⎜− 1 + αξ + a(1 − c)hx ⎟ ⎜ k and f (X, Y ) = ⎜ ⎟. ⎝ θ a(1 − c)x(t − 1)y(t − 1) ⎠ 1 + αξ + a(1 − c)hx(t − 1) ⎛

ξ



1 + αξ + a(1 − c)hx∗ ⎟ ⎠,

0



0 −1

dη(ϑ , μ)φ(ϑ).

(13)

η(ϑ , μ) = (τ ∗ + μ)[M(τ ∗ + μ)δ(ϑ) − C (τ ∗ + μ)δ(ϑ + 1)], where δ is the Dirac delta function and therefore,

⎧ ∗ (τ + μ)M(τ ∗ + μ), ⎪ ⎨ 0, η(ϑ , μ) = ⎪ ⎩ ∗ −(τ + μ)C (τ ∗ + μ),

ϑ = 0, ϑ ∈ (−1, 0), ϑ = −1.

For φ ∈ C 1 ([−1, 0], R2 ), define the operator A(μ) as

⎧ dφ(ϑ) ⎪ ⎪ ⎨ dϑ , A(μ)φ(ϑ) =  0 ⎪ ⎪ ⎩ dη(ζ , μ)φ(ζ ), 

ϑ ∈ [−1, 0), ϑ = 0.

ϑ ∈ [−1, 0), h(μ, φ), ϑ = 0.

0,

⎞ rφ12 (0) a(1 − c)φ1 (0)φ2 (0) − ⎟ ⎜− k 1 + αξ + a(1 − c)hφ1 (0) ⎟ ⎜ with h(μ, φ) = ⎜ ⎟. ⎠ ⎝ θ a(1 − c)φ1 (−1)φ2 (−1) 1 + αξ + a(1 − c)hφ1 (−1) ⎛

Appendix B. Proof of Theorem 4.1





R(μ)φ(ϑ) =

Remark. The delay τ will not affect the stability of the system, when Eq. (9) has no real root under the condition (A2 − 2C − B2 ) > 0, C 2 > D2 . The system is asymptotically stable for τ ≥ 0.

rx∗ [ahk(1 − c) − (1 + αξ ) − 2ah(1 − c)x∗ ] k[1 + αξ + a(1 − c)hx∗ ]

Lμ (φ) =

and

)) | Thus, d(Redλ(τ τ =τ ∗ ,λ=iω∗ > 0. τ Therefore, the transversality condition holds and a Hopf bifurcation occurs at τ = τ ∗ .

⎜ M(τ ) = ⎝

By the Riesz representation theorem, there exists a 2×2 matrix function η(ϑ , μ), ϑ ∈ [−1, 0], whose elements are functions of bounded variation such that, for φ ∈ C

−1

since (A2 − 2C − B2 ) > 0.



Lμ (φ) = (τ ∗ + μ)[M(τ ∗ + μ)φ(0) + N(τ ∗ + μ)φ(−1)].

In fact, we can choose

ω∗2 (D − AB) − CD 2π j + ∗ . ω D2 + ω∗2 B2

)) | We now verify that the transversality condition d(Redλ(τ τ =τ ∗ > 0 τ holds for j = 0. This will show that there exists at least one eigenvalue with a positive real part for τ > τ ∗ . Moreover, the conditions for the existence of Hopf bifurcation are then satisfied and a periodic solution exists. Differentiating (5) with respect to τ , it follows that



For φ = (φ1 , φ2 )T ∈ C ([−1, 0], R2 ), we define

(14)

Then the system (12) is equivalent to

u˙ (t) = A(μ)ut + R(μ)ut ,

(15)

where ut (ϑ) = u(t + ϑ) for ϑ ∈ [−1, 0]. For ψ ∈ C 1 ([0, 1], (R2 )∗ ), define

⎧ dψ(s) ⎪ ⎪ , ⎨− ds ∗ A ψ(s) =  0 ⎪ ⎪ ⎩ dη(t, 0)ψ(−t),

s ∈ (0, 1], s = 0.

−1

and a bilinear inner product

¯ 0)φ(0) −

ψ(s), φ(ϑ) = ψ(



0 −1

 ϑ ζ =0

¯ − ϑ)dη(ϑ)φ(ζ )dζ , (16) ψ(ζ

where η(ϑ) = η(ϑ , 0). Then A(0) and A∗ are adjoint operators. From the discussion in Appendix A, we know that ±iω∗ τ ∗ are eigenvalues of A(0). Thus, they are also eigenvalues of A∗ . We first need to compute the eigenvectors of A(0) and A∗ corresponding to +iω∗ τ ∗ and −iω∗ τ ∗ respectively. Suppose that q(ϑ) = ∗ ∗ ∗ ∗ (1, q2 )T eiω τ ϑ , (ϑ ∈ [−1, 0]) and q∗ (s) = K1¯ (1, q∗2 )eiω τ s , (s ∈ [0, 1]) ∗ are the eigenvectors of A(0) and A corresponding to the eigenvalues iω∗ τ ∗ and −iω∗ τ ∗ , respectively where,

q2 =

θ rx∗ [akh(1 − c) − (1 + αξ ) − 2ah(1 − c)x∗ ] k[d(1 + αξ + a(1 − c)hx∗ ) − θ ξ )] iω∗ θ [1 + αξ + a(1 − c)hx∗ ] , − d[1 + αξ + a(1 − c)hx∗ ] − θ ξ x∗ [akh(1 − c) − (1 + αξ ) − 2ah(1 − c)x∗ ] θ [1 + αξ − a(1 − c)hξ ](k − x∗ )eiω∗ τ ∗ iω∗ k[1 + αξ + a(1 − c)hx∗ ] − , θ r[1 + αξ − a(1 − c)hξ ](k − x∗ )eiω∗ τ ∗

q∗2 = −

B. Sahoo, S. Poria / Mathematical Biosciences 261 (2015) 62–73

and

(1) (1) + W11 (−1)z¯z + W02 (−1)

K = 1 + q2 q¯2 ∗ + τ ∗ q¯2 ∗ e−iω τ

θ r[1 + αξ − a(1 − c)hξ ](k − x∗ ) × d . + q 2 k[1 + αξ + a(1 − c)hx∗ ] ∗



We define, z(t) = q∗ , ut , W (t, ϑ) = ut (ϑ) − 2Re{z(t)q(ϑ)}. (17) On the center manifold C0 , we have W (t, ϑ) = W (z(t), z¯ (t), ϑ)

z2 + W11 (ϑ)z¯z 2 z3 z¯ 2 + W02 (ϑ) W30 (ϑ) + · · · 2 6

where, W (z, z¯ , ϑ) = W20 (ϑ)

z˙ (t) = iω∗ τ ∗ z + q¯ ∗ (0)f (0, W (z, z¯ , 0) + 2Re{zq(ϑ)}) =def iω∗ τ ∗ z + q¯ ∗ (0)f0 (z, z¯ ). iω∗ τ ∗ z(t) + g(z, z¯ ),

where

z2 z2 z¯ z¯ 2 + · · · (19) g(z, z¯ ) = q¯ ∗ (0)f0 (z, z¯ ) = g20 + g11 z¯z + g02 + g21 2 2 2 It follows from (17) and (18) that

ut (ϑ) = W (t, ϑ) + 2Re{z(t)q(t)}

∗ ∗ ∗ ∗ + (1, q2 )T eiω τ ϑ z + (1, q¯ 2 )T e−iω τ ϑ z¯ + · · ·

(20)

It follows that

⎜− τ = (1, q¯ ∗2 ) ⎜ ⎜ ⎝ K¯ That is,

⎞ ru21t (0) a(1 − c)u1t (0)u2t (0) − ⎟ k 1 + αξ + a(1 − c)hu1t (0) ⎟ ⎟. ⎠ θ a(1 − c)u1t (−1)u2t (−1) 1 + αξ + a(1 − c)hu1t (−1) ⎛

− ⎜ ⎜ g(z, z¯ ) = 1, q¯ ∗2 ⎜ ⎝ K¯

τ∗ 

where

P=

a(1 − c) . 1 + αξ

g20

2τ ∗ = K¯

g11 =

2τ ∗ K¯

g02 =

2τ ∗ K¯

g21 =

 

∗ ∗ q2 a(1 − c) r q2 q¯2 ∗ θ a(1 − c)e−2iω τ + , − + k 1 + αξ 1 + αξ  

Re(q2 )a(1 − c) Re(q2 )q¯2 ∗ θ a(1 − c) r + + , − k 1 + αξ 1 + αξ 

 ∗ ∗ q¯2 a(1 − c) r q¯2 q¯2 ∗ θ a(1 − c)e2iω τ + − , + k 1 + αξ 1 + αξ  a(1 − c) r (1) (1) (0) + 2W20 (0) − − 4W11 k 1 + αξ

2τ ∗ K¯  (2) (2) (1) (1) × 2W11 (0) + W20 (0) + q¯2 W20 (0) + 2q2 W11 (0)  2ha(1 − c) q¯ ∗ θ a(1 − c) (2q¯2 + q2 ) + 2 1 + αξ 1 + αξ   ∗ ∗ ∗ ∗ ( 1 ) ( 2 ) × 2e−iω τ q2 W11 (−1) + W11 (−1) + eiω τ

  2ah(1 − c)  ∗ ∗ (1) (2) × q¯2 W20 (−1) + W20 (−1) − (2q2 + q¯2 )e−iω τ . 1 + αξ From (15) and (17), we have

W˙ = u˙ t − z˙ − z¯˙ q¯ =

  τ∗ z2 z¯ 2 (1) (1) (1) g(z, z¯ ) = − (0) + W11 (0)z¯z + W02 (0) + · · · z + z¯ + W20 2 2 K¯

r z2 z¯ 2 (1) (1) (1) × z + z¯ + W20 (0) + W11 (0)z¯z + W02 (0) + · · · k 2 2

2 z z¯ 2 (2) (2) (2) + P q2 z + q¯2 z¯ + W20 (0) + W11 (0)z¯z + W02 (0) + · · · 2 2 

2 z z¯ 2 (1) (1) (1) × 1 − Ph z + z¯ + W20 (0) + W11 (0)z¯z + W02 (0) + · · · 2 2 

2 ¯ z z ( 1 ) ( 1 ) + ··· − ··· + P 2 h2 z2 + z¯ 2 + 2z¯z + 2W11 (0)z2 z¯ + 2W20 2

2 τ∗ ∗ z ∗ ∗ ∗ ∗ (1) + (−1) q¯2 θ P × ze−iω τ + z¯ eiω τ + W20 2 K¯

⎧ ⎨AW − 2Re{q¯ ∗ (0)f0 q(ϑ)},

ϑ ∈ [−1, 0),

⎩AW − 2Re{q¯ ∗ (0)f q(ϑ)} + f , 0 0

ϑ = 0.

=def AW + H(z, z¯ , ϑ),



ru21t (0) Pu1t (0)u2t (0) − ⎟ k 1 + Phu1t (0) ⎟ ⎟, θ Pu1t (−1)u2t (−1) ⎠ 1 + Phu1t (−1)

Now g(z, z¯ ) can be written in the form

Comparing the coefficients with (19), we have



z2 z¯ 2 = W20 (ϑ) + W11 (ϑ)z¯z + W02 (ϑ) 2 2





z2 ∗ ∗ ∗ ∗ (2) (−1) zq2 e−iω τ + z¯ q¯2 eiω τ + W20 2 2 ¯ z (2) (2) (−1)z¯z + W02 (−1) + · · · + W11 2 

z2 ∗ ∗ ∗ ∗ (1) (1) × 1 − Ph ze−iω τ + z¯ eiω τ + W20 (−1) + W11 (−1)z¯z 2  z¯ 2 ∗ ∗ ∗ ∗ (1) (−1) + · · · + P2 h2 z2 e−2iω τ + z¯ 2 e2iω τ 2z¯z + W02 2  z2 z¯  ∗ ∗ ∗ ∗ (1) (1) 4W11 + (−1)e−iω τ + 2W20 (−1)e−iω τ − · · · . 2

(18)

z and z¯ are local co-ordinate for central manifold C0 in the direction of q∗ and q¯ ∗ . Note that W is real if ut is real. We only consider real solutions. For solution ut ∈ C0 of (17), since μ = 0, we have

g(z, z¯ ) = q¯ ∗ (0)f (0, ut ) ⎛

z¯ 2 + ··· 2

×

In the remaining part of this section, we use the same notations as in Gopalsamy and He [59]. We first compute the coordinates to describe the center manifold C0 at μ = 0.

We rewrite this equation as z˙ (t) =

71

where, H(z, z¯ , ϑ) = H20 (ϑ)

z2 z¯ 2 + H11 (ϑ)z¯z + H02 (ϑ) + · · · 2 2

(21) (22)

Substituting the corresponding series into (21) and comparing the coefficients, we obtain

(A − 2iω∗ τ ∗ )W20 (ϑ) = −H20 (ϑ), AW11 (ϑ) = −H11 (ϑ).

(23)

From (21), we know that for ϑ ∈ [−1, 0)

H(z, z¯ , ϑ) = −q¯ ∗ (0)f0 q(ϑ) − q∗ (0)f¯0 (ϑ) = −g(z, z¯ )q(ϑ) − g¯ (z, z¯ )q¯ (ϑ).

(24)

Comparing the coefficients with (22), we get

H20 (ϑ) = −g20 q(ϑ) − g¯ 20 (ϑ)

(25)

and H11 (ϑ) = −g11 q(ϑ) − g¯ 11 (ϑ).

(26)

From (23) and (25) and the definition of A, it follows that

W˙ 20 (ϑ) = 2iω∗ τ ∗ W20 (ϑ) + g20 q(ϑ) + g¯ 02 q¯ (ϑ). ∗ ∗ Note that q(ϑ) = (1, q2 )T eiω τ ϑ , hence

W20 (ϑ) =

ig20

ω∗ τ ∗

ig¯ 02 ∗ ∗ ∗ ∗ ∗ ∗ q(0)eiω τ ϑ + q¯ (0)e−iω τ ϑ + E1 e2iω τ ϑ . 3ω ∗ τ ∗

Similarly, from (23) and (26) we obtain

72

B. Sahoo, S. Poria / Mathematical Biosciences 261 (2015) 62–73

W11 (ϑ) = −

ig11

ω∗ τ ∗

ig¯ 11 ∗ ∗ ∗ ∗ q(0)eiω τ ϑ + ∗ ∗ q¯ (0)e−iω τ ϑ + E2 ,

ωτ

where E1 = (E1(1), E1(2)) ∈ R2 and E2 = (E2(1), E2(2)) ∈ R2 are constant vectors, computed as: ⎛ rx∗ [ahk(1 − c) − (1 + αξ ) − 2ah(1 − c)x∗ ] ∗ ⎜2iω − k[1 + αξ + a(1 − c)hx∗ ] ⎜ E1 = 2 ⎜ ⎝ θ r[1 + αξ − a(1 − c)hξ ](k − x∗ ) −2iω∗ τ ∗ e k[1 + αξ + a(1 − c)hx∗ ] ⎛

q2 a(1 − c) r ⎜ − k − 1 + αξ ⎜ ×⎜ ⎝ q2 θ a(1 − c)e−2iω∗ τ ∗

d

θ



ξ

⎞−1

1 + αξ + a(1 − c)hx∗ ⎟ ⎟ ⎟ ⎠ ∗ ∗ 2iω∗ + d − de−2iω τ

⎞ ⎟ ⎟ ⎟ ⎠

1 + αξ

and ⎛ ⎜− ⎜ E2 = 2 ⎜ ⎝

rx∗ [ahk(1 − c) − (1 + αξ ) − 2ah(1 − c)x∗ ] k[1 + αξ + a(1 − c)hx∗ ] −

θ r[1 + αξ − a(1 − c)hξ ](k − x∗ ) k[1 + αξ + a(1 − c)hx∗ ]



r Re(q2 )a(1 − c) ⎜− k − 1 + αξ ⎜ ×⎜ ⎝ θ a(1 − c)Re(q2 )

d

θ



ξ

⎞−1

1 + αξ + a(1 − c)hx∗ ⎟ ⎟ ⎟ ⎠ 0

⎞ ⎟ ⎟ ⎟. ⎠

1 + αξ

Consequently, gij can be expressed by the parameters and delay τ ∗ . Thus we can compute the following values:

|g02 |2 g21 , g11 g20 − 2|g11 |2 − + 2ω τ 3 2 Re(C1 (0)) =− , Re(λ (τ ∗ )) = 2Re(C1 (0)), i

C1 (0) =

μ2 β2

T2 = −





Im(C1 (0)) + μ2 Im(λ (τ ∗ ))

ω∗ τ ∗

.

These expressions give a description of the bifurcating periodic solutions in the center manifold of system (3) at critical values τ ∗ . References [1] M.W. Sabelis, P.C.J. Van Rijn, When does alternative food promote biological pest control?, IOBC WPRS Bull. 29 (2006) 195. [2] G.R. Huxel, K. McCann, Food web stability: the influence of trophic flows across habitats, Am. Natural. 152 (3) (1998) 460–469. [3] G.R. Huxel, K. McCann, G.A. Polis, Effects of partitioning allochthonous and autochthonous resources on food web stability, Ecol. Res. 17 (4) (2002) 419–432. [4] B. Sahoo, S. Poria, Disease control in a food chain model supplying alternative food, Appl. Math. Model. 37 (8) (2013) 5653–5663. [5] B. Sahoo, S. Poria, The chaos and control of a food chain model supplying additional food to top-predator, Chaos, Solitons Fractals 58 (2014) 52–64. [6] P.D.N. Srinivasu, B.S.R.V. Prasad, M. Venkatesulu, Biological control through provision of additional food to predators: a theoretical study, Theor. Popul. Biol. 72 (1) (2007) 111–120. [7] P.D.N. Srinivasu, B.S.R.V. Prasad, Role of quantity of additional food to predators as a control in predator–prey systems with relevance to pest management and biological conservation, Bull. Math. Biol. 73 (10) (2011) 2249–2276. [8] B. Sahoo, Disease control through provision of alternative food to predator: a model based study, Int. J. Dyn. Control (2014), doi: 10.1007/s40435-014-0099-0. [9] M. Haque, D. Greenhalgh, When a predator avoids infected prey: a model-based theoretical study, Math. Med. Biol. 27 (1) (2010) 75–94. [10] L.N. Guin, M. Haque, P.K. Mandal, The spatial patterns through diffusion-driven instability in a predator–prey model, Appl. Math. Model. 36 (5) (2012) 1825–1841. [11] R.J. Orth, The importance of sediment stability in seagrass communities, Ecol. Mar. Benthos 6 (1977) 281–300. [12] A.W. Stoner, Species-specific predation on amphipod crustacea by the pinfish Lagodon rhomboides: mediation by macrophyte standing crop, Mar. Biol. 55 (1979) 201–207. [13] J.F. Savino, R.A. Stein, Predator-prey interaction between largemouth bass and bluegills as influenced by simulated, submersed vegetation, Trans. Am. Fisheries Soc. 111 (1982) 255–266. [14] J.F. Savino, R.A. Stein, Behavioural interactions between fish predators and their prey: effects of plant density, Anim. Behav. 37 (1989) 311–321. [15] O. Anderson, Optimal foraging by largemouth bass in structured environments, Ecology 65 (1984) 851–861.

[16] C.H. Ryer, Pipeflsh foraging: effects of fish size, prey size and altered habitat complexity, Mar. Ecol. Prog. Ser. 48 (1988) 37–45. [17] S.C. Pennings, Predator-prey interactions in opisthobranch gastropods: effects of prey body size and habitat complexity, Mar. Ecol. Prog. Ser. 62 (1990) 95–101. [18] J.H. Grabowski, Habitat complexity disrupts predator-prey interactions but not the trophic cascade on oyster reefs, Ecology, 85 (4) (2004) 995–1004. [19] L.S. Luckinbill, Coexistence in laboratory populations of Paramecium aurelia and its predator Didinium nasutum, Ecology 54 (6) (1973) 1320–1327. [20] N. MacDonald, Biological Delay Systems: Linear Stability Theory, Cambridge University Press, 1989. [21] Y. Song, S. Yuan, J. Zhang, Bifurcation analysis in the delayed Leslie–Gower predator–prey system, Appl. Math. Model. 33 (11) (2009) 4049–4061. [22] X. Yan, C. Zhang, Asymptotic stability of positive equilibrium solution for a delayed prey–predator diffusion system, Appl. Math. Model. 34 (1) (2010) 184–199. [23] X. Yang, L. Chen, J. Chen, Permanence and positive periodic solution for the single-species nonautonomous delay diffusive models, Comput. Math. Appl. 32 (4) (1996) 109–116. [24] M. Haque, S. Sarwardi, S. Preston, E. Venturino, Effect of delay in a Lotka–Volterra type predator–prey model with a transmissible disease in the predator species, Math. Biosci. 234 (1) (2011) 47–57. [25] G. Liu, J. Yan, Existence of positive periodic solutions for neutral delay Gause-type predator–prey system, Appl. Math. Model. 35 (12) (2011) 5741–5750. [26] Z. Du, Y. Lv, Permanence and almost periodic solution of a Lotka–Volterra model with mutual interference and time delays, Appl. Math. Model. 37 (3) (2013) 1054– 1068. [27] Y. Chen, F. Zhang, Dynamics of a delayed predator–prey model with predator migration, Appl. Math. Model. 37 (3) (2013) 1400–1412. [28] S. Sarwardi, M. Haque, P.K. Mandal, Ratio-dependent predator-prey model of interacting population with delay effect, Nonlinear Dyn. 69 (2012) 817–836. [29] J. Caperon, Time lag in population growth response of isochrysis galbana to a variable nitrate environment, Ecology 50 (1969) 188–192. [30] P.A. Abrams, L.R. Ginzburg, The nature of predation: prey dependent, ratio dependent or neither?, Trends Ecol. Evol. 15 (8) (2000) 337–341. [31] C.S. Holling, Some characteristics of simple types of predation and parasitism, Can. Entomol. 91 (7) (1959) 385–398. [32] I.J. Winfield, The influence of simulated aquatic macrophytes on the zooplankton consumption rate of juvenile roach, Rutilus rutilus, rudd, Scardinius erythrophthalmus, and perch, Perca fluviatilis, J. Fish Biol. 29 (1986) 37–48. [33] M. Kot, Elements of Mathematical Ecology, Cambridge University Press, (2001). [34] L.A. Berndt, S.D. Wratten, P. Hassan, Effects of buckwheat flowers on leafroller (Lepidoptera: Tortricidae) parasitoids in a New Zealand vineyard, Agricul. Forest Entomol. 4 (2002) 39–45. [35] S. Wratten, L. Berndt, G. Gurr, J. Tylianakis, P. Fernando, R. Didham, Adding floral diversity to enhance parasitoid fitness and efficacy, in: Proceedings of the First International Symposium on Biological Control of Arthropods, Honolulu, Hawaii, 2002, pp. 211–214. [36] G.M. Gurr, S.D. Wratten, FORUM’Integrated biological control’: a proposal for enhancing success in biological control, Int. J. Pest Manag. 45 (2) (1999) 81–84. [37] D. Perdikis, D. Lykouressis, L. Economou, The influence of temperature, photoperiod and plant type on the predation rate of Macrolophus pygmaeus on Myzus persicae, BioControl 44 (1999) 281–289. [38] D. Perdikis, D. Lykouressis, Macrolophus pygmaeus (Hemiptera: Miridae) population parameters and biological characteristics when feeding on eggplant and tomato without prey, J. Econ. Entomol. 97 (2004) 1291–1298. [39] J.T. Margaritopoulos, J.A. Tsitsipis, D.C. Perdikis, Biological characteristics of the mirids Macrolophus costalis and Macrolophus pygmaeus preying on the tobacco form of Myzus persicae (Hemiptera: Aphididae), Bull. Entomol. Res. 93 (2003) 39–45. [40] B. Vandekerkhove, P. De Clercq, Pollen as an alternative or supplementary food for the mirid predator Macrolophus pygmaeus, Biol. Control 53 (2010) 238–242. [41] P. De Clercq, Culture and natural enemies on factitious foods and artificial diets, Encyclopedia of Entomology I, 2008, pp. 1133–1136. [42] J.D. Harwood, K.D. Sunderland, W.O.C. Symondson, Prey selection by linyphiid spiders: molecular tracking of the effects of alternative prey on rates of aphid consumption in the field, Mol. Ecol. 13 (2004) 3549–3560. [43] J.D. Harwood, J.J. Obrycki, The role of alternative prey in sustaining predator populations, in: Proceedings of Second International Symposium on Biological Control of Arthropods II, 2005, pp. 453–462. [44] S. Toft, The quality of aphids as food for generalist predators: implications for natural control of aphids, Europ. J. Entomol. 102 (3) (2005) 371–383. [45] B.D. Hassard, Theory and Applications of Hopf Bifurcation, vol. 41, CUP Archive, 1981. [46] G.W. Harrison, Comparing predator-prey models to Luckinbill’s experiment with Didinium and Paramecium, Ecology 76 (2) (1995) 357–374. [47] G.W. Salt, Predator and prey densities as controls of the rate of capture by the predator Didinium nasutum, Ecology 55 (1974) 434–439. [48] E. Reükauf, Zur Biologie von Didinium nasutum Stein, J. Comp. Physiol. A: Neuroethol. Sens. Neural Behav. Physiol. 11 (4) (1930) 689–701. [49] C. Jost, S.P. Ellner, Testing for predator dependence in predator-prey dynamics: a non-parametric approach, Proc. R. Soc. Lond. B: Biol. Sci. 267 (2000) 1611–1620. [50] H.M. Butzel, A.B. Bolten, The relationship of the nutritive state of the prey organism Paramecium aurelia to the growth and encystment of Didinium nasutum, J. Eukaryot. Microbiol. 15 (2) (1968) 256–258. [51] D. Pauly, V. Christensen, J. Dalsgaard, R. Froese, F. Torres, Fishing down marine food webs, Science 279 (1998) 860–863.

B. Sahoo, S. Poria / Mathematical Biosciences 261 (2015) 62–73 [52] J.E. Duffy, Biodiversity loss, trophic skew and ecosystem functioning, Ecol. Lett. 6 (2003) 680–687. [53] R.A. Myers, B. Worm, Rapid worldwide depletion of predatory fish communities, Nature 423 (2003) 280–283. [54] B.S. Halpern, E.T. Borer, E.W. Seabloom, J.B. Shurin, Predator effects on herbivore and plant stability, Ecol. Lett. 8 (2005) 189–194. [55] A. Dobson, D. Lodge, J. Alder, G.S. Cumming, J. Keymer, J. McGlade et al., Habitat loss, trophic collapse, and the decline of ecosystem services, Ecology 87 (2006) 1915–1924.

73

[56] B.G. Veilleux, An Analysis of the predatory interaction between paramecium and didinium, J. Anim. Ecol. 48 (1979) 787–803. [57] R.V. Culshaw, S. Ruan, A delay-differential equation model of HIV infection of CD4+ T-cells, Math. Biosci. 165 (1) (2000) 27–39. [58] J.A. Dieudonné, Foundations of Modern Analysis, Academic Press, New York, 1960, p. 286. [59] K. Gopalsamy, X. He, Delay-independent stability in bidirectional associative memory networks, IEEE Trans. Neural Netw. 5 (6) (1994) 998–1002.

Effects of additional food in a delayed predator-prey model.

We examine the effects of supplying additional food to predator in a gestation delay induced predator-prey system with habitat complexity. Additional ...
1MB Sizes 1 Downloads 6 Views