Biocontrol in an impulsive predator–prey model

MBS 7526

No. of Pages 14, Model 5G

5 September 2014 Mathematical Biosciences xxx (2014) xxx–xxx 1

Contents lists available at ScienceDirect

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

Biocontrol in an impulsive predator–prey model

3 6

Q1

7

Division of Mathematics, University of Dundee, Dundee DD1 4HN, UK

8 9 10

a r t i c l e

1 3 2 3 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32

Alan J. Terry

i n f o

Article history: Received 24 October 2013 Received in revised form 17 August 2014 Accepted 18 August 2014 Available online xxxx AMS subject classifications (2010): 34K12 34K45 34K60 92D40

Q2

Keywords: Biological pest control Predator–prey model Impulsive releases Pest-eradication solution Economic threshold Beddington–DeAngelis functional response Holling type II functional response

a b s t r a c t We study a model for biological pest control (or ‘‘biocontrol’’) in which a pest population is controlled by a program of periodic releases of a fixed yield of predators that prey on the pest. Releases are represented as impulsive increases in the predator population. Between releases, predator-pest dynamics evolve according to a predator–prey model with some fairly general properties: the pest population grows logistically in the absence of predation; the predator functional response is either of Beddington–DeAngelis type or Holling type II; the predator per capita birth rate is bounded above by a constant multiple of the predator functional response; and the predator per capita death rate is allowed to be decreasing in the predator functional response and increasing in the predator population, though the special case in which it is constant is permitted too. We prove that, when the predator functional response is of Beddington–DeAngelis type and the predators are not sufficiently voracious, then the biocontrol program will fail to reduce the pest population below a particular economic threshold, regardless of the frequency or yield of the releases. We prove also that our model possesses a pest-eradication solution, which is both locally and globally stable provided that predators are sufficiently voracious and that releases occur sufficiently often. We establish, curiously, that the pest-eradication solution can be locally stable whilst not being globally stable, the upshot of which is that, if we delay a biocontrol response to a new pest invasion, then this can change the outcome of the response from pest eradication to pest persistence. Finally, we state a number of specific examples for our model, and, for one of these examples, we corroborate parts of our analysis by numerical simulations. Ó 2014 Elsevier Inc. All rights reserved.

34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55

56 57

1. Introduction

58

Biological pest control (or ‘‘biocontrol’’) is the control of a pest population by the release or conservation of biological agents [26]. These biological agents can also be called ‘‘biocontrol agents’’. Biocontrol agents include pathogens, parasitoids, and predators. Pathogens can be released directly or indirectly. Indirect release can be achieved by releasing, into a healthy wild pest population, a population of the pest that has been bred and infected with a pathogen in a laboratory [43]. The infected pests then spread disease into the healthy wild population. This strategy can work in controlling pests because, when a healthy pest becomes infected, it may no longer be able to function as a pest and may be more likely to die [42]. Regarding the use of parasitoids and predators as biocontrol agents, we offer examples: in greenhouses growing tomatoes and cucumbers, the parasitoid Encarsia formosa has long been used to control the whitefly pest Trialeurodes vaporariorum [6,23], and in the California citrus industry, the predatory vedalia

59 60 61 62 63 64 65 66 67 68 69 70 71 72 73

E-mail address: [email protected]

ladybird beetle Rodolia cardinalis has long been used to control the cottony-cushion scale pest Icerya purchasi [7,15]. Despite these and other successes, biocontrol can have disastrous consequences when it is poorly planned, an example being the introduction of the cane toad in Australia [24]. An informative review of biocontrol and its history is given in [22]. If a biocontrol agent can be quickly and cheaply produced, and if it can be released without causing disruption to the target ecosystem other than by disrupting the ability of individuals in the target pest population to act as pests, then the periodic release of the agent amounts to a sustainable control program. In such a situation, it is natural to ask how often the releases should occur and how much of the agent should be dispersed on each release? Of course, the answers depend on the goal of the program, which could be complete eradication of the target pest population, or the long-term reduction of this population, or its temporary reduction. Mathematical modelling has a role to play in determining the frequency and yield of releases of a biocontrol agent in particular circumstances. Hitherto there have been modelling studies looking at: the periodic release of predators [17,27,23,18,21,20]; the

http://dx.doi.org/10.1016/j.mbs.2014.08.009 0025-5564/Ó 2014 Elsevier Inc. All rights reserved.

Q1 Please cite this article in press as: A.J. Terry, Biocontrol in an impulsive predator–prey model, Math. Biosci. (2014), http://dx.doi.org/10.1016/ j.mbs.2014.08.009

74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94

MBS 7526

No. of Pages 14, Model 5G

5 September 2014 Q1 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141

2

A.J. Terry / Mathematical Biosciences xxx (2014) xxx–xxx

periodic release of pests infected by a disease [42]; the periodic release of predators and infected pests [28]; the periodic release of infected pests combined with periodic applications of pesticides [43]; and the periodic release of predators combined with periodic applications of pesticides [44]. The release of a biocontrol agent may be an event of short duration, especially in comparison to the time that passes between consecutive events of this kind. A similar comment holds if the events are applications of pesticides. Therefore, it is typical in modelling studies, including all of those cited in the previous paragraph, to represent the releases of the biocontrol agent or the applications of pesticides as impulsive events, with populations being governed by systems of differential equations between consecutive events. The resulting impulsive differential equation models can be studied using Floquet theory and comparison methods, thereby making it possible to find conditions on the period and strengths of the events so as to drive the pest extinct as time tends to infinity (see the papers cited in the previous paragraph for examples of such calculations). In view of the observations above, our aim in this paper is to study an impulsive predator–prey model, where the prey is assumed to be a pest and where predators are released periodically as a biocontrol program to suppress the prey (pest) population. Although we are not the first to consider such a problem (for example, see [17,27,23,18,21,20]), our model contains novel features. In particular, the predator dynamics are treated in quite a general way. More specifically, the per capita birth and death rates for the predator are represented as functions that satisfy a small number of biologically feasible properties, and the predator functional response is assumed to take either the Holling type II form or the Beddington–DeAngelis form. We outline the rest of this paper. In Section 2, we state our model and comment on basic features of its solutions. In Section 3, we prove that, under certain conditions, the pest population will persist, potentially at a level sufficient to cause economic harm, regardless of the frequency or yield of the predator release events. In Section 4, we show that our model possesses a pest-eradication solution and we establish conditions for the local and global stability of this solution. In Section 5, we prove, curiously, that the pesteradication solution can be locally stable when it is not globally stable, and we comment on the practical implications of this result. Specifically, we mention the implication that, if action against a new pest invasion is delayed, then we may miss our window of opportunity to control it. In Section 6, we show how our work has a general applicability by stating various examples for our model. For one of these examples, we corroborate parts of our analysis by numerical simulations in Section 7. In Section 8, we offer concluding remarks.

142

2. The model

143

Before we state the model to be studied, we describe the assumptions on which it is based. We begin with the following set of assumptions:

144 145 146 147 148 149 150 151 152 153 154 155 156

(A1) N ¼ NðtÞ is the population of a pest at time t P 0. (A2) The pest has a predator, where P ¼ PðtÞ is the population of this predator at time t P 0. (A3) We assume that, initially (time t ¼ 0), the pest population is positive, that is, Nð0Þ > 0. (A4) To control the pest, a fixed positive quantity Q of the predator is released at times nT, where n 2 f0 [ Zþ g and T is a positive constant. The first release is modelled by assuming that the initial predator population satisfies Pð0Þ P Q ; the condition Pð0Þ > Q therefore indicates that the predator is already naturally present when the first release occurs.

Subsequent releases are modelled as impulsive increases in the predator population. We shall refer to the releases as predator release events or release events. We shall refer to Q as the release yield and to T as the inter-release period. (A5) Between predator release events, the predator and pest populations evolve according to a predator–prey model, which is described below. The predator–prey model mentioned in assumption (A5) is: dN dt dP dt

  ¼ rN 1  NK  PFðN; PÞ; ¼ P fBðFðN; PÞ; PÞ  ½DðFðN; PÞÞ þ mPg;

158 159 160 161 162 163 164

165

) ð1Þ

167

where we now describe the terms in this model. The term rNð1  ðN=KÞÞ represents logistic growth for the pest population in the absence of predation (r and K are positive constants, respectively known as the intrinsic growth rate and environmental carrying capacity of the pest). The term FðN; PÞ is the predator functional response, that is, the prey capture rate per predator. We will suppose that: ) FðN; PÞ ¼ 1þb1aN ; Nþb2 P

168 169 170 171 172 173 174

175

a and b1 are positive constants; and b2 is a non-negative constant: ð2Þ

177

If b2 > 0, then the predator functional response FðN; PÞ is of Beddington–DeAngelis type [4,9,44], whilst if b2 ¼ 0, then the predator functional response is of Holling type II [31]. These are both common forms for predator functional responses in predator–prey models. Indeed, if there can be said to be a standard form, then it would arguably be Holling type II [31]. The term BðFðN; PÞ; PÞ in (1) is the per capita birth rate for the predator (that is, the reproduction rate of the total predator population, divided by the predator population), for which we make some general assumptions: 9 BðFðN; PÞ; PÞ is a continuous autonomous non-negative function of > > > > > FðN; PÞ P 0 and P P 0:It has continuous first order partial > = derivatives in N and P: It is increasing ðnot necessarily strictlyÞ in > > > FðN; PÞ: We assume that there is a positive constant v such that > > > ; BðFðN; PÞ; PÞ 6 vFðN; PÞ for FðN; PÞ P 0 and P P 0:

178

ð3Þ

190

The assumptions in (3) are made in order to help ensure that solutions to our biocontrol model will exist and have desirable properties, and also in order to reflect certain requirements that are biologically necessary, sensible, or reasonable. Thus, it is biologically necessary for the per capita birth rate for the predator, BðFðN; PÞ; PÞ, to be non-negative. It is sensible to expect BðFðN; PÞ; PÞ to be increasing in the predator functional response FðN; PÞ, since the latter represents the prey capture rate and predators are more likely to be able to reproduce successfully if they are able to catch a sufficient quantity of prey. It is sensible to allow BðFðN; PÞ; PÞ to not necessarily be strictly increasing in the predator functional response or prey capture rate FðN; PÞ, since for all sufficiently small values of the prey capture rate, predators may not have sufficient energy resources to reproduce successfully [5,39]. Finally, it is sensible to expect BðFðN; PÞ; PÞ to be limited by the prey capture rate FðN; PÞ (by being bounded above by vFðN; PÞ for some positive constant v), since prey consumption provides the energy required for predator reproduction. Note that, from the assumptions in (3), we have 0 6 BðFðN; PÞ; PÞ 6 vFðN; PÞ for FðN; PÞ P 0 and P P 0, from which we must have BðFðN; PÞ; PÞ ¼ 0 when FðN; PÞ ¼ 0, that is, an individual predator will not reproduce when its prey capture rate is zero, as we might expect. The per capita birth rate for the predator may depend on the predator population P not only through the functional response FðN; PÞ but also directly, which is why we have used the notation

191

Q1 Please cite this article in press as: A.J. Terry, Biocontrol in an impulsive predator–prey model, Math. Biosci. (2014), http://dx.doi.org/10.1016/ j.mbs.2014.08.009

157

179 180 181 182 183 184 185 186 187

188

192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214

MBS 7526

No. of Pages 14, Model 5G

5 September 2014 Q1 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236

237

A.J. Terry / Mathematical Biosciences xxx (2014) xxx–xxx

BðFðN; PÞ; PÞ to represent it. Direct dependence on P may arise, through mating interactions, when the predator is a sexually reproducing species [38]. Specific forms for BðFðN; PÞ; PÞ that contain direct dependence on P are stated in Section 6 below. The possibility that the per capita birth rate for the predator is directly dependent on P is permitted in our model but is not essential to it. We may ignore this direct dependence, by writing BðFðN; PÞ; PÞ simply as BðFðN; PÞÞ, without invalidating the assumptions in (3) or invalidating our analysis that depends on (3). The assumption of autonomy in (3) (that is, no explicit dependence on time t) is reasonable if predator reproduction is determined only by the prey capture rate and also possibly by mating interactions, that is, seasonal fluctuations are negligible. This could occur for species of predator that may reproduce at any time of year. The term DðFðN; PÞÞ þ mP in (1) is the per capita death rate for the predator. Here m is a non-negative constant. We allow either m ¼ 0 or m > 0 for the sake of generality. When m ¼ 0, then the term mP is zero and is irrelevant. When m > 0, then the term mP increases in the predator population and represents mortality that may arise through aggressive encounters between predators [2,12,38]. For the function DðFðN; PÞÞ, we make some general assumptions:

the model, it could also be written in terms of a delta function (for examples of impulsive differential equation models written in terms of delta functions, see [29,30]). In view of our assumptions above, model (5) has four qualitatively different cases:

For model (5), a unique solution exists for t P 0, and it satisfies positivity (NðtÞ P 0 and PðtÞ P 0 for t P 0) and upper boundedness (there exist positive constants, N and P, such that NðtÞ 6 N and PðtÞ 6 P for t P 0). Moreover, NðtÞ and PðtÞ satisfy strict positivity, that is, NðtÞ > 0 and PðtÞ > 0 for t P 0. All of this can be established by standard results (for example, see Section 3 in [44], and see also the books by Lakshmikantham et al. [16] and Bainov and Simeonov [3]). Three special cases of model (5) have been studied in [21]. Models similar to model (5), possessing fairly general properties but not completely overlapping model (5), have been studied in [20,17,18].

275

9 > > > > FðN; PÞ P 0: It has continuous first order partial derivatives in N > > = and P: We assume Dð0Þ ¼ l > 0 and that DðFðN; PÞÞ is decreasing > > ðnot necessarily strictlyÞ in FðN; PÞ: Also DðFðN; PÞÞ P m > 0 for > > > > ; FðN; PÞ P 0; for a constant m with 0 < m 6 l:

3. Biocontrol may fail

286

The Economic Injury Level (EIL) of a pest population can be defined as the population size or density at which the economic damage caused by the pest equals the cost of managing the population [32]. In practice, we do not wish for a pest population to even approach its EIL, and so action is generally taken when it reaches a smaller level called its Economic Threshold (ET) [33,11]. (Note, however, that with some types of pest invasion, action is taken to eradicate the pest as soon as the invasion has been detected. In such cases, we could define the ET to be zero.) The biocontrol program in model (5) is stronger if the release yield Q is larger (more predators dispersed with each release event) or if the inter-release period T is smaller (more frequent release events). The following theorem shows that, under some circumstances, the biocontrol program can fail to reduce the pest population below an economic threshold regardless of how large Q is, or how small T is.

287

Theorem 1. In the biocontrol model in (5), suppose a < b2 r. Also suppose that the pest has an economic threshold U N . If   U N < K 1  ba2 r , then the biocontrol program will ultimately fail,

303

that is, NðtÞ > U N for all time t sufficiently large.

306

Proof. First note, in view of the assumption that a < b2 r, that we must have b2 > 0. Using that b2 > 0, and using positivity for NðtÞ and strict positivity for PðtÞ, we have:

307

DðFðN; PÞÞ is a continuous autonomous positive function of

239

ð4Þ

240

As with the assumptions in (3), the assumptions in (4) are made in order to help ensure that solutions to our biocontrol model will exist and have desirable properties, and also in order to reflect certain requirements that are biologically necessary, sensible, or reasonable. Thus, it is biologically sensible for the per capita death rate for the predator, DðFðN; PÞÞ þ mP, to be positive, and (4) ensures this. It is sensible to expect the per capita death rate for the predator to be decreasing in the predator functional response FðN; PÞ, because, the more food a predator has, the less likely it will be to die from lack of food [39]. It is sensible to expect the per capita death rate for the predator to not necessarily always be strictly decreasing in the predator functional response FðN; PÞ, because, for instance, increases in the predator functional response above a certain level will not change the fact that a predator has enough to eat [39]. The assumption of autonomy in (4) is reasonable if predator mortality is not significantly influenced by seasonal fluctuations. An example of such a fluctuation would be where predators are significantly more susceptible to death during winter. We treat any such fluctuations as negligible. Collecting together assumptions (A1) to (A5), and (1)–(4), we are now in a position to state our biocontrol model:

241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260

261

9 9   dN > > > > = ¼ rN 1  NK  PFðN; PÞ; > > dt > t > 0; t – nT; > > > > dP > > > > ¼ P fBðFðN; PÞ; PÞ  ½DðFðN; PÞÞ þ mPg; ; > > dt > =   NðnTÞ ¼ NðnT Þ; t ¼ nT;  > > PðnTÞ ¼ PðnT Þ þ Q ; > > > > > Nð0Þ > 0; Pð0Þ P Q ; > > > > > r; K; m are constants where r > 0; K > 0; m P 0; > > ; F; B; and D satisfy; respectively; ð2Þ; ð3Þ; and ð4Þ;

263

ð5Þ

264

where n 2 Zþ , and where nT  denotes the time instantaneously before the time nT. Note that, due to the presence of impulses in

265

3

(i) (ii) (iii) (iv)

m > 0; m > 0; m ¼ 0; m ¼ 0;

b2 b2 b2 b2

> 0; ¼ 0; > 0; ¼ 0.

267 268 269 270 271 272 273 274

aNP aNP aN 6 ¼ fort P 0: 1 þ b1 N þ b2 P b2 P b2

ð6Þ

Using (2), (6), and the first equation in model (5), we have, for t > 0, except at the times of impulsive releases of the predator:

   

dN N aN a N   : P rN 1  ¼ rN 1  dt K b2 b2 r K

277 278 279 280 281 282 283 284 285

288 289 290 291 292 293 294 295 296 297 298 299 300 301 302

304 305

308 309

310

312 313 314

315

317

But then, bearing in mind that NðtÞ is not changed by each of the impulsive release events and that Nð0Þ > 0 from (5), we have NðtÞ P N 1 ðtÞ for t P 0, where N 1 ð0Þ ¼ Nð0Þ > 0 and where, for t > 0,

 

dN1 a N1  : ¼ rN1 1  b2 r dt K

276

ð7Þ 318 319 320

321

ð8Þ

Q1 Please cite this article in press as: A.J. Terry, Biocontrol in an impulsive predator–prey model, Math. Biosci. (2014), http://dx.doi.org/10.1016/ j.mbs.2014.08.009

266

323

MBS 7526

No. of Pages 14, Model 5G

5 September 2014 Q1

4

A.J. Terry / Mathematical Biosciences xxx (2014) xxx–xxx

327

Using the assumption that a < b2 r, and either by solving for N1 ðtÞ directly or by a standard phase portrait argument, we deduce from   (8) that N 1 ðtÞ ! K 1  ba2 r > 0 as t ! 1. Hence, since NðtÞ P N1 ðtÞ   for t P 0, and since U N < K 1  ba2 r by assumption, then NðtÞ > U N

328

for all time t sufficiently large.

324 325 326

329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367

h

We note in the first sentence of the proof of Theorem 1 that we must have b2 > 0. It follows that the predator functional response is of Beddington–DeAngelis type in the theorem. Consequently, an upper bound can be placed on the total predation rate that depends only on the pest population, as we see in (6). This limits the impact of predation on the pest population, regardless of how quickly predators are released. Combining the limited impact of predation with the condition a < b2 r, which limits the efficiency of the prey capture rate per predator, causes the failure of the biocontrol program in the theorem. The idea that the periodic release of a predator may fail to control a pest when the predator is subject to a Beddington–DeAngelis functional response has also been observed in [21,20]. It is debatable whether the total predation rate would be limited in real life in the way that it is in (6). Nevertheless, the Beddington–DeAngelis functional response compares well with other types of functional response when fitted to real-world data sets [31]. Moreover, biocontrol programs, in the real world, have not always met with the successes that have been anticipated or hoped for [15], a fact that is certainly consistent with Theorem 1. In any case, Theorem 1 does not hold when the functional response is of Holling type II, as happens when b2 ¼ 0. Suppose that model (5) is considered an appropriate representation for a particular real-world application of a biocontrol program, but the available data suggests that the program will fail because the conditions of Theorem 1 will hold. Then it may be possible to avoid failure of the program by changing the biological conditions that hold, specifically by selectively breeding the predators, in a controlled environment, until individuals with increased voracity (such that a > b2 r) have appeared [8]. These individuals could then be bred until they have produced a population that is large enough to use in a biocontrol program. We will see in the next section that, when a > b2 r, then a sufficiently strong biocontrol program will always eradicate the pest. Remark 1. Under the conditions of Theorem 1, the pest population remains, in the long run, above the economic threshold U N . However, we have found by numerical simulation that the biocontrol program may nevertheless maintain the pest population somewhat below its environmental carrying capacity K (results not shown).

368

4. Successful biocontrol

369

At the start of the previous section, we mentioned the concept of an economic threshold as a maximum tolerable level for a pest population. It is certainly possible to prevent a pest population from exceeding an economic threshold if a biocontrol program eradicates it. Therefore, in this section, we seek criteria for pest eradication in model (5). To this end, we begin by seeking to understand the dynamics of the model when the pest has already been eradicated. Thus, setting NðtÞ  0, and using (2)–(4), we see that the first two equations in model (5) decouple and we are led to consider the properties of the submodel

370 371 372 373 374 375 376 377 378

379

dP dt

¼ ðl þ mPÞP;

t > 0;

PðnTÞ ¼ PðnT  Þ þ Q ; 381

Pð0Þ P Q;

382

þ

9 t – nT; > =

t ¼ nT;

> ;

By a T-periodic solution to submodel (9), we mean a solution which is periodic, with period T. The following lemma establishes that submodel (9) has an attracting T-periodic solution.

383

Lemma 1. For t P 0, submodel (9) has a positive T-periodic solution e PðtÞ such that

386

e ¼ PðtÞ

385

387

388

 lt

lP e

l þ mP ð1  elt Þ

;

t 2 ½0; TÞ;

ð10Þ

390 391

e PðTÞ ¼ P ;

ð11Þ

where

P ¼

393 394

8 > <

Q 1elT

395

if m ¼ 0; rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi n offi

1 2 4Qml l þ ð l Þ þ ; Qm  Qm   l T 2m 1e

> :

;

if m > 0: ð12Þ

397

e Moreover, for every solution PðtÞ of submodel (9), jPðtÞ  PðtÞj ! 0 as t ! 1.

398

Proof. We first consider the case m > 0. Thus, assume m > 0. Then, we begin by solving (9) on the interval ½nT; ðn þ 1ÞT for any n 2 f0 [ Zþ g, that is, we solve:

400

dP ¼ ðl þ mPÞP; dt

t 2 ðnT; ðn þ 1ÞTÞ;

399

401 402

403

ð13Þ

405

with initial condition PðnTÞ and final condition Pððn þ 1ÞTÞ ¼ Pððn þ 1ÞT  Þ þ Q . Solving on ½nT; ðn þ 1ÞTÞ, by separation of variables and partial fractions, we find:

406 407 408

409

lðtnTÞ

PðtÞ ¼

lPðnTÞe

;

l þ mPðnTÞð1  elðtnTÞ Þ

t 2 ½nT; ðn þ 1ÞTÞ:

ð14Þ



Hence, from Pððn þ 1ÞTÞ ¼ Pððn þ 1ÞT Þ þ Q , we have:

411 412

413

lT

Pððn þ 1ÞTÞ ¼ Q þ

lPðnTÞe

l þ mPðnTÞð1  elT Þ

:

ð15Þ

Setting Pn ¼ PðnTÞ, then (15) is equivalent to a one-dimensional map:

Pnþ1 ¼ f ðPn Þ ¼ Q þ

lPn elT l þ mPn ð1  elT Þ

:

ð16Þ

415 416 417

418 420

b of the map must satisfy P b ¼ f ð PÞ. b Setting Fixed points P b in (16), and solving the resulting equation for P, b Pn ¼ P nþ1 ¼ P yields a unique positive fixed point P , which is found to satisfy (12) (where m > 0). Since P is a positive fixed point of (16), then it follows from (14) and (16) that (10) and (11) is a T-periodic solution to (9), assuming Pð0Þ ¼ P  . (Note by (16) that we must have P > Q , a fact which we mention here because it will be needed in the proof of Theorem 3 below.) Finally, it is easy to show, by a simple induction argument or by a standard cobwebbing argument, that the orbit of P0 ¼ Pð0Þ P Q under the map in (16) will always tend to P  . But then, in view of (14) and (16), it follows trivially by an -argument that e jPðtÞ  PðtÞj ! 0 as t ! 1, where PðtÞ is any solution of submodel

421

(9). This completes the proof for the case m > 0. For the case m ¼ 0, we may use a very similar argument as that used for the case m > 0. h

435

^ pital’s rule, it is easy to see that Remark 2. By L’Ho

438

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi " #) 1 4Qml Q 2 Qm  l þ ¼ ðQm  lÞ þ : 2m 1  elT 1  elT

(

ð9Þ lim

m!0þ

where n 2 Z . Q1 Please cite this article in press as: A.J. Terry, Biocontrol in an impulsive predator–prey model, Math. Biosci. (2014), http://dx.doi.org/10.1016/ j.mbs.2014.08.009

384

422 423 424 425 426 427 428 429 430 431 432 433 434 436 437

439

441

MBS 7526

No. of Pages 14, Model 5G

5 September 2014 Q1 442 443 444 445 446 447 448 449 450

451

453 454

455

457

5

A.J. Terry / Mathematical Biosciences xxx (2014) xxx–xxx

e to submodel (9) in the case Therefore, in Lemma 1, the solution PðtÞ m > 0 tends, in the limit as m ! 0þ , to the solution to submodel (9) in the case m ¼ 0. From our observations so far in this section, model (5) has a peste e is defined in Lemma 1. We where PðtÞ eradication solution, ð0; PðtÞÞ, now prove a theorem on the local and global stability of this pesteradication solution. e is Theorem 2. In model (5), the pest-eradication solution, ð0; PðtÞÞ, locally stable (LS) provided that

Z

T

0

e a PðtÞ dt > rT: e 1 þ b2 PðtÞ

ð17Þ

e is globally asymptotically stable (GAS) provided that Further, ð0; PðtÞÞ

Z 0

e a PðtÞ

T

e 1 þ b1 K þ b2 PðtÞ

dt > rT:

ð18Þ

can be used to examine the local stability of fixed points in linear systems with periodic coefficients [14]. To apply this theory to the linear system in (20), we first construct the fundamental matrix UðtÞ of this linear system that is obtained from assuming identity matrix initial conditions and from solving over the time interval 0 6 t 6 T [45]. Thus, for 0 6 t 6 T, we have:

459 460

Proof. The proof is the same in each of the following four cases: (i) m > 0; b2 > 0; (ii) m > 0; b2 ¼ 0; (iii) m ¼ 0; b2 > 0; (iv) m ¼ 0; b2 ¼ 0.

462

465

NðtÞ ¼ xðtÞ;

466

in model (5), where xðtÞ and yðtÞ are small perturbations. By expanding each of the first two equations in model (5) as a Taylor e series about ð0; PðtÞÞ, using (2)–(4), and neglecting quadratic and higher order terms in xðtÞ and yðtÞ, and also by linearising the impulsive perturbation conditions in model (5), we find that the linearisation of model (5) about the pest-eradication solution is:

463

467 468 469 470 471

472

e þ yðtÞ; PðtÞ ¼ PðtÞ

 ¼ xðtÞ r 

 ;

ð19Þ

9 > > =

dxðtÞ ae P ðtÞ dt P ðtÞ 1þb2 e   > t > 0; dyðtÞ > e ¼ xðtÞWðtÞ  yðtÞ l þ 2m PðtÞ ;; dt  xðnTÞ ¼ xðnT  Þ; t ¼ nT; yðnTÞ ¼ yðnT  Þ; 475

478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498

9 > > > > > t – nT; > > = > > > > > > > ; ð20Þ

474 476

where n 2 Zþ , and where

x1 ðtÞ x2 ðtÞ ; UðtÞ ¼ y1 ðtÞ y2 ðtÞ

ð22Þ

where ðx1 ðtÞ; y1 ðtÞÞ and ðx2 ðtÞ; y2 ðtÞÞ are solutions, still yet to be derived, of the linear system in (20) with the following initial conditions:

502 503 504

507 508 509 510

511

x1 ð0Þ ¼ 1;

x2 ð0Þ ¼ 0;

y1 ð0Þ ¼ 0;

y2 ð0Þ ¼ 1:

513

Solving, we obtain:

514

  @ dP

WðtÞ ¼ :

@N dt ðN;PÞ¼ð0;e P ðtÞÞ

ð21Þ

We do not consider the linearisation of the initial conditions in model (5) because they are not used in determining local stability in what follows. The coefficient of xðtÞ in the second equation in (20), namely WðtÞ, has been left in the unevaluated form in (21), since it cannot be evaluated explicitly without knowing the explicit forms for the functions defined in (3) and (4). Note, however, that it does exist, in view of the assumptions in (2)–(4). Moreover, it is not needed in an explicit form in what follows. e is T-periodic. Therefore, the coefficient of Now the function PðtÞ xðtÞ in the first equation in (20), and the coefficient of yðtÞ in the second equation in (20), are both also T-periodic. Although the coefficient of xðtÞ in the second equation in (20) has not been evaluated in an explicit form, we may nevertheless deduce from the form that it takes in (21), and from the assumptions in (2)–(4), including the assumption of autonomy in (3) and (4), that it is T-periodic as well. Thus, all of the coefficients of the linear system in (20) are T-periodic. e The pest-eradication solution, ð0; PðtÞÞ, of model (5) will be locally stable if the fixed point (or equilibrium) ðx ; y Þ ¼ ð0; 0Þ of the linearised model in (20) is locally stable [14]. Floquet theory

ð23Þ

 Z th i  e du for t 2 ½0; T; l þ 2m PðuÞ y2 ðtÞ ¼ exp  0 " # ! Z t e a PðuÞ du for t 2 ½0; T: r x1 ðtÞ ¼ exp e 0 1 þ b2 PðuÞ

ð24Þ ð25Þ 517

There is no need to calculate the exact form of y1 ðtÞ since it is not required in the analysis that follows. Next, we calculate the ‘‘Floquet multipliers’’, which are defined as the eigenvalues of the ‘‘monodromy matrix’’ UðTÞ [45], that is, they are the solutions of the following eigenvalue problem:

detðUðTÞ  kIÞ ¼ 0;

ð26Þ

where I is the 2  2 identity matrix. From (22), (23), and (26), we have:

det

x1 ðTÞ  k

0

y1 ðTÞ

y2 ðTÞ  k

¼ 0:

ð27Þ

Thus, from (27), there are two Floquet multipliers, k1 and k2 , which (using (25) and (24)) we find to be:

k1 ¼ x1 ðTÞ ¼ exp

Z

"

#

 Z k2 ¼ y2 ðTÞ ¼ exp 

T

h

1 :¼ exp

0

T

521 522

523 525 526 527

530 531 532

535 536

 e l þ 2m PðuÞ du : i

ð29Þ

e According to Floquet theory, the pest-eradication solution, ð0; PðtÞÞ, is locally stable if the absolute values of all Floquet multipliers are less than unity. Since we clearly have jk2 j < 1 from (29), local stability therefore depends on whether jk1 j < 1, which holds provided that (17) holds. This completes the local stability part of the proof. e is Now we prove that the pest-eradication solution, ð0; PðtÞÞ, globally asymptotically stable provided that condition (18) is satisfied. Thus, assume (18) holds. Then for any sufficiently small  > 0, we can write:

(

520

ð28Þ

0

Z

519

533

!

e a PðuÞ du ; r e 1 þ b2 PðuÞ

T

518

528



0



) !

e  Þ að PðtÞ dt r e  Þ 1 þ b1 ðK þ Þ þ b2 ð PðtÞ

< 1:

538 539 540 541 542 543 544 545 546 547

548

ð30Þ 550

e e Next note that PðtÞ is T-periodic and PðtÞ > 0 for t 2 ½0; T (see

551

e in Lemma 1). Hence for all  > 0 small enough, the definition of PðtÞ e we have PðtÞ   > 0 for t 2 ½nT; ðn þ 1ÞT for n 2 f0 [ Zþ g. Given

552

this, and bearing in mind that (30) holds for all sufficiently small  > 0, we now choose  > 0 so that (30) holds and so that e   > 0 for t 2 ½nT; ðn þ 1ÞT for n 2 f0 [ Zþ g. PðtÞ

554

Q1 Please cite this article in press as: A.J. Terry, Biocontrol in an impulsive predator–prey model, Math. Biosci. (2014), http://dx.doi.org/10.1016/ j.mbs.2014.08.009

501

515

To study the local stability of the pest-eradication solution, e ð0; PðtÞÞ, we linearise about it by setting

461

500

505



x2 ðtÞ ¼ 0; 458

499

553 555 556

MBS 7526

No. of Pages 14, Model 5G

5 September 2014 Q1

6

A.J. Terry / Mathematical Biosciences xxx (2014) xxx–xxx

608 557 558 559 560 561 562 563 564 565

566

From model (5), observe by positivity that dN 6 rNð1  ðN=KÞÞ dt for t > 0, except at the times of the predator release events. Bearing in mind that NðtÞ is not changed by each release event, we then have by a standard comparison argument that NðtÞ 6 N 1 ðtÞ for 1 t P 0, where N 1 ð0Þ ¼ Nð0Þ > 0 and where dN dt ¼ rN 1 ð1  ðN 1 =KÞÞ for t > 0. Clearly, N 1 ðtÞ ! K as t ! 1. We may conclude that NðtÞ < K þ  for all t large enough – say, for t P t 1 . Using positivity, we then have 0 6 NðtÞ < K þ  for t P t 1 . Next, observe from model (5) and positivity that

9 > t – nT; > =

dP P ðl þ mP ÞP; t > 0; dt PðnTÞ ¼ PðnT  Þ þ Q ; t ¼ nT; Pð0Þ P Q;

568 569 570 571 572 573 574 575 576 577

578

e   for all t large enough – say, for t P t  . that PðtÞ > PðtÞ 2  Let t3 be the least multiple of T that is greater than maxft 1 ; t2 g. Then t3 ¼ n T for some integer n > 0. Also, from our argument so e   > 0 for far, we have both 0 6 NðtÞ < K þ  and PðtÞ > PðtÞ tP

Hence, by (5), we have:

"

e  Þ dNðtÞ að PðtÞ 6 NðtÞ r  e  Þ dt 1 þ b1 ðK þ Þ þ b2 ð PðtÞ

9 > >  ; t > t 3 ; t – nT; > > =

#

> > > > ;

NðnTÞ ¼ NðnT  Þ; 06

Nðt 3 Þ

< K þ ;

580

ð32Þ

581

where n > n . Integrating (32) for t 2 ½nT; ðn þ 1ÞT, for any particular n P n , we obtain:

582

583

NðtÞ 6 NðnTÞ

Z t( r  exp nT

585 586

587

590 591 592 593 594 595

596

601 602 603 604 605 606 607

ð34Þ  Nðt 3 Þ ðnn Þ

From (34), we have, by induction, that NðnTÞ 6 1 for n P n . Therefore, since 0 6 Nðt 3 Þ < K þ  by (32), and since 1 < 1 by (30), we have NðnTÞ ! 0 as n ! 1. e We already know by our choice for  > 0 that PðuÞ   > 0 for u 2 ½nT; ðn þ 1ÞT, where n P n . Then, for t 2 ½nT; ðn þ 1ÞT for n P n , we have

Z t( r nT

600

ð33Þ

¼ NðnTÞ1:

) Z ðnþ1ÞT e  Þ að PðuÞ du 6 rdu ¼ rT: e nT  Þ 1 þ b1 ðK þ Þ þ b2 ð PðuÞ ð35Þ

598 599

) ! e  Þ að PðuÞ du : e  Þ 1 þ b1 ðK þ Þ þ b2 ð PðuÞ

e From (33), and using the T-periodicity of PðtÞ and (30), we have: ) ! Z ðnþ1ÞT ( e  Þ að PðuÞ Nððn þ 1ÞTÞ 6 NðnTÞ exp r du e nT  Þ 1 þ b1 ðK þ Þ þ b2 ð PðuÞ

589

t 2 ½0; TÞ;

ð36Þ

By (33) and (35), and the fact that NðtÞ satisfies strict positivity (as mentioned in Section 2), we have 0 < NðtÞ 6 NðnTÞ expðrTÞ for t 2 ½nT; ðn þ 1ÞT, for n P n . Then, since we have seen that NðnTÞ ! 0 as n ! 1, we must have NðtÞ ! 0 as t ! 1. e as t ! 1. This is equivalent to Next, we prove that PðtÞ ! PðtÞ e  1 < PðtÞ < PðtÞ e þ 1 proving that, for any 1 > 0, we have PðtÞ for all sufficiently large t. Thus, we first choose any 1 > 0. Now suppose that A is a positive constant, and define, for t P 0, e 2 ðtÞ, where, for t 2 ½0; T, we have: the T-periodic function P

ð37Þ

where

P2 ¼

613 614

8 > > > < > > > :

615

Q ; if m ¼ 0; 1  eAT sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi " ffi# 1 4QmA 2 Qm  A þ ðQm  AÞ þ ; 2m 1  eAT

if m > 0: ð38Þ

617

e 2 ðtÞ is equivalent to PðtÞ, e e is Note that, when A ¼ l, then P where PðtÞ defined in Lemma 1. In fact, for A sufficiently close to l, then

618

e 2 ðtÞ < PðtÞ e þ P

1 2

for t P 0:

ð39Þ

We may choose A > 0 sufficiently close to l to ensure that (39) holds if we let A ¼ DðFð0 ; 0ÞÞ  vFð0 ; 0Þ for 0 > 0 sufficiently small. This is because, by the assumptions in (2)–(4), DðFðN; 0ÞÞ  vFðN; 0Þ is continuous in N, and DðFð0; 0ÞÞ  vFð0; 0Þ ¼ Dð0Þ  0 ¼ l > 0. Thus, we let A ¼ DðFð0 ; 0ÞÞ  vFð0 ; 0Þ where 0 > 0 is chosen to be sufficiently small to ensure that A > 0 and that (39) holds. Since NðtÞ ! 0 as t ! 1, there exists an integer n2 > 0 such that NðtÞ < 0 for t P n2 T :¼ t4 . Combining this with the positivity of N, we see that 0

0 6 NðtÞ <  for t P

t 4

¼

n2 T:

ð40Þ

Using (40), the positivity of P, and the assumptions in (2)–(4), we find:

619

620 622 623 624 625 626 627 628 629 630 631

632 634 635 636

637

0

BðFðN; PÞ; PÞ  DðFðN; PÞÞ 6 vFðN; 0Þ  DðFðN; 0ÞÞ < vFð ; 0Þ  DðFð0 ; 0ÞÞ ¼ A < 0 for t P t 4 ¼ n2 T: ð41Þ

639

Using the first equation in (31), (41), and the positivity of P, we have

640 641

642

dP ðl þ mPÞP 6 6 ðA þ mPÞP for t 2 ðnT; ðn þ 1ÞTÞ; n P n2 : dt ð42Þ

644

From (5), (42), the strict positivity and upper boundedness (by P) of P at t ¼ t 4 (strict positivity and upper boundedness are mentioned in Section 2), and from trivially adapting Lemma 3.4 in [44] to shift the time origin, we obtain P 1 ðtÞ 6 PðtÞ 6 P 2 ðtÞ where P 1 ðtÞ is the solution to

645

dP1 ¼ ðl þ mP1 ÞP1 ; t > t 4 ¼ n2 T; dt P1 ðnTÞ ¼ P1 ðnT  Þ þ Q ; t ¼ nT; 0<

P1 ðt 4 Þ

where n >

¼

n2 ,

Pðt 4 Þ

0<

> > ;

6 P;

¼

Pðt 4 Þ

6 P;

646 647 648 649

650

ð43Þ 652

and where P2 ðtÞ is the solution to

dP2 ¼ ðA þ mP2 ÞP2 ; t > t 4 ¼ n2 T; dt P2 ðnTÞ ¼ P2 ðnT  Þ þ Q ; t ¼ nT; P2 ðt 4 Þ

9 > t – nT; > =

653

9 > t – nT; > = > > ;

654

ð44Þ 656

where, again, n > n2 . Trivially modifying the argument in the proof of Lemma1, we e find from (43) that P1 ðtÞ ! PðtÞ as t ! 1. Since A > 0, we can also

657

trivially modify the argument in the proof of Lemma 1 to find, from e 2 ðtÞ as t ! 1, where P e 2 ðtÞ is the T-periodic (44), that P 2 ðtÞ ! P

660

function of time t P 0 that satisfies (36)–(38) where A ¼ DðFð0 ; 0ÞÞ  vFð0 ; 0Þ. Collecting our observations, and in

662

Q1 Please cite this article in press as: A.J. Terry, Biocontrol in an impulsive predator–prey model, Math. Biosci. (2014), http://dx.doi.org/10.1016/ j.mbs.2014.08.009

610 611

e 2 ðTÞ ¼ P ; P 2

ð31Þ

> > ;

where n 2 Zþ . Hence by standard comparison results for impulsive differential equations [16,3], we have PðtÞ P P 1 ðtÞ for t P 0, where P1 ðtÞ is the solution to (9). But by Lemma 1, we see that e jP 1 ðtÞ  PðtÞj ! 0 as t ! 1. Combining our observations, we deduce

t3 .

AP 2 eAt e 2 ðtÞ ¼ P ; A þ mP2 ð1  eAt Þ

658 659

661 663

MBS 7526

No. of Pages 14, Model 5G

5 September 2014 Q1

7

A.J. Terry / Mathematical Biosciences xxx (2014) xxx–xxx

724 665

particular using (39), we deduce that, for all t sufficiently large, e  1 < PðtÞ < P e 2 ðtÞ þ ð1 =2Þ < PðtÞ e þ 1 . This completes then PðtÞ

666

the proof.

667

682

Note that the local stability condition for pest eradication in (17) is independent of the environmental carrying capacity K of the pest, whilst the global stability condition in (18) is not. However, the global stability condition becomes an increasingly close approximation to the local stability condition when K becomes smaller, which is sensible because, reasoning heuristically, the smaller K is, the smaller the pest population will be, so the closer it will be to eradication, so the more likely it is that local stability of the pest-eradication solution will ensure its global stability. (Global stability of the pest-eradication solution necessarily implies its local stability, of course.) We derived conditions for pest eradication in Theorem 2 by studying the local and global stability of the pest-eradication solution. A more direct approach to deriving conditions for pest eradication is simply to ask for dN=dt < 0 for all sufficiently large time t, which we do in the following theorem:

683

Theorem 3. In model (5), assume a > b2 r. Let W ¼ P   Q , where P

684

Wb1 K is given in (12). Let I ¼ 1þb22b and J ¼ K½ðabb21rÞWr . Then the pestr 1

664

668 669 670 671 672 673 674 675 676 677 678 679 680 681

685 686

687 689

h

e is globally asymptotically stable (GAS) if eradication solution, ð0; PðtÞÞ, either 2

J  I > 0;

690

or

693

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi J  I2 6 0 and ðI2  JÞ  I < 0:

691

694 695 696 697

698 700

ð45Þ

ð46Þ

Proof. The proof is the same in each of the following four cases: (i) m > 0; b2 > 0; (ii) m > 0; b2 ¼ 0; (iii) m ¼ 0; b2 > 0; (iv) m ¼ 0; b2 ¼ 0. 

First, recall from the proof of Lemma 1 that P > Q . Hence: 

W ¼ P  Q > 0:

ð47Þ

703

Next, let I and J be as defined in the statement of the theorem. We prove that, if either (45) or (46) holds, then, for N 2 ½0; K, we have

706

ðN þ IÞ2 þ ðJ  I2 Þ > 0:

707

Clearly (48) holds for N 2 ½0; K if J  I2 > 0, that is, if (45) holds. Now suppose J  I2 6 0 in (48). Then I2  J P 0 and we can

701 702

704

708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723

ð48Þ

rearrange (48) to ðN þ IÞ2 > I2  J, knowing that the right hand side of this latter inequality is greater than or equal to zero. Hence qffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffi N þ I > I2  J or N þ I <  I2  J . But N þ I > I2  J will hold qffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffi for N 2 ½0; K provided that I2  J  I < 0, whilst N þ I <  I2  J qffiffiffiffiffiffiffiffiffiffiffiffi will hold for N 2 ½0; K provided that  I2  J  I > K. Therefore, qffiffiffiffiffiffiffiffiffiffiffiffi (48) will hold for N 2 ½0; K if J  I2 6 0 and either I2  J  I < 0 or qffiffiffiffiffiffiffiffiffiffiffiffi  I2  J  I > K. The combined requirements that J  I2 6 0 and qffiffiffiffiffiffiffiffiffiffiffiffi I2  J  I < 0 are equivalent to (46). The combined requirements qffiffiffiffiffiffiffiffiffiffiffiffi that J  I2 6 0 and  I2  J  I > K cannot hold. To see why they qffiffiffiffiffiffiffiffiffiffiffiffi cannot hold, notice that they imply I2  J < K  I, which, with J  I2 6 0 or I2  J P 0, implies that K  I > 0 or K þ I < 0. The latter inequality implies that b2 W < b1 K  1, which is impossible since W > 0 by (47). Thus, we have established that, if either (45) or (46) holds, then (48) holds for N 2 ½0; K, that is,

ðN þ IÞ2 þ ðJ  I2 Þ > 0 for N 2 ½0; K:

ð49Þ

726

e is globWe now prove that the pest-eradication solution, ð0; PðtÞÞ, ally asymptotically stable (GAS) when (49) holds. Using the definitions of I and J, we find that (49) is equivalent to:

727



 N ð1 þ b1 N þ b2 WÞ  aW < 0 for N 2 ½0; K: r 1 K

ð50Þ

Also, since a > b2 r by assumption, then

  N  a 6 b2 r  a < 0 for N 2 ½0; K: b2 r 1  K

729

730 732 733

734

ð51Þ

736

From (50) and (51), the following quantity exists and is positive:

737

(  )  r 1  NK ð1 þ b1 N þ b2 WÞ  aW     : d :¼ min N2½0;K b2 r 1  NK  a

740

Given that d > 0, and using (47), the following quantity exists and is positive:

  1 minfW; dg: c :¼ 2

ð53Þ

Since W > 0 by (47), then:

W  c > 0:

ð55Þ

Using (54), it follows from (55) that HðNÞ < 0 for N P K. Then we will have HðNÞ < 0 for all N > 0 if we have HðNÞ < 0 for N 2 ð0; K. Using (54), we will have HðNÞ < 0 for N 2 ð0; K if, for N 2 ð0; K,

 

  N N a >r 1 ð1 þ b1 N þ b2 WÞ  aW: c b2 r 1  K K

741 742

743 745 746

ð54Þ

In view of (54), we may formulate a well-defined function HðNÞ of N > 0 as follows:

    N aNðW  cÞ  : HðNÞ ¼ rN 1  K 1 þ b1 N þ b2 ðW  cÞ

738

ð52Þ

ð56Þ

747 749 750 751

752 754 755 756 757

758 760

But (56) is true for N 2 ð0; K in view of how we have defined c in (53) and in view of (51). Thus, HðNÞ < 0 for all N > 0. Now we saw in the proof of Theorem 2, in the paragraph e containing (31), that PðtÞ P P1 ðtÞ for t P 0, where P 1 ðtÞ ! PðtÞ as

761

e is defined in Lemma 1. It is clear from the t ! 1, where PðtÞ e e that PðtÞ is least on ½0; T at time T  , that is, at the definition of PðtÞ time momentarily before T. From this, and from the T-periodicity of e e e e Þ ¼ PðtÞ and the fact that PðtÞ satisfies (9), we have PðtÞ P PðT

765

e PðTÞ  Q ¼ P  Q ¼ W > 0 for t P 0. Then, for c > 0 as defined in (53), we have PðtÞ P P1 ðtÞ > W  c for all t sufficiently large, say for t P t P 0. Using that PðtÞ > W  c for t P t  , and also using positivity, the first equation in (5), (2), and (54), we have, for t P t  , except at the times of the predator release events, that dN=dt 6 HðNÞ, where HðNÞ is defined in (55). Also, by the positivity and upper boundedness of N from the penultimate paragraph of Section 2, we have 0 6 Nðt Þ < N. Therefore, bearing in mind that NðtÞ is not changed by each release event, we may apply a standard comparison argument: NðtÞ 6 N 1 ðtÞ for t P t , where 0 6 N 1 ðt Þ ¼ Nðt Þ < N and where, for t > t , we have dN 1 =dt ¼ HðN 1 Þ. Now clearly Hð0Þ ¼ 0, and, from observations above, it is obvious that HðN 1 Þ < 0 for N 1 > 0. But then a standard phase portrait argument reveals that N 1 ðtÞ ! 0 as t ! 1. Hence, and also bearing in mind the positivity of N, we have 0 6 NðtÞ 6 N 1 ðtÞ ! 0 as t ! 1. In other words, NðtÞ ! 0 as t ! 1. e Given that NðtÞ ! 0 as t ! 1, we must have PðtÞ ! PðtÞ as t ! 1 by the argument in the last part of the proof of Theorem 2. In summary, we have shown that the pest-eradication solution, e ð0; PðtÞÞ, is globally asymptotically stable (GAS), as required. h

769

Q1 Please cite this article in press as: A.J. Terry, Biocontrol in an impulsive predator–prey model, Math. Biosci. (2014), http://dx.doi.org/10.1016/ j.mbs.2014.08.009

728

762 763 764

766 767 768

770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789

MBS 7526

No. of Pages 14, Model 5G

5 September 2014 Q1

8

A.J. Terry / Mathematical Biosciences xxx (2014) xxx–xxx

793

It is not obvious, at first inspection, if the global stability conditions for pest eradication in Theorems 2 and 45 can hold. We now present a corollary showing that these conditions can and do hold for all sufficiently strong biocontrol programs, provided a > b2 r.

794

Corollary 1. In model (5), assume a > b2 r.

790 791 792

(i) Suppose m ¼ 0. Then for all sufficiently large release yields Q > 0, or all sufficiently small inter-release periods T > 0: (A) condition (18) in Theorem 2 will hold; (B) condition (45) in Theorem 45 will hold if, additionally, b2 ¼ 0; (C) condition (46) in Theorem 45 will hold if, additionally, b2 > 0. (ii) Suppose m > 0. Then for all sufficiently small inter-release periods T > 0: (A) condition (18) in Theorem 2 will hold; (B) condition (45) in Theorem 45 will hold if, additionally, b2 ¼ 0; (C) condition (46) in Theorem 45 will hold if, additionally, b2 > 0.

795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811

812 813 814

In all six of these cases, the pest will be eradicated. Proof. Proof of (i). First let W ¼ P  Q as in Theorem 45, and recall that W > 0 from (47). By assumption in part (i), m ¼ 0. But then, using (12), we lT

815 816 817 818 819 820 821 822 823 824 825 826

827 829

830 831 832 833 834 835 836 837 838 839 840 841 842

843

Qe have W ¼ P  Q ¼ 1e lT . Clearly, then, for all sufficiently large release yields Q > 0, or all sufficiently small inter-release periods T > 0; W will be arbitrarily large. Hence, if the criteria for pest eradication can be shown to hold for all W sufficiently large, then they will hold for all sufficiently large release yields Q > 0, or all sufficiently small inter-release periods T > 0. Our task, then, is simply to show that the conditions in parts (A), (B), and (C) of part (i) of the corollary hold for all W sufficiently large. We begin by showing that the condition in part (A) of part (i) of the corollary (that is, condition (18)) holds for W sufficiently large. e Now we have seen in the proof of Theorem 45 that PðtÞ PW>0 for t P 0. Hence:

Z 0

e a PðtÞ

T

e 1 þ b1 K þ b2 PðtÞ

dt P

Z

T

0

aW aWT : dt ¼ 1 þ b1 K þ b2 W 1 þ b1 K þ b2 W

Hence (18) will hold provided that ðaWTÞ=ð1 þ b1 K þ b2 WÞ > rT, which, in view of the assumption in the corollary that a > b2 r, holds if W > ðr þ rb1 KÞ=ða  b2 rÞ. Thus, (18) holds for W sufficiently large, as required. Next we show that part (B) of part (i) of the corollary holds for W sufficiently large, that is, we show that (45) holds for W sufficiently large when b2 ¼ 0. To this end, recall the definitions of I and J in Theorem 45. Setting b2 ¼ 0 in these definitions reveals that (45) holds for W sufficiently large, as required. Finally we show that part (C) of part (i) of the corollary holds for W sufficiently large, that is, we show that (46) holds for W sufficiently large when b2 > 0. Thus, assume b2 > 0. Also, for reasons that will become clear, let us define GðWÞ as follows: 2 b2 rW 2

2

845

GðWÞ :¼

846

848

Since we have assumed that b2 > 0, then certainly b2 r > 0, and so GðWÞ P 0 for W sufficiently large, say for W P X where X > 0. qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Now the conditions in (46) are J  I2 6 0 and ðI2  JÞ  I < 0.

849

We argue that these conditions hold for W such that:

847

850

852

þ ½2b2 rð1 þ b1 KÞ  4ab1 K W þ rð1 þ b1 KÞ :

ð57Þ

2



b1 K  1 r W > max X; ; b2 a  b2 r

 :

ð58Þ

Notice that the lower bound on W in (58) is well-defined in view of the assumption in the corollary that a > b2 r, and is also positive. Thus, assume W satisfies (58). Then GðWÞ P 0 where GðWÞ is defined in (57). But GðWÞ P 0 is equivalent to J  I2 6 0, so the first condition in (46) holds. Next, from (58), we have W > r=ða  b2 rÞ, which is equivalent to J > 0 given the assumption that a > b2 r. From J > 0, we have

853

I2  J < I 2 :

ð59Þ

862

From (58), we have W > ðb1 K  1Þ=b2 , which is equivalent to I > 0. Since we have already established that J  I2 6 0, then we have

863

I2  J P 0, which, in combination with (59) and I > 0, leads to qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðI2  JÞ < I, or, equivalently, ðI2  JÞ  I < 0. But the latter

865

inequality is the second condition in (46). Therefore, we have shown that, for W large enough to satisfy (58), then the conditions in (46) hold when b2 > 0. This completes the proof of part (i). Proof of (ii). The argument is the same as in part (i), except that W will only necessarily be arbitrarily large for all sufficiently small interrelease periods T > 0, and not necessarily for all sufficiently large release yields Q > 0. In fact, when m > 0, then it is trivial to see, from (12), that W ¼ P  Q can be made arbitrarily large for all sufficiently small inter-release periods T > 0. However, for fixed T, then W may be bounded above for all Q > 0, and in particular W may not be large enough to satisfy the lower bounds imposed on it in part (i). To see this, observe, using (12) and the assumption m > 0, that

867

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  l 2 2mW þ l ¼ Q þ Q 2 þ þ 2ZQ ; m m

ð60Þ



   where Z ¼ lð1 þ elT Þ = mð1  elT Þ > 0. Since Q ; l=m, and Z are all positive, then, from (60), we have

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  l 2 l 2mW þ l Qþ < Q þ þZ ¼Zþ : m m m

ð61Þ

From (61), we find:

Z lð1 þ elT Þ W< ¼ : 2 2mð1  elT Þ

855 856 857 858 859

860

864

866 868 869 870 871 872 873 874 875 876 877 878 879 880

881

883 884 885

886

888 889

890

ð62Þ 892

Inequality (62) is entirely independent of the parameter K. If all model parameters except K and Q are given, with a > b2 r and m > 0, then K can be chosen large enough that the lower bounds on W in parts (A), (B), and (C) of part (i) all exceed the upper bound on W in (62). Hence, regardless of the value of Q > 0, none of these lower bounds on W can hold, and the arguments in part (i) fail. This does not in itself prove, given a > b2 r; m > 0, and a particular fixed T, that pest eradication is impossible for all sufficiently large release yields Q. We have merely shown that this cannot be established by the methods in part (i). h

893

Remark 3. In model (5), when a < b2 r, then by Theorem 1, pest eradication can never occur. However, when a > b2 r, then by Corollary 1, pest eradication will always occur for all sufficiently strong biocontrol programs – high release yield, or short interrelease period, in the case m ¼ 0; short inter-release period in the case m > 0. Note that a > b2 r holds automatically when b2 ¼ 0, that is, when the functional response is of Holling type II.

903

Next we present a corollary showing that the pest-eradication solution is locally stable for all sufficiently strong biocontrol programs, provided a > b2 r.

910

Corollary 2. In model (5), assume a > b2 r.

913

Q1 Please cite this article in press as: A.J. Terry, Biocontrol in an impulsive predator–prey model, Math. Biosci. (2014), http://dx.doi.org/10.1016/ j.mbs.2014.08.009

854

894 895 896 897 898 899 900 901 902

904 905 906 907 908 909

911 912

MBS 7526

No. of Pages 14, Model 5G

5 September 2014 Q1

(i) If m ¼ 0, then for all sufficiently large release yields Q > 0, or all sufficiently small inter-release periods T > 0, the local stability condition for pest eradication in (17) will hold. (ii) If m > 0, then for all sufficiently small inter-release periods T > 0, the local stability condition for pest eradication in (17) will hold.

914 915 916 917 918 919 920

921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952

953

Proof. The result is true by an argument very similar to that used, in the proof of Corollary 1, to prove that condition (18) in Theorem 2 will hold. h From Corollaries 1 and 2, we know that, when a > b2 r, then, for all sufficiently small inter-release periods T > 0: (i) the local stability condition for pest eradication in (17) will hold; (ii) the global stability condition for pest eradication in (18) will hold; (iii) the global stability condition for pest eradication in (45) will hold when b2 ¼ 0; and (iv) the global stability condition for pest eradication in (46) will hold when b2 > 0. It is consistent with intuition that pest eradication should be associated with sufficiently small inter-release periods T, because, when the release yield Q is held fixed and T is reduced, then the rate at which predators are released per unit time is increased. What happens when the inter-release period T is not ‘‘sufficiently small’’? This question is worth asking because, in the real-world, there will be a limit to how quickly releases can be carried out. For example: it takes time and resources for a new batch of predators to be bred and prepared for release, and people cannot spend all of their time dispersing predators because they will have other things to do as well. In fact, it is practical and economical to seek to minimise the frequency of predator release events, whilst still ensuring control of the pest population. Therefore, given any of the pest eradication conditions in (17), (18), (45), or (46), we are motivated to ask: for any particular release yield Q, is there a limit to how large we can choose the inter-release period T whilst still ensuring that the given condition holds? The answer in each case is: yes, there is such a limit. We prove this in the remainder of this section. Thus, we first consider condition (17). From Lemma 1, we have e PðtÞ > 0 for t 2 ½0; T and also P > 0. Then we may bound above the left hand side of (17) as follows:

Z 0

T

e a PðtÞ dt 6 e 1 þ b2 PðtÞ

957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972

Z

T

0

Z

<

955 956

9

A.J. Terry / Mathematical Biosciences xxx (2014) xxx–xxx

0

e ¼ a PðtÞdt

Z

T

0 T

aP elt dt ¼

 lt

alP e dt l þ mP ð1  elt Þ

  aP 

l

 1  elT :

ð63Þ

Thus, when (17) holds, then by (63) we have ððaP Þ=lÞ   1  elT > rT. But the latter inequality cannot hold for any T sufficiently large, because, from (12), we trivially have, for all T suffi  ciently large, that ððaP Þ=lÞ 1  elT < rT, whether m ¼ 0 or m > 0. Hence, condition (17) does not hold for T sufficiently large. By a very similar argument, condition (18) does not hold for T sufficiently large. Finally, we show that, if T is sufficiently large, then neither (45) nor (46) can hold. Indeed, as T ! 1, then from (12), we find P  ! Q , so W ¼ P   Q ! 0. Hence, for T sufficiently large, we have J < 0 where J is defined in Theorem 45. But then J  I2 < 0 where I is defined in Theorem 45. Hence (45) does not hold. Since qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi J  I2 < 0, then I2  J > 0, so that ðI2  JÞ > 0. Hence if I 6 0, then (46) cannot hold. If I > 0, then, since we have said J < 0, we must qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi have ðI2  JÞ > I2 ¼ I, so that ðI2  JÞ  I > 0, so that (46) cannot hold. In summary, neither (45) nor (46) can hold for T sufficiently large, as required.

Proving that our stability conditions for the pest-eradication solution do not hold for sufficiently large inter-release periods T is not the same as proving that pest eradication is impossible when the inter-release period T is large. After all, we have not proven that our stability conditions are the only such conditions. Nevertheless, it is consistent with intuition that there would be a limit to how large we can choose the inter-release period T whilst still ensuring pest eradication. Indeed, if the inter-release period T is very large, then we are basically doing nothing to control the pest the vast majority of the time. Under such circumstances, it would be unreasonable to expect pest eradication. From Corollaries 1 and 2, and our arguments after Corollary 2, we can say that, when a > b2 r, then our local and global stability conditions for the pest-eradication solution hold for sufficiently small values of the inter-release period T but do not hold for sufficiently large values of T. Bearing in mind our comment above that it is economical to seek to minimise the frequency of predator release events whilst still ensuring control of the pest population, this suggests a practical method for estimating a value for T in a real-world application. We begin with a very small value of T, so that the global stability conditions will hold; then we increase T in small increments until we reach a first value T 0 where none of these conditions still hold; and finally we choose T to be the value incrementally before T 0 . Thus, the largest value of T that ensures global stability from at least one of the conditions can be used in the application, assuming it is large enough to be practical.

973

5. Local stability does not imply global stability

999

The conditions for local and global stability of the pest-eradication solution in Theorem 2 are not the same. Also, it is difficult to say, at first glance, how the local stability condition in Theorem 2 is related to the global stability conditions for pest eradication in Theorem 45. The following theorem reveals that the pest can persist in model (5) even when the pest-eradication solution is locally stable.

1000

Theorem 4. In model (5), the following can be true simultaneously:

1007

(i) the local stability condition for the pest-eradication solution (condition (17)) holds; (ii) there exists a positive constant LN and solutions ðN; PÞ such that NðtÞ P LN for t P 0.

1008 1012

Proof. We begin by deriving conditions such that part (ii) in the statement of the theorem holds. To do this, we first assume that Nð0Þ 6 K. Then NðtÞ is bounded above by K for t P 0. This is easily seen by a simple bounding and comparison argument such as the one used in the paragraph before the paragraph containing (31). Using that NðtÞ 6 K for t P 0, and using positivity, then, between predator release events, we find, from model (5), that

1015



dP vaKP vaK 6  mP 6 P dt 1 þ b1 K þ b2 P 1 þ b1 K

for t 2 ½nT; ðn þ 1ÞTÞ;

976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998

1001 1002 1003 1004 1005 1006

1009 1013 1010 1014 1011

1016 1017 1018 1019 1020 1021

1024

n 2 f0 [ Zþ g;

1025

ð64Þ

where M ¼ ððvaKÞ=ð1 þ b1 KÞÞ  m. But then, applying the impulsive release condition at t ¼ ðn þ 1ÞT, we have

Pððn þ 1ÞTÞ ¼ Pððn þ 1ÞT  Þ þ Q 6 PðnTÞeMT þ Q

1026 1028 1029 1030

1031

for n 2 f0 [ Zþ g: ð65Þ

1033

Using (65), it is easily established by induction that PðnTÞ 6   Pð0ÞenMT þ Q 1 þ eMT þ . . . þ eðn1ÞMT for n 2 Zþ . Now Sn ¼ 1þ

1034

Q1 Please cite this article in press as: A.J. Terry, Biocontrol in an impulsive predator–prey model, Math. Biosci. (2014), http://dx.doi.org/10.1016/ j.mbs.2014.08.009

975

1022

 m :

Hence, by a standard comparison argument, we see that

PðtÞ 6 PðnTÞeMðtnTÞ

974

1035

MBS 7526

No. of Pages 14, Model 5G

5 September 2014 Q1 1036 1037 1038

10

A.J. Terry / Mathematical Biosciences xxx (2014) xxx–xxx

eMT þ . . . þ eðn1ÞMT is a geometric progression, with first term 1 and common ratio eMT . If M < 0, then 0 < eMT < 1 and Sn is bounded above by S1 ¼ 1=ð1  eMT Þ > 0. Therefore, assuming

1039 1041 1042 1043

1044 1046 1047

M¼

vaK 1 þ b1 K

 m < 0;

ð66Þ þ

we will have PðnTÞ 6 Pð0Þ þ Q =ð1  eMT Þ for n 2 f0 [ Z g. Using this, and again using (66), we then have, by (64), that:

PðtÞ 6 Pð0Þ þ Q =ð1  e

MT

Þ ¼: U P

for t P 0:

0 < LN < K

1051

and

1052

dN

> 0 for t > 0; t – nT; n 2 Zþ : dt N¼LN

1055 1056

1057

ð68Þ

ð69Þ

In view of (67) and positivity, and assuming (68) holds, we will have

(69) if dN > 0 for P 2 ½0; U P , which holds if dt N¼L

1059 1060

which holds if

1061

1063

  r 1  LKN ð1 þ b1 LN Þ Q    ; Pð0Þ þ < ¼ U P ð1  eMT Þ a  b2 r þ b2KrLN

1064

provided that

1065 1067

a > b2 r:

ð70Þ

ð71Þ

1071 1073

Nð0Þ 2 ½LN ; K;

1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089

ð72Þ

will exist and will satisfy NðtÞ P LN for t P 0, since (69) will hold and since NðtÞ is not changed at the times of the release events (that is, at the times nT where n 2 f0 [ Zþ g). In other words, we have derived conditions such that part (ii) in the statement of the theorem holds. Therefore, the theorem will be true if (17), (66), (68), (70), (71), and (72) can hold simultaneously. But we may easily satisfy these conditions simultaneously if we choose parameters and initial conditions according to the following procedure: (P1) Choose b2 P 0; r > 0; a > b2 r. (P2) Choose l > 0; m P 0; Q > 0; T > 0 such that (17) holds. This is possible for T small enough in view of Corollary 2. (P3) Choose K > 0. Choose LN so that 0 < LN < K, and choose m so that 0 < m 6 l. Choose Pð0Þ P Q ; v > 0. Choose b1 > 0 large enough to satisfy (66) and (70). This is clearly possible. (P4) Choose Nð0Þ to satisfy (72). h

1090

1091 1092 1093 1094 1095 1096 1097 1098

The assumptions in model (5) are sufficiently general to permit many special cases, or specific examples. Five examples are stated in Table 1. The examples for model (5) stated in the first four rows of Table 1 are such that:

1111

BðFðN; PÞ; PÞ ¼ vFðN; PÞ:

ð73Þ

In words, Eq. (73) says that the per capita birth rate for the predator, BðFðN; PÞ; PÞ, is proportional to the predator functional response, or prey capture rate, FðN; PÞ. This seems like a reasonable assumption when one bears in mind that predator reproduction is the conversion of prey biomass into new predators. In the model in the bottom row of Table 1, the per capita birth rate for the predator takes the form

 P : BðFðN; PÞ; PÞ ¼ vFðN; PÞ dþP

1070

1074

1110

1101 1102 1103 1104 1105 1106 1107 1108 1109

1112 1113 1114 1115

1116 1118

Remark 4. From Theorem 4, it is clear that the pest-eradication e can be locally stable whilst not being globally solution, ð0; PðtÞÞ, stable. This suggests that there may be circumstances in which delaying a response to a new pest invasion can change the outcome of the response from pest eradication to pest persistence – the window of opportunity to control the pest may be lost by the delay. Indeed, we demonstrate such circumstances by simulation in Fig. 3 below.

ð74Þ

1120 1121 1122 1123 1124 1125

1128

Notice that the per capita birth rate in (74) is equal to the per capita birth rate in (73), scaled by the term P=ðd þ PÞ where d is a positive constant. This scaling has little impact on the birth rate when the predator population P is large, but reduces the birth rate when P is small. Thus, the scaling reflects the following idea: if the predators reproduce sexually, then individuals may struggle to reproduce when the population is small, because the population density may also then be small, causing individuals to struggle to find a mate. Another form for the per capita birth rate for the predator that could be used in model (5), and that is qualitatively different to the forms in (73) and (74), is as follows:

1129

8 0; 0 6 FðN; PÞ 6 /1 ; > > < n   h io FðN;PÞ/1 2 p b cos BðFðN; PÞ; PÞ ¼ 1 þ /2 /1 ; /1 6 FðN; PÞ 6 /2 ; 2 > > : b; FðN; PÞ P /2 ;

1140

ð75Þ

1142

where /1 ; /2 , and b are all positive constants, with /1 < /2 . The form in (75) reflects several ideas:

1143

 when the predator functional response, or prey capture rate, FðN; PÞ, is sufficiently small (less than or equal to /1 ), then predators do not have a high enough level of energy intake to be able to reproduce;  when the prey capture rate FðN; PÞ reaches a sufficient level (more than /1 ), then predators become able to reproduce, and their rate of reproduction then grows with the prey capture rate since prey consumption fuels their reproduction;  when the prey capture rate FðN; PÞ reaches a high level (/2 ), then the per capita birth rate for predators reaches a plateau level, since there may be natural limits to how quickly individuals can reproduce which cannot be overcome simply by eating more.

1145

Q1 Please cite this article in press as: A.J. Terry, Biocontrol in an impulsive predator–prey model, Math. Biosci. (2014), http://dx.doi.org/10.1016/ j.mbs.2014.08.009

1119

1126



Combining our observations, we see that, if (66), (68), (70), and (71) hold, then any trajectory with initial conditions ðNð0Þ; Pð0ÞÞ, where

1069

6. Examples

1100

N

  LN aU P > r 1 ; K 1 þ b1 LN þ b2 U P

1068

1099

Now let us ask for LN such that

1048 1050

1054

ð67Þ

Remark 5. When parameters and initial conditions in model (5) are chosen according to the procedure stated at the end of the proof of Theorem 4, then the pest and predator will co-exist for all sufficiently large time t. In fact, by part (ii) of Theorem 4, the pest population will satisfy NðtÞ P LN > 0 for t P 0. In addition, from the proof of Theorem 2, we have PðtÞ P P 1 ðtÞ for t P 0 where e e P1 ðtÞ is a function such that P 1 ðtÞ ! PðtÞ as t ! 1 (where PðtÞ is defined in Lemma 1), and from the proof of Theorem 45, we have e PðtÞ P W for t P 0, where W ¼ P  Q is a positive constant (where P  is defined in Lemma 1). Hence, PðtÞ P W=2 > 0 for all t sufficiently large.

1130 1131 1132 1133 1134 1135 1136 1137 1138 1139

1144

1146 1147 1148 1149 1150 1151 1152 1153 1154 1155 1156 1157

MBS 7526

No. of Pages 14, Model 5G

5 September 2014 Q1

Table 1 Five examples for model (5). Where a parameter appears in the second column of the table, it is assumed to be positive. Assumptions in model (5)

Dynamics between release events    aNP  dN ¼ rN 1  NK  1þb , dt 1N h  i vaN dP ¼ P 1þb1 N  l dt

b2 ¼ 0, BðFðN; PÞ; PÞ ¼ vFðN; PÞ, DðFðN; PÞÞ ¼ l; m ¼ 0

(Rosenzweig–MacArthur model [40])     dN ¼ rN 1  NK  1þbaNP , dt 1 Nþb2 P h  i vaN dP ¼ P 1þb1 Nþb2 P  l dt

b2 > 0, BðFðN; PÞ; PÞ ¼ vFðN; PÞ, DðFðN; PÞÞ ¼ l; m ¼ 0

(DeAngelis model [40])    aNP  dN ¼ rN 1  NK  1þb , dt 1N h i  vaN dP ¼ P 1þb1 N  l  mP dt

b2 ¼ 0, BðFðN; PÞ; PÞ ¼ vFðN; PÞ, DðFðN; PÞÞ ¼ l; m > 0

(Bazykin model [40])     dN ¼ rN 1  NK  1þbaNP , dt 1 Nþb2 P h i  vaN dP ¼ P 1þb1 Nþb2 P  l  mP dt

b2 > 0, BðFðN; PÞ; PÞ ¼ vFðN; PÞ, DðFðN; PÞÞ ¼ l; m > 0

b2 ¼ 0, BðFðN; PÞ; PÞ ¼ vFðN; PÞ



P dþP

DðFðN; PÞÞ ¼ l; m ¼ 0



,

(original Beddington– DeAngelis model [13])    aNP  dN ¼ rN 1  NK  1þb , dt 1N h   i vaN dP P ¼ P 1þb1 N dþP  l dt (proposed by Verdy; see the model defined by Eqs. (33) and (36) in [41])

1

0

1

2

3

4

5

6

7

F (N,P) Fig. 1. Plots of the birth and death functions that are defined, respectively, in (75) and (78). The parameter choices used for these plots were: /1 ¼ 1:5; /2 ¼ 5; b ¼ 2:5; h1 ¼ 1:5; h2 ¼ 5; l ¼ 1:5; m ¼ 0:5.

1160 1161 1162 1163 1164 1165 1166 1167

1168

A per capita birth function for predators, with the same qualitative features as in (75), has been defined in Eq. (35) in [39]. An example of (75) is plotted in Fig. 1. As noted above, the per capita birth rate for predators in (74) is a scaled version of the birth rate in (73), where the scaling term, P=ðd þ PÞ, accounts for the possibility that sexually-reproducing predators may struggle to find a mate when the predator population P is small. Similarly, the birth rate in (75) may be adapted to account for such a possibility by scaling it by the term P=ðd þ PÞ, thereby yielding: 8 0; 0 6 FðN;PÞ 6 /1 ; > > n  h io  > < FðN;PÞ/1 2 p P ; /1 6 FðN;PÞ 6 /2 ; 1 þ /2 /1 dþP 2 BðFðN;PÞ;PÞ ¼ bcos >   > > P :b ; FðN;PÞ P /2 ; dþP

1170

ð76Þ

1171

The birth rate in (76) is consistent with the assumptions in model (5). We have not seen such a birth rate appear in any previous modelling studies. All of the examples for model (5) in Table 1 are such that:

1172 1173 1174

In words, Eq. (77) says that the per capita death rate for the predator, in the absence of intraspecific competition that directly causes mortality, is constant. Thus, when (77) holds, the per capita death rate for the predator is not directly dependent on the predator functional response or prey capture rate. However, it would seem reasonable to assume that predator mortality does depend directly on the prey capture rate, since prey consumption fuels all of the life processes of the predator, including, most fundamentally, staying alive. Accordingly, we define a form for the predator death function DðFðN; PÞÞ that is arguably more realistic than (77), based on the following assumptions: first, for all sufficiently small values of the prey capture rate FðN; PÞ, the death function DðFðN; PÞÞ is constant and at its maximum value, reflecting the idea that, for all sufficiently small prey capture rates, a predator suffers its highest risk of mortality from starvation or weakness brought on by lack of food; second, the death function DðFðN; PÞÞ falls across a positive range of values for FðN; PÞ, since a predator will be less likely to die, from starvation or weakness brought on by lack of food, as its prey capture rate increases above a certain value; third, the death function DðFðN; PÞÞ is constant and at its least value for all FðN; PÞ sufficiently large, reflecting the idea that increases above a certain level in the prey capture rate (and therefore, in the prey consumption rate) of an individual predator may have little impact on its short-term chance of death. Based on these assumptions, we define DðFðN; PÞÞ as follows:

1178 1179 1180 1181 1182 1183 1184 1185 1186 1187 1188 1189 1190 1191 1192 1193 1194 1195 1196 1197 1198 1199 1200 1201 1202

1203

1205

where h1 ; h2 ; l, and m are all positive constants, with h1 < h2 and m < l. Note that (78) is consistent with the assumptions in model

1206

(5). A death function for predators, with the same qualitative features as in (78), has been defined in Eq. (36) in [39]. An example of (78) is plotted in Fig. 1. Various examples for model (5), including (but not limited to) those stated in Table 1, can be constructed by combining one of the birth functions in (73)–(75), or (76), with one of the death functions in (77) or (78), and by assuming that the parameter b2 is either zero or positive, and that the parameter m is either zero or positive. In the next section, we shall explore, by simulation, the example defined as follows:

1208

model ð5Þ subject to ð76Þ and ð78Þ; where the parameters b2 and m are assumed to be positive:

 ð79Þ

1207 1209 1210 1211 1212 1213 1214 1215 1216 1217

1218 1220

7. Simulations

1221

In this section, we illustrate, by numerical simulation results, parts of our analysis from Sections 4 and 5, for the example of model (5) that is defined in (79). Our simulation results were created using MATLAB [1]. Parameter values were chosen with the single goal of corroborating our analytical results; we have not used real-world data but it would certainly be possible to seek to use such data in future work.

1222

7.1. Successful biocontrol

1229

We numerically explore our conditions for pest eradication from Section 4. To be specific, we explore the condition for local stability of the pest-eradication solution in (17), and the conditions for global stability of the pest-eradication solution in (18) and (46). From Corollaries 1 and 2, we know that, when a > b2 r and b2 > 0, conditions

1230

Q1 Please cite this article in press as: A.J. Terry, Biocontrol in an impulsive predator–prey model, Math. Biosci. (2014), http://dx.doi.org/10.1016/ j.mbs.2014.08.009

1175 1177

ð78Þ

B(F(N,P)) D(F(N,P))

2

1159

ð77Þ

l; 0 6 FðN; PÞ 6 h1 ; n  o FðN;PÞh1 2 p DðFðN; PÞÞ ¼ m þ ðl  mÞ cos 2 h2 h1 ; h1 6 FðN; PÞ 6 h2 ; > > : m; FðN; PÞ P h2 ;

3

0

DðFðN; PÞÞ ¼ l:

8 > > <

4

1158

11

A.J. Terry / Mathematical Biosciences xxx (2014) xxx–xxx

1223 1224 1225 1226 1227 1228

1231 1232 1233 1234

MBS 7526

No. of Pages 14, Model 5G

5 September 2014 Q1 1235 1236 1237 1238 1239 1240 1241 1242 1243 1244 1245 1246 1247 1248 1249 1250 1251 1252 1253 1254

12

A.J. Terry / Mathematical Biosciences xxx (2014) xxx–xxx

(17), (18), and (46) hold for sufficiently small inter-release periods T, and from the discussion in Section 4 after Corollary 2, we know that these three conditions do not hold for sufficiently large inter-release periods T. Therefore, assuming a > b2 r and b2 > 0, we may, for each of the three conditions in (17), (18), and (46), begin with a sufficiently small value of T, so that the condition in question will necessarily hold, and we may increase T in small increments until we reach a first value where it does not hold, and then we may choose T to be the value incrementally before this first value. By so doing, we may, for each condition, acquire a value for T which satisfies the condition and is maximal in some sense, and which we may thus call T max . Seeking such a T max for each condition is sensible – in the discussion in Section 4 after Corollary 2, we mentioned that it is practical and economical to seek to minimise the frequency of predator release events (which is equivalent to maximising T), whilst still ensuring control of the pest population. In Fig. 2(a), for the model in (79), we plot T max against the release yield Q, for each of the conditions in (17), (18), and (46), where T max is determined by the method described in the previous paragraph, and where model parameters are chosen as follows:

the local stability condition is open to debate, because we have seen in Theorem 4 that the pest can persist even when the pesteradication solution is locally stable. Interestingly, in Figs. 2(a)–(c), as Q grows, T max appears to be approaching a maximum value for the local and global stability conditions. If so, then this suggests that, no matter how many predators Q are dispersed on each release event, the events must always occur with some minimum frequency in order for the pest population to be eradicated locally or globally. Perhaps letting Q tend to infinity in the local and global stability conditions would yield some analytical insight into this issue – we leave this as an option for future work.

1280

7.2. Local but not global stability

1292

We show, by numerical simulation results, that there can be circumstances in which delaying a response to a new pest invasion can change the outcome of the response from pest eradication to pest persistence (see Remark 4). To do this, we first envisage the following sequence of events:

1293

l ¼ 0:8; m ¼ 0:05;

(E1) A new pest invasion occurs in a region at time t0 < 0. Let N ¼ NðtÞ be the pest population at time t, where t P t 0 . Since the pest invades the region at time t 0 , we assume Nðt0 Þ > 0. (E2) When the pest invades the region, it has no predators there. The pest population grows logistically until time t ¼ 0, that is, dN=dt ¼ rNð1  ðN=KÞÞ for t 2 ðt0 ; 0Þ. (E3) A predator is introduced to control the pest, where this occurs, from time t P 0 onwards, through a biocontrol program as in the model in (79), with the same values for the intrinsic growth rate r and environmental carrying capacity K of the pest as in (E2). We assume Nð0Þ ¼ Nð0 Þ where 0 is the time instantaneously before the time t ¼ 0. From (79), the first predator release event occurs at time t ¼ 0, so we assume Pð0Þ ¼ Q since predators are absent before the first release event.

1298

1281 1282 1283 1284 1285 1286 1287 1288 1289 1290 1291

1294 1295 1296 1297

1255

K ¼ 5;

b1 ¼ 0:5;

b2 ¼ 1;

a ¼ b2 r þ 0:1;

1257

ð80Þ

1258

and the value of r is shown under the figure. In exactly the same way, except for slightly changing the value of r, we created Figs. 2(b) and (c), where we state under each of these figures the value of r used in it. To evaluate the integrals in (17) and (18), we used the inbuilt function in MATLAB called ‘‘quad’’, which uses recursive adaptive Simpson quadrature [1]. In Fig. 2(a), we see that, for each particular Q > 0, the global stability condition in (18) (shown as ‘‘global stability 1’’) gives a larger value for T max than the global stability condition in (46) (shown as ‘‘global stability 2’’). However, in Fig. 2(c), the opposite is true. Moreover, in Fig. 2(b), we see that the global stability condition that gives a larger value for T max depends on Q, with (46) giving the larger value for T max for Q small, and (18) giving the larger value for T max for Q large. We may deduce from these observations that neither of the global stability conditions is always necessarily more practical or useful than the other. It would be of interest to compare the two conditions rigorously, but we leave this as a suggestion for future work. From Figs. 2(a)–(c), we see that, for each particular Q > 0, the local stability condition gives a larger value for T max than either of the global stability conditions. This is not surprising, since global stability implies local stability. Unfortunately, the practical value of

1259 1260 1261 1262 1263 1264 1265 1266 1267 1268 1269 1270 1271 1272 1273 1274 1275 1276 1277 1278 1279

12

6

local stability global stability 1 global stability 2

8 T

max

4

0

6

local stability global stability 1 global stability 2

100

Q

200

(a) r = 0.1

300

0

local stability global stability 1 global stability 2

4 T

max

2

0

max

2

0

100

Q

200

(b) r = 0.25

300

0

0

100

200

300

Q

(c) r = 0.4

Fig. 2. Plots of the maximal inter-release period T max versus the release yield Q, where T max is calculated, by the method in the first two paragraphs in Section 7.1, to satisfy the condition for the local stability of the pest-eradication solution in (17) (shown as ‘‘local stability’’) and the conditions for the global stability of the pest-eradication solution in (18) and (46), shown respectively as ‘‘global stability 1’’ and ‘‘global stability 2’’.

Q1 Please cite this article in press as: A.J. Terry, Biocontrol in an impulsive predator–prey model, Math. Biosci. (2014), http://dx.doi.org/10.1016/ j.mbs.2014.08.009

1300 1301 1302 1303 1304 1305 1306 1307 1308 1309 1310 1311 1312 1313 1314

From (E1) to (E3), we note that t 0 is the amount of time that passes between the new pest invasion and the start of the biocontrol program. Hence t0 measures the response time to the new pest invasion. Now let us choose values for the model parameters, the quantity LN , and the initial condition Pð0Þ in accordance with the sequence of events in (E1) to (E3), and in accordance with the procedure stated at the end of the proof of Theorem 4. Note that the

4 T

1299

1315 1316 1317 1318 1319 1320 1321 1322

MBS 7526

No. of Pages 14, Model 5G

5 September 2014 Q1

13

A.J. Terry / Mathematical Biosciences xxx (2014) xxx–xxx

N(t) P(t)

40

30

300 250 N(t) P(t)

200 150

20

100

10

50

0 −t 0

0

10

t

20

30

(a) t0 = 7.46, N (0) = 1.3633 < LN

0

−t0

−10

0

10

t

20

30

40

(b) t0 = 23, N (0) = 154.0301 > LN

Fig. 3. The two plots shown here were created in the same way (described in Section 7.2), with the only difference being the response time t0 from the new pest invasion (at time t ¼ t 0 ) to the first predator release event (at time t ¼ 0). Under each plot, we state the value chosen for t0 . Also under each plot, we state the value of the pest population at time t ¼ 0, namely Nð0Þ. If the response time is sufficiently small (plot (a)), then the pest is eradicated. If the response time is larger (plot (b)), then the pest persists, growing to a level very near its carrying capacity, namely K ¼ 300.

1323 1324 1325

1326

parameters h1 ; h2 , and d are not restricted by these criteria, but they are necessarily positive and we must also have h1 < h2 . Our choices are as follows:

b2 ¼ 1; r ¼ 0:35; a ¼ 10:35;

9 > > > =

/1 ¼ 1; /2 ¼ 2; b1 ¼ 200:

;

h1 ¼ 1; h2 ¼ 2; d ¼ 5;

l ¼ 1; m ¼ 0:05; Q ¼ 30; T ¼ 2; K ¼ 300; LN ¼ 150; m ¼ 0:5; Pð0Þ ¼ Q ¼ 30; v ¼ 0:8; b ¼ 0:8; > > >

1328

ð81Þ

1329

Next, observe that, for a new pest invasion, the initial pest population may be very small relative to the environmental carrying capacity K of the pest in the region that has been invaded. Thus, since we have chosen K ¼ 300 in (81), a reasonable choice for Nðt0 Þ (and one which we make) is Nðt0 Þ ¼ 0:1. From the sequence of events in (E1) to (E3), we know that the pest population NðtÞ will grow logistically from Nðt 0 Þ at time t ¼ t0 towards K until time t ¼ 0. In particular, since we have Nðt0 Þ ¼ 0:1 < LN ¼ 150 < 300 ¼ K, then NðtÞ will grow until it reaches LN and will then continue growing towards K, assuming the response time t0 is large enough. Numerical simulation reveals that Nð0Þ ¼ 154:0301 > LN when t 0 ¼ 23. Suppose we choose t 0 ¼ 23. Then we will have Nð0Þ 2 ½LN ; K. Then, in view of our method for choosing model parameters, the quantity LN , and initial condition Pð0Þ in (81), we see that all of the conditions in the procedure at the end of the proof of Theorem 4 have been satisfied for the model in (79). But then Theorem 4 holds, so that, by part (ii) of Theorem 4, we have NðtÞ P LN for t P 0. These observations are confirmed by simulation in Fig. 3(b); in performing the simulation, we used the MATLAB solver for systems of ordinary differential equations called ‘‘ode23t’’ before the first release event, as well as between consecutive release events. In fact, Fig. 3(b) shows not only that NðtÞ P LN for t P 0, but that the biocontrol program fails to have an impact on the pest population of any practical value. In the absence of the program, the pest population would tend to its carrying capacity K ¼ 300; by zooming in on part of Fig. 3(b) (results not shown), the pest population tends to a solution with small amplitude oscillations marginally below the carrying capacity. Fig. 3(a) was created in the same way as Fig. 3(b), except that the response time t0 was set to t 0 ¼ 7:46, giving Nð0Þ ¼ 1:3633 < LN . Theorem 4 is no longer guaranteed to hold for the model in (79), and indeed it does not. However, model parameters have been chosen according to the procedure at the end of the proof of Theorem 4, so that the pest-eradication solution is locally stable. Given this local stability, and the fact that the biocontrol

1330 1331 1332 1333 1334 1335 1336 1337 1338 1339 1340 1341 1342 1343 1344 1345 1346 1347 1348 1349 1350 1351 1352 1353 1354 1355 1356 1357 1358 1359 1360 1361 1362 1363 1364

program in Fig. 3(a) begins at a time when the prey population is still small compared to its carrying capacity K, the biocontrol program causes the prey population to tend to zero. We have found (results not shown) that if the response time is fractionally larger than in Fig. 3(a) (specifically, when t 0 ¼ 7:47), then the long-term model dynamics are the same as in Fig. 3(b). Thus, a very small additional delay in the response time can yield a profound difference in the long-term dynamics of the pest population. There is an obvious practical conclusion: in some circumstances, if we wish to control a new pest invasion, we should take action as quickly as possible. Such a conclusion highlights the importance of careful monitoring for new pest invasions, and of being prepared for swift action in case such an invasion is detected.

1365

8. Discussion

1378

We have studied a model for biological pest control (or ‘‘biocontrol’’) in which a pest population is controlled by a program of periodic releases of a fixed yield of predators that prey on the pest. Between releases, predator-pest dynamics were assumed to satisfy a predator–prey model with some fairly general properties. For example, the predator functional response was either of Beddington–DeAngelis type or Holling type II. We proved that, when the predator functional response was of Beddington–DeAngelis type and the predators were not sufficiently voracious, then the biocontrol program would fail to reduce the pest population below a particular economic threshold, regardless of the frequency or yield of the releases. We proved also that our model possessed a pest-eradication solution, which was both locally and globally stable provided that predators were sufficiently voracious and that releases occurred sufficiently often. We established, curiously, that the pest-eradication solution could be locally stable whilst not being globally stable, the upshot of which was that, if we were to delay a biocontrol response to a new pest invasion, then this could change the outcome of the response from pest eradication to pest persistence. We stated a number of specific examples for our model, and, for one of these examples, we corroborated parts of our analysis by numerical simulations. We list a few options for extending this work. Our biocontrol results have focused on the stability of the pest-eradication solution, but in many practical applications, it is enough to maintain a pest population below some given economic threshold rather than requiring its eradication. We could explore such applications by modifying our model to allow the times of the predator release events to be state-dependent, occurring whenever the pest popula-

1379

Q1 Please cite this article in press as: A.J. Terry, Biocontrol in an impulsive predator–prey model, Math. Biosci. (2014), http://dx.doi.org/10.1016/ j.mbs.2014.08.009

1366 1367 1368 1369 1370 1371 1372 1373 1374 1375 1376 1377

1380 1381 1382 1383 1384 1385 1386 1387 1388 1389 1390 1391 1392 1393 1394 1395 1396 1397 1398 1399 1400 1401 1402 1403 1404 1405 1406 1407

MBS 7526

No. of Pages 14, Model 5G

5 September 2014 Q1

14

A.J. Terry / Mathematical Biosciences xxx (2014) xxx–xxx

1421

tion reaches the given economic threshold [19]. Chaotic behaviour has been observed in numerous impulsively perturbed dynamical systems [29,10,42]; it would be of theoretical interest to explore the possibility of such behaviour in our model. An explicit spatial component could be added to our model by extending it to a patchy environment – some information on patch models may be found in [37,35,36,34]. Alternatively, a spatial component could be added by including diffusion terms for the pest and predator – note that a reaction-diffusion–advection model of biocontrol has previously been considered in [25]. Finally, our model could be modified to account for seasonal fluctuations in the pest population, by replacing, with a time-periodic function, the term that represents logistic growth for the pest population in the absence of predation.

1422

References

1423 1424 1425 1426 1427 1428 1429 1430 1431 1432 1433 1434 1435 1436 1437 1438 1439 1440 1441 1442 1443 1444 1445 1446 1447 1448 1449 1450 1451 1452 1453 1454 1455 1456 1457 1458 1459 1460 1461 1462

[1] http://www.mathworks.co.uk/products/matlab/. [2] P.A. Abrams, L.R. Ginzburg, The nature of predation: prey dependent, ratio dependent or neither? Trends Ecol. Evol. 15 (2000) 337–341. [3] D.D. Bainov, P.S. Simeonov, Impulsive Differential Equations: Periodic Solutions and Applications, Longman, London, 1993. [4] J.R. Beddington, Mutual interference between parasites or predators and its effect on searching efficiency, J. Animal Ecol. 44 (1975) 331–340. [5] J.R. Beddington, M.P. Hassell, J.H. Lawton, The components of arthropod predation: II. The predator rate of increase, J. Animal Ecol. 45 (1976) 165–185. [6] R. Boukadida, S. Michelakis, The use of Encarsiaformosa in integrated programs to control the whitefly Trialeurodesvaporariorum Westw. (Hom., Aleyrodidae) on greenhouse cucumber, J. Appl. Entomol. 118 (1994) 203–208. [7] L.E. Caltagirone, R.L. Doutt, The history of the vedalia beetle importation to California and its impact on the development of biological control, Annu. Rev. Entomol. 34 (1989) 1–16. [8] J.K. Conner, Artificial selection: a powerful tool for ecologists, Ecology 84 (2003) 1650–1660. [9] D.L. DeAngelis, R.A. Goldstein, R.V. O’Neill, A model for trophic interaction, Ecology 56 (1975) 881–892. [10] B. Deng, Equilibriumizing all food chain chaos through reproductive efficiency, Chaos 16 (043125) (2006). [11] W. Gao, S. Tang, The effects of impulsive releasing methods of natural enemies on pest control and dynamical complexity, Nonlinear Anal.: Hybrid Syst. 5 (2011) 540–553. [12] S.A. Gourley, Y. Kuang, A stage structured predator–prey model and its dependence on maturation delay and death rate, J. Math. Biol. 49 (2004) 188– 200. [13] M. Haque, A detailed study of the Beddington–DeAngelis predator–prey model, Math. Biosci. 234 (2011) 1–16. [14] G. Iooss, D. Joseph, Elementary Stability and Bifurcation Theory, Springer, New York, 1980. [15] P. Kindlmann, A.F.G. Dixon, Insect predator–prey dynamics and the biological control of aphids by ladybirds, in: R. van Driesche (Ed.), First international symposium on biological control of arthropods, USDA Forest Service, Morgantown, West Virginia, USA, 2003, pp. 118–124. [16] V. Lakshmikantham, D.D. Bainov, P.S. Simeonov, Theory of Impulsive Differential Equations, World Scientific, Singapore, 1989. [17] L. Mailleret, F. Grognard, Optimal release policy for prophylactic biological control, in: Positive Systems, Lecture Notes in Control and Information Sciences (LNCIS), vol. 341, Springer, Berlin, 2006, pp. 89–96.

1408 1409 1410 1411 1412 1413 1414 1415 1416 1417 1418 1419 1420

[18] L. Mailleret, F. Grognard, Global stability and optimisation of a general impulsive biological control model, Math. Biosci. 221 (2009) 91–100. [19] L. Nie, Z. Teng, L. Hu, J. Peng, Existence and stability of periodic solution of a predator–prey model with state-dependent impulsive effects, Math. Comput. Simul. 79 (2009) 2122–2134. [20] S. Nundloll, L. Mailleret, F. Grognard, Influence of intrapredatory interferences on impulsive biological control efficiency, Bull. Math. Biol. 72 (2010) 2113– 2138. [21] S. Nundloll, L. Mailleret, F. Grognard, Two models of interfering predators in impulsive biological control, J. Biol. Dynam. 4 (2010) 102–114. [22] D. Orr, Biological control and integrated pest management, in: R. Peshin, A.K. Dhawan (Eds.), Integrated Pest Management: Innovation-Development Process, Springer, Netherlands, 2009, pp. 207–239. [23] G. Pang, F. Wang, L. Chen, Extinction and permanence in delayed stagestructure predator–prey system with impulsive effects, Chaos Solit. Fract. 39 (2009) 2216–2224. [24] B.L. Phillips, G.P. Brown, J.K. Webb, R. Shine, Invasion and the evolution of speed in toads, Nature 439 (7078) (2006). 803–803. [25] N. Sapoukhina, Y. Tyutyunov, R. Arditi, The role of prey taxis in biological control: a spatial theoretical model, Am. Natural. 162 (2003) 61–76. [26] P.A. Shah, J.K. Pell, Entomopathogenic fungi as biological control agents, Appl. Microbiol. Biotechnol. 61 (2003) 413–423. [27] R. Shi, L. Chen, Staged-structured Lotka–Volterra predator–prey models for pest management, Appl. Math. Comput. 203 (2008) 258–265. [28] R. Shi, X. Jiang, L. Chen, A predator–prey model with disease in the prey and two impulses for integrated pest management, Appl. Math. Modell. 33 (2009) 2248–2256. [29] B. Shulgin, L. Stone, Z. Agur, Pulse vaccination strategy in the SIR epidemic model, Bull. Math. Biol. 60 (1998) 1123–1148. [30] R.R.L. Simons, S.A. Gourley, Extinction criteria in stage-structured population models with impulsive culling, SIAM J. Appl. Math. 66 (2006) 1853–1870. [31] G.T. Skalski, J.F. Gilliam, Functional responses with predator interference: viable alternatives to the Holling type II model, Ecology 82 (2001) 3083–3092. [32] V.M. Stern, R.F. Smith, R. van den Bosch, K.S. Hagen, The integrated control concept, Hilgardia 29 (1959) 81–101. [33] S. Tang, G. Tang, R.A. Cheke, Optimum timing for integrated pest management: modelling rates of pesticide application and natural enemy releases, J. Theor. Biol. 264 (2010) 623–638. [34] A.J. Terry, Control of Pests and Diseases, Ph.D. thesis, University of Surrey, 2009. [35] A.J. Terry, Impulsive culling of a structured population on two patches, J. Math. Biol. 61 (2010) 843–875. [36] A.J. Terry, Pulse vaccination strategies in a metapopulation SIR model, Math. Biosci. Eng. 7 (2010) 457–479. [37] A.J. Terry, Dynamics of a structured population on two patches, J. Math. Anal. Appl. 378 (2011) 1–15. [38] A.J. Terry, Prey resurgence from mortality events in predator–prey models, Nonlinear Anal.: Real World Appl. 14 (2013) 2180–2203. [39] A.J. Terry, A predator–prey model with generic birth and death rates for the predator, Math. Biosci. 248 (2014) 57–66. [40] P. Turchin, Complex Population Dynamics: A Theoretical/Empirical Synthesis, Princeton University Press, Princeton and Oxford, 2003. [41] A. Verdy, Modulation of predator–prey interactions by the Allee effect, Ecol. Modell. 221 (2010) 1098–1107. [42] L. Wang, L. Chen, J. Nieto, The dynamics of an epidemic model for pest control with impulsive effect, Nonlinear Anal.: Real World Appl. 11 (2010) 1374–1386. [43] H. Zhang, L. Chen, P. Georgescu, Impulsive control strategies for pest management, J. Biol. Syst. 15 (2007) 235–260. [44] H. Zhang, P. Georgescu, L. Chen, On the impulsive controllability and bifurcation of a predator–pest model of IPM, BioSystems 93 (2008) 151–171. [45] D. Zwillinger, Handbook of Differential Equations, Academic Press, New York, 1989.

1463 1464 1465 1466 1467 1468 1469 1470 1471 1472 1473 1474 1475 1476 1477 1478 1479 1480 1481 1482 1483 1484 1485 1486 1487 1488 1489 1490 1491 1492 1493 1494 1495 1496 1497 1498 1499 1500 1501 1502 1503 1504 1505 1506 1507 1508 1509 1510 1511 1512 1513 1514 1515 1516 1517 1518 1519 1520 1521 1522 1523 1524 1525

Q1 Please cite this article in press as: A.J. Terry, Biocontrol in an impulsive predator–prey model, Math. Biosci. (2014), http://dx.doi.org/10.1016/ j.mbs.2014.08.009