Dynamics for a non-autonomous predator-prey system with generalist predator

Journal Pre-proof Dynamics for a non-autonomous predator-prey system with generalist predator

Dingyong Bai, Jianshe Yu, Meng Fan, Yun Kang

PII:

S0022-247X(19)31088-1

DOI:

https://doi.org/10.1016/j.jmaa.2019.123820

Reference:

YJMAA 123820

To appear in:

Journal of Mathematical Analysis and Applications

Received date:

11 March 2019

Please cite this article as: D. Bai et al., Dynamics for a non-autonomous predator-prey system with generalist predator, J. Math. Anal. Appl. (2020), 123820, doi: https://doi.org/10.1016/j.jmaa.2019.123820.

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Published by Elsevier.

2

Dynamics for a non-autonomous predator-prey system with generalist predator

3

Dingyong Baia , Jianshe Yua , Meng Fanb , and Yun Kang∗c

1

a

School of Mathematics and Information Science, Guangzhou University, Guangzhou 510006, PR China; & Center for Applied Mathematics, Guangzhou University, Guangzhou 510006, PR China b School of Mathematics and Statistics, Northeast Normal University, Jilin 130024, PR China c Science and Mathematics Faculty, College of Integrative Sciences and Arts, Arizona State University, USA

4 5 6 7

8 9 10 11 12 13 14 15 16 17 18

Abstract: In this paper, we study complex dynamics of a non-autonomous predator-prey system with a Holling type II functional response and predator being generalist. We first studied basic dynamics such as boundedness, positive invariance, permanence, non-persistence and globally asymptotic stability. We also provide sufficient conditions for the existence, uniqueness and globally asymptotic stability of positive periodic solutions and boundary periodic solutions of the proposed model when parameters are periodic. We give the integral conditions to prove the extinction of prey and predator and globally asymptotic stability of boundary periodic solutions, and we show that these conditions have more reasonable biological interpretation than those expressed by supremum and infimum of parameters in some literature. We also perform numerical simulations to complement our analytical results and obtain more insights. Some of our numerical simulations indicate that periodic system may promote or suppress the permanence of its autonomous version with parameters being the averages of periodic parameters.

19

21

Keywords: Predator-prey system; persistence; non-permanence; globally asymptotic stability; periodic solutions

22

1

20

23 24 25 26 27 28 29 30 31 32

Introduction

Population dynamics of predator and prey have always been one of the central topics in mathematical ecology. Since the first prey-predator model was introduced by Volterra and Lotka [7, 8], many indepth studies on this topic have been done through developing different mathematical models of prey-predator interactions to obtain useful biological insights (e.g., [1, 2, 3, 4, 5, 6]). Many of these literature work have been focused on the constant environment, e.g., the life history parameters are constant. However, constant environment can’t reflect the reality of nature because seasonal effects of weather, food supply, mating habits, hunting or harvesting seasons, etc, would directly affect the important life history parameters being time-dependent. For this reason, time-dependent predator-prey systems have been widely investigated. A relatively early and representative work is by Cushing [9], in which the dynamics of a classical prey-predator periodic system was studied, and interesting results on the existence of periodic ∗ Corresponding author: Yun Kang. E-mail addresses of all authors: [email protected](Y.K), [email protected](J.Y), [email protected](D.B), [email protected](M.F).

1

1 2 3 4 5 6 7 8 9 10 11 12 13 14

solutions were provided. Since then, more attention has been gained on studying the dynamics of timedependent ecological models. For example, Samanta [10] studied the persistence and periodic solutions for a two-species nonautonomous Lotka-Volterra model in a polluted environment. Fan and Kuang [11] studied the permanence, extinction, global asymptotic stability for a nonautonomous predator-prey system with a Beddington-DeAngelis functional response, and obtained some results on the existence, uniqueness and stability of a positive (almost) periodic solution and a boundary (almost) periodic solution. Recently, Marv´ a et.al. [12] gave a simple geometrical condition for the existence of periodic solutions of planar periodic systems and applied it to a Holling’s type II predator-prey periodic model as well as a classic periodic competition model. We also refer to [13, 14] for the studies on the dynamics of ratiodependent non-autonomous systems and density-dependent non-autonomous systems. In [11, 13, 14], all parameters are continuous functions defined on R and bounded below by positive constants (except that the parameter α(t) of the Beddington-DeAngelis functional response is nonnegative in [11]), and the conditions for the permanence, extinction, global asymptotic stability of positive bounded solutions and periodic boundary solutions are expressed by supremum and infimum of parameters.

15 16 17

18

19 20 21 22 23 24 25 26 27

In this paper, we study the following nonautonomous prey-predator model with a Holling Type II functional response:     dx1 a(t)x2 (t) x1 (t) = x1 (t) r1 (t) 1 − − , (1.1) dt K1 (t) 1 + h(t)a(t)x1 (t)     dx2 e(t)a(t)x1 (t) x2 (t) = x2 (t) r2 (t) 1 − − d(t) + , (1.2) dt K2 (t) 1 + h(t)a(t)x1 (t) whose autonomous version is shown below:   dx1 = x1 r1 1 − dt   dx2 = x2 r2 1 − dt

  x1 ax2 − , K1 1 + hax1   x2 eax1 −d+ , K2 1 + hax1

(1.3) (1.4)

where x1 and x2 represent the population densities of prey and predator at time t, respectively; parameters ri , Ki (i = 1, 2), a, e, h, d have their specific ecological interpretations: r1 and r2 denote the intrinsic growth rate of species x1 and x2 , respectively, in the absence of other species; K1 and K2 are the carrying capacities of species x1 and x2 individuals, respectively; d stands for the dead rate of predator due to hunting or attacking all potential prey resources; a is the capturing efficiency of a predator and h is the predator handling time. In [15], Kang and Fewell studied a host-parasite co-evolutionary model, in which the relationship between host (prey) and parasite (predator) is described by (1.3)-(1.4). The ecological dynamics of (1.3)-(1.4), including the boundedness, permanence, stability of boundary and interior equilibria, and extinction of one species, were performed in [15].

28 29 30 31 32 33 34 35

In the system (1.3)-(1.4), both species x1 and x2 follow the Logistic growth in the absence of the other species. This indicates that predator species x2 is generalist and has alternative food resources in the absence of prey species x1 . In natural, many predators are generalist, and their preys consist of different species (see, for example, [16, 17, 18, 19, 20]). The study on prey-predator models involving generalist predators has attracted wide attentions, and provided additional biological insights on dynamical outcomes for models with predator being generalist versus specialist. For example, the recent work of Kang et.al. [15, 17] shows that models with generalist predator could exhibit much more complicate dynamics 2

1 2 3

and more likely to have "top down" regulation by comparing to the similar models with specialist predators. We also refer to [21, 22, 23, 24, 25] for more references on the study of prey-predator models with generalist predators.

4 5 6 7 8 9 10 11

In this paper, we study the dynamics of non-autonomous (1.1)-(1.2). We assume that all parameters r1 (t), r2 (t), K1 (t), K2 (t), d(t), a(t), h(t), e(t) are continuous functions defined on R, and bounded above by positive constants, in which r1 (t), r2 (t), K1 (t) and K2 (t) are bounded below by positive constants, and other parameter functions d(t), a(t), h(t) and e(t) are allowed to be nonnegative unless otherwise stated. We explore how time-dependent environment may contribute to the different dynamical outcomes of the corresponding system in the constant environment. We use differential inequalities, Brouwer fixed point theorem, coincidence degree theory and Lyapunov functions to prove our main results.

12 13 14

Comparing to the works in literature (see, for example [11, 13, 14]), our main results provide the following improvements: 1. The dynamics of a non-autonomous predator-prey system with Holling type II functional response and predator being generalist are considered. Some parameters, such as a(t), h(t), e(t) and d(t) are nonnegative function. The nonnegativity of those parameters are more reasonable than requiring them be bounded below by positive constants.

15 16 17 18

2. The conditions given to prove the extinction of prey and predator and globally asymptotic stability of boundary periodic solutions are presented by integral forms, which reflect the effects of the long term predation behaviors on the number of species, and has more reasonable biological interpretation than those for the corresponding autonomous system and those expressed by supremum and infimum of parameters for non-autonomous system in some literature (see, for example [11, 13, 14]).

19 20 21 22 23

3. Numerical simulations indicates that periodic system may promote or suppress the permanence of its autonomous version with parameters being the averages of periodic parameters.

24 25

26 27 28 29 30 31 32

The remaining sections of this paper is organized as follows. In Section 2, we provide basic dynamics of (1.1)-(1.2), such as boundedness, positive invariance, permanence, non-persistence and globally asymptotic stability. In Section 3, we explore sufficient conditions for the existence, uniqueness and globally asymptotic stability of positive periodic solutions and boundary periodic solutions of (1.1)-(1.2) with periodic parameters. In Section 4, we perform some interesting numerical simulations to complement our analytical results. Finally, we conclude our findings in Section 5 and provide the detailed proofs of our theoretical results in the last section.

33

34

35 36 37 38 39 40 41

2

Global dynamics of our proposed model

In this section, we first discuss the positive invariance and permanence of System (1.1)-(1.2), and then consider the global attractivity of the bounded positive solution of the proposed system. For the biological reasons, we only consider the solutions (x1 (t), x2 (t)) of (1.1)-(1.2) with nonnegative initial values, i.e., x1 (0) ≥ 0 and x2 (0) ≥ 0. For convenience of readers, we provide the following definitions. Definition 2.1. The solution set of system (1.1)-(1.2) is said to be ultimately bounded if there exist B > 0 such that for every solution (x1 (t), x2 (t)) of (1.1)-(1.2), there exists T > 0 such that ||(x(t), y(t))|| ≤ B for all t > T , where B is independent of particular solution while T may depend on the solution. 3

Definition 2.2. System (1.1)-(1.2) is said to be permanent if there exist positive constants δ, Δ with 0 < δ < Δ such that min lim inf xi (t) ≥ δ, max lim sup xi (t) ≤ Δ i∈{1,2} t→+∞

1 2

3 4 5

6 7

for all solutions of (1.1)-(1.2) with positive initial values. System (1.1)-(1.2) is said to be nonpersistent if there is a positive solution (x1 (t), x2 (t)) of (1.1)-(1.2) satisfying min lim inf xi (t) = 0. i∈{1,2} t→+∞

Definition 2.3. A bounded nonnegative solution (x∗1 (t), x∗2 (t)) of (1.1)-(1.2) is said to be globally asymptotically stable (or globally attractive) if any other solution (x1 (t), x2 (t)) of (1.1)-(1.2) with positive initial values satisfies lim (|x1 (t) − x∗1 (t)| + |x2 (t) − x∗2 (t)|) = 0. t→+∞

For a bounded continuous function f (t) on R, we use the notations f u := supt∈R f (t), f l := inf t∈R f (t). Define   e(t)a(t)K1u 2 (t) Mx2 = sup K r2 (t) r2 (t) − d(t) + 1+h(t)a(t)K1u , t∈R   (2.1) e(t)a(t)mx1 K2 (t) 1 (t) mx1 = inf K (r (t) − a(t)M ), m = inf (t) − d(t) + r , 1 x x 2 2 2 r1 (t) r2 (t) 1+h(t)a(t)mx t∈R

8

i∈{1,2} t→+∞

1

t∈R

and Γ := {(x1 , x2 ) ∈ R2 | mx1 ≤ x1 ≤ K1u , mx2 ≤ x2 ≤ Mx2 }. We define Condition (H1) as following ⎧ inf (r1 (t) − a(t)Mx2 ) > 0, ⎪ ⎪ ⎪ t∈R   ⎨ e(t)a(t)K1u > 0, inf r2 (t) − d(t) + 1+h(t)a(t)K u 1 t∈R   ⎪ ⎪ e(t)a(t)mx1 ⎪ ⎩ inf r (t) − d(t) + > 0. t∈R

9 10 11 12

2

(2.2)

(H1)

1+h(t)a(t)mx1

which implies that mx1 is smaller than K1u and hence mx2 is smaller than Mx2 . We would like to point out that the third inequality in (H1) can’t guarantee that Mx2 > 0 (i.e., the second inequality) holds. For example, we take constant coefficients K1 = 4, r1 = 0.04, K2 = r2 = 1, h = a = e = 1 and d = 1.9. 1 By a simple calculation, we get Mx2 = −0.1 < 0, mx1 = 14 > K1 and mx2 = 30 > 0.

13 14 15 16 17 18 19

Notes: Condition (H1) is important for system (1.1)-(1.2) being permanent. It suggests that (i) at any time t, the death rate d(t) of predator x2 due to hunting or attacking all potential prey resources is smaller than the sum of its intrinsic growth rate r2 (t) in the absence of prey x1 and the benefits by predating its prey x1 ; (ii) at any time t, the intrinsic growth rate r1 (t) of prey x1 in the absence of predator x2 is large such that inf (r1 (t) − a(t)Mx2 ) > 0 holds. t∈R

20 21 22

23 24

Now we provide the following lemma which is very important in the proof of the basic dynamical properties (e.g., permanence) of System (1.1)-(1.2). Lemma 2.1. Let p(t) and q(t) be continuous functions defined on R and bounded above and below by positive constants. If the positive function u(t) satisfies u (t) ≤ p(t)u(t) − q(t)u2 (t),

4

t ∈ [t0 , +∞),

1

p(t) p(t) then lim sup u(t) ≤ sup p(t) q(t) . Moreover, u(t) ≤ sup q(t) for all t ∈ [t0 , +∞) if 0 < u(t0 ) ≤ sup q(t) . While

2

if u(t) satisfies

t→+∞

t∈R

t∈R

t∈R

u (t) ≥ p(t)u(t) − q(t)u2 (t), 3

4 5 6 7

p(t) . t∈R q(t)

then lim inf u(t) ≥ inf t→+∞

p(t) t∈R q(t)

t ∈ [t0 , +∞), p(t) . t∈R q(t)

Moreover, u(t) ≥ inf

for all t ∈ [t0 , +∞) if u(t0 ) ≥ inf

Then we have the following theorem: Theorem 2.2. (The basic dynamic properties) Both the nonnegative and positive cones of R2 are positively invariant for (1.1)-(1.2). If Condition (H1) holds, then the set Γ is positively invariant for (1.1)-(1.2), and system (1.1)-(1.2) is permanent and the set Γ defined by Γ := {(x1 , x2 ) ∈ R2 | mx1 −  ≤ x1 ≤ Mx1 + , mx2 −  ≤ x2 ≤ Mx2 + }

8 9 10

(2.3)

is an ultimately bounded region of (1.1)-(1.2), where  > 0 sufficiently small such that mx1 −  > 0 and mx2 −  > 0. Notes: We would like to point out the following three implications from the theorem above:

11

1. If we take the methods in [11, 13, 14], then we have to assume that all parameters are bounded below by positive constants and satisfy the following condition (H1∗ ) ⎧ ˜ x > 0, ⎪ r l − au M 2 ⎪ ⎨ 1 el a l K u l u r2 − d + 1+hu au1K u > 0, (H1∗ ) 1 ⎪ l l ⎪ e a m ˜ x 1 ⎩ r l − du + > 0. u u 12 13 14 15

K2u r2l



u u

u



2

1+h a m ˜ x1

l

e a K K ˜ x ). By comparison, we can see that Condir2u − dl + 1+hl al K1 u , m ˜ x1 = ru1 (r1l − au M 2 1 1 tion (H1) is more flexible than (H1∗ ). In Theorem 2.2, not only parameter functions d(t), e(t), a(t) and h(t) are allowed to be nonnegative but also r2 (t) − d(t) could be negative at some time t or even in some period of time (see Example 4.1).

˜x = where M 2

16 17 18

2.

If the condition inf (r2 (t) − d(t)) > 0 holds, then the second and third inequalities of Condition t∈R

(H1) naturally hold, and it can be simplified as inf (r1 (t) − a(t)Mx2 ) > 0. t∈R

19 20 21 22 23 24 25 26 27 28 29

3. If all parameters in System (1.1)-(1.2) are positive constants, then Theorem 2.2 provides the following permanence results for the corresponding autonomous version of (1.3)-(1.4): eaK1 1 r2 (i) r2 > d and raK > r2 − d + 1+haK ; or 2 1 (ii) r2 < d and Mx∗2

r1 r2 aK2 m∗x1

> r2 − d +

eaK1 1+haK1

> r2 − d +

eam∗ x1 1+ham∗ x1

> 0.

where and are defined as (2.1) by all parameters being replaced by positive constants. Please note that these conditions are different to those given by Kang and Fewell (see Theorem 3.1, [15]). In [15], their permanent conditions for autonomous system (1.3)-(1.4) are given as eaK1 (i) r2 + 1+haK > d > r2 ; or 1 r1 K2 (ii) a > r2 (r2 − d) > 0, which were obtained by the results in Hutson ([26], Theorem 2.5 and its corollary).

30

31 32

The following theorem gives conditions under which System (1.1)-(1.2) is non-permanent (i.e., one species goes extinct in the system). 5

1

2 3 4 5

Theorem 2.3. (1)

If the following condition holds,  +∞  e(t)a(t)K1u r2 (t) − d(t) + dt = −∞ 1 + h(t)a(t)K1u 0

(2.4)

then the generalist predator x2 will be extinct; (2) Assume a(t) and h(t) are bounded below by positive constants, and inf (r2 (t) − d(t)) > 0, then prey t≥0

x1 will be extinct for system (1.1)-(1.2) if either 1 (t) (i) sup(h(t)a(t)K1 (t)) < 1 and m ˜ x2 > sup ra(t) , or t≥0

t≥0

2

1 (t)+1] inf (h(t)a(t)K1 (t)) > 1 and m ˜ x2 > sup r1 (t)[h(t)a(t)K . 4h(t)a2 (t)K1 (t)

6

(ii)

7

where m ˜ x2 = inf

8

Notes: If either sup(h(t)a(t)K1 (t)) > 1 or inf (h(t)a(t)K1 (t)) < 1, then the second part of Theorem 2.3

9

fails. However, we have the following extinct condition of integral form for prey x1 .

t≥0

K2 (t) (r2 (t) t≥0 r2 (t)

t≥0

− d(t)). t≥0

t≥0

10

11

Corollary 2.4. If inf (r2 (t) − d(t)) > 0 and the following equality holds t≥0



+∞



a(t)m ˜ x2 r1 (t) − 1 + h(t)a(t)K1u

0

 dt = −∞.

12

then prey x1 goes extinct for System (1.1)-(1.2).

13

Notes: We would like to point out the following two implications from the theorem above:

(2.5)

14 15 16 17 18 19

1. Based on Theorem 2.3, we can get the following extinct conditions of one species for the corresponding autonomous system (1.3)-(1.4). eaK1 (1) species x2 will be extinct if r2 + 1+haK < d, which is consistent with the extinct condition of species 1 x2 given by Kang and Fewell [15] (see Theorem 3.3 in [15]); r1 1) (2) species x1 will be extinct if (i) 1 − rd2 > r1 (1+haK , or (ii) haK1 < 1 and 1 − rd2 > aK , or (iii) aK2 2 2

20 21

(1+haK1 ) haK1 > 1 and 1 − rd2 > r14ha . The last two conditions are same, respectively, with the second and 2K K 1 2 third extinct condition of species x1 given by Kang and Fewell [15] (see Theorem 3.3 in [15]).

22 23 24 25 26 27 28 29 30

2. Condition (2.4) and (2.5) have more reasonable biological interpretation than those expressed by supremum and infimum of parameters in the case all parameters are bounded below by positive constants, and those for the corresponding autonomous system (1.3)-(1.4). Condition (2.4) indicates that for a long period of time the total sum of predation benefits and the growth of predator x2 is far less than the total mortality in predator x2 , the predator x2 will be extinct for system (1.1)-(1.2). Condition (2.5) shows that the long term effects of the predation behavior to prey x1 being larger than its growth rate leads e(t)a(t)K1u to the extinction of prey. Condition (2.4) and (2.5) also indicate that r2 (t) − d(t) + 1+h(t)a(t)K u and r1 (t) −

a(t)m ˜ x2 1+h(t)a(t)K1u

1

are allowed to change their signs (see Example 4.2).

31

32 33 34

Now we discuss the global stabilities of a positive solution of System (1.1)-(1.2). First, we introduce the following lemma from [27], which can be employed in establishing the globally asymptotic stability of System (1.1)-(1.2). 6

1 2

3 4

Lemma 2.5. Let t0 be a real number and f be a nonnegative function defined on [t0 , +∞) such that f is integrable on [t0 , +∞) and is uniformly continuous on [t0 , +∞). Then lim f (t) = 0. t→+∞

Theorem 2.6. Let (x∗1 (t), x∗2 (t)) be a bounded positive solution of System (1.1)-(1.2). Then (x∗1 (t), x∗2 (t)) is globally asymptotically stable if Condition (H1) holds and the following inequalities hold for all t ∈ R,   ⎧ h(t)a2 (t)Mx2 r1 (t) e(t)a(t) ⎨ inf K > 0, − (1+h(t)a(t)m 2 − (1+h(t)a(t)m 2 (t) ) ) 1 x x 1 1 t∈R  . (2.6) a(t) ⎩ inf r2 (t) − > 0, K2 (t) 1+h(t)a(t)mx 1

t∈R

5

6 7

8 9 10

11 12 13 14

3

Dynamics with periodic coefficients

In the section, we assume that all of parameters in System (1.1)-(1.2) are ω-periodic in t, and we study the existence of positive periodic solutions and the global attractivity of boundary periodic solutions. Theorem 3.1. If Condition (H1) holds, then System (1.1)-(1.2) has at least one positive ω-periodic solution (x∗1 (t), x∗2 (t)) ∈ Γ. If in addition, the assumption (2.6) in Theorem 2.6 holds, then System (1.1)-(1.2) has a unique positive ω-periodic solution that is globally asymptotically stable. Remark 3.1. If all the parameters in system (1.1)-(1.2) are positive constants, then the positive periodic solution obtained in Theorems 3.1 degenerates to a trivial positive periodic solution, i.e., the positive equilibrium E ∗ (x∗1 , x∗2 ) of the corresponding constant coefficients System (1.3)-(1.4). Thus, from Theorems 3.1 and Theorem 2.6, we can conclude that E ∗ is globally asymptotically stable for System (1.3)-(1.4) if



 ⎧ eamx1 r1 eaK1 1 1 ⎪ r > r > 0, > − d + − d + 2 2 ⎪ aK r 1+haK r 1+ham 2 2 1 2 x ⎨ 1 ⎪ ⎪ ⎩

15 16

r1 K1 r2 K2

− −

ha2 Mx2 ea (1+hamx1 )2 − (1+hamx1 )2 a 1+hamx1 > 0.

> 0,

r2 r1 r1 a 1 Since K r2 − d + − > 0 is equivalent to > + 1+ham aK r K1 r 2 x1 2 2 2 can be rewritten as follows ⎧ eamx1 > d, r2 + 1+ham ⎪ ⎪ x1 ⎨ haMx2 r1 e aK1 > (1+hamx1 )2 + (1+hamx1 )2 ,  ⎪ ⎪ ⎩ r2 > 1 + 1 r2 − d + eaK1 . aK2

K1

eaK1 1+haK1

 , thus, the condition above

1+haK1

r1

17

where Mx2 and mx1 are defined as (2.1) by all parameters being replaced by positive numbers.

18

Now, we use a continuation theorem in coincidence degree theory to establish basic criteria for the same problem but in terms of the averages of the related parameters over an interval of the common period. Let X, Y be normed vector spaces, L : DomL ⊂ X→Y be a linear mapping, N : X→Y be a continuous mapping. The mapping L will be called a Fredholm mapping of index zero if dimKerL = codim ImL < +∞ and ImL is closed in Y . If L is a Fredholm mapping of index zero and there exist continuous projections P : X→X and Q : Y →Y such that ImP = KerL, ImL = KerQ = Im(I − Q), it follows that L|DomL ∩ KerP : (I − P )X→ImL is invertible. We denote the inverse of that map by L−1 P . ¯ ¯ If Ω is an open bounded subset of X, the mapping N will be called L-compact on Ω if QN (Ω) is bounded ¯ and KP−1 (I − Q)N : Ω→X is compact. Since ImQ is isomorphic to KerL, there exists an isomorphism J : ImQ→KerL. The following lemma is from Gains and Mawhin [28].

19 20 21 22 23 24 25 26 27

7

1 2 3 4

¯ Suppose Lemma 3.2. Let L be a Fredholm mapping of index zero and N be L-compact on Ω. (1) for each λ ∈ (0, 1), every solution x of Lx = λN x is such that x ∈ / ∂Ω; (2) QN x = 0 for each x ∈ ∂Ω ∩ KerL and the Brouwer degree deg(JQN, Ω ∩ KerL, 0) = 0. ¯ Then the operator equation Lx = N x has at least one solution lying in DomL ∩ Ω. For a continuous and periodic function g(t) with period ω, we denote 1 ω g(t)dt. A(g) := ω 0

5

For convenience, let S1 = S3 =

A(r1 ) A(r1 /K1 ) ,

S2 =

A(r1 −aS2 exp{2ωA(r2 )}) , A(r1 /K1 )

S4 =

  eaS1 exp{2ωA(r1 )} A r2 −d+ 1+haS exp{2ωA(r )} 1

A



A(r2 /K2 )

1

eaS3 exp{−2ωA(r1 )} r2 −d+ 1+haS 3 exp{−2ωA(r1 )}

A(r2 /K2 )

, 

.

Theorem 3.3. If the conditions stated in (H2) hold,  ⎧ eaS1 exp{2ωA(r1 )} ⎪ A r − d + 2 ⎪ 1+haS1 exp{2ωA(r1 )} > 0, ⎪ ⎨ A (r1 − aS2 exp{2ωA(r2 )}) > 0, ⎪  ⎪ ⎪ eaS3 exp{−2ωA(r1 )} ⎩ A r2 − d + 1+haS > 0, 3 exp{−2ωA(r1 )} 6

(H2)

which are equivalent to A(r1 )A(r2 /K2 ) A(a)

 eaS1 exp{2ωA(r1 )} A(r2 − d) + A 1 + haS1 exp{2ωA(r1 )}   eaS3 exp{−2ωA(r1 )} > 0, A(r2 − d) + A 1 + haS3 exp{−2ωA(r1 )} 

> >

7

then System (1.1)-(1.2) has at least one positive ω-periodic solution.

8

Note: If all the parameters in System (1.1)-(1.2) are positive constants, then Condition (H1) is same as (H2) as ω→0.

9 10

11 12 13

Now we consider the existence and global stability of boundary ω-periodic solutions of (1.1)-(1.2). Since ri (t), Ki (t)(i = 1, 2) are positive ω-periodic functions, System (1.1)-(1.2) always has boundary ω-periodic solution (x∗1 (t), 0), where x∗1 (t)

14

15 16

 = exp



ω 0



 

t+ω

r1 (s)ds − 1 t

 t  −1 r1 (s) exp − r1 (τ )dτ ds , K1 (s) s

(3.1)

which is the unique positive ω-periodic solution of Logistic equation (see Theorem 2.1 in [14, 29])   r1 (t)x1 (t) . x1 (t) = x1 (t) r1 (t) − K1 (t) ω In addition, if 0 (r2 (s) − d(s))ds > 0, (1.1)-(1.2) has another boundary ω-periodic solutions (0, x∗2 (t)), where  t   ω    t+ω  −1 r2 (s) ∗ exp − x2 (t) = exp (r2 (s) − d(s))ds − 1 (r2 (τ ) − d(τ ))dτ ds , (3.2) K2 (s) 0 t s 8

1

which is the unique positive ω-periodic solutions of Logistic equation   r2 (t)x2 (t) . x2 (t) = x2 (t) (r2 (t) − d(t)) − K2 (t)

2

Theorem 3.4. The following two statements are true:

3 4 5

6 7 8 9

1. Boundary periodic solution (x∗1 (t), 0) is globally asymptotically stable for system (1.1)-(1.2) if the following inequality holds  ω e(t)a(t)K1u r2 (t) − d(t) + dt < 0. (3.3) 1 + h(t)a(t)K1u 0 ω 2. Assume a(t), h(t) are bounded below by positive constants, and 0 (r2 (s) − d(s))ds > 0. If one of the following conditions is satisfied, then (0, x∗2 (t)) is globally asymptotically stable for system (1.1)-(1.2). a(t)m ˜∗

e(t)a(t)K u

x2 1 (i) r1 (t) + 1+h(t)a(t)K t ∈ [0, ω]; u < 1+h(t)a(t)K u , 1 1 (ii) h(t)a(t)K1 (t) < 1, r1 (t) − a(t)m ˜ ∗x2 < 0, t ∈ [0, ω];

r1 (t)[h(t)a(t)K1 (t)+1]2 , 4h(t)a2 (t)K1 (t) (r2 (t)−d(t))K2 (t) inf . r2 (t) t∈[0,ω]

(iii)

h(t)a(t)K1 (t) > 1, m ˜ ∗x2 >

11

where

m ˜ ∗x2

12

If one condition of the three conditions (i), (ii) and (iii) in the second part of Theorem 3.4 isn’t satisfied, then the second part of Theorem 3.4 fails. However, similar to Corollary 2.4, we have the following condition in an integral form for the globally asymptotical stability of (0, x∗2 (t)). ω Corollary 3.5. Assume that 0 (r2 (s) − d(s))ds > 0. Boundary periodic solution (0, x∗2 (t)) is globally asymptotically stable if  ω a(t)m ˜ ∗x2 r1 (t) − dt < 0. (3.4) 1 + h(t)a(t)K1u 0

10

13 14

15 16

17

=

t ∈ [0, ω].

Notes: We would like to point out the following four implications from the theoretic results above:

18 19 20 21 22 23 24 25 26

1. If all the parameters in System (1.1)-(1.2) are positive constants, then the boundary periodic  ∗ ∗ solutions (x1 (t), 0) and (0, x2 (t)) degenerate to the boundary equilibria (K1 , 0) and 0, (1 − rd2 )K2 of the autonomous version of (1.3)-(1.4), respectively. Theorem 3.4 indicates that eaK1 (1) (K1 , 0) is globally asymptotically stable for system (1.3)-(1.4) if r2 + 1+haK < d; 1

 d (2) 0, (1 − r2 )K2 is globally asymptotically stable if either (i) r1 r2 (1 + haK1 ) + ear2 K1 < a(r2 − d)K2 , or r1 (ii) haK1 < 1 and 1 − rd2 > aK , or 2 (iii) haK1 > 1 and 1 −

d r2

>

r1 (1+haK1 )2 4ha2 K1 K2 .

27 28 29 30

2. Similar to (2.4) and (2.5), Condition (3.3) and (3.4) have more reasonable biological interpretations than those expressed by supremum and infimum of parameters in the case all parameters are bounded below by positive constants.

31 32 33 34

3. Either Condition (3.3) or (3.4) implies that System (1.1)-(1.2) is not permanent. However, if all parameters are replaced by their averages A(·) on the periodic interval [0, ω], the corresponding autonomous system (1.3)-(1.4) may be permanent (See Example 4.6).

35

9

a(t)m ˜∗

e(t)a(t)K u

2

x2 1 4. Conditions (3.3) and (3.4) imply that r2 (t)−d(t)+ 1+h(t)a(t)K u and r1 (t)− 1+h(t)a(t)K u are allowed 1 1 to change their signs on [0, ω](See Example 4.5 and 4.6).

3

4

1

4 5 6

Numerical simulations

In this section, we perform different numerical scenarios to verify our theoretical results. We also use numerical simulations to illustrate that the periodic system may promote or suppress the permanence of its autonomous version with parameters being the averages of the corresponding periodic parameters.

7 8

Recall that Γ is defined as follows: Γ := {(x1 , x2 ) ∈ R2 | mx1 −  ≤ x1 ≤ Mx1 + , mx2 −  ≤ x2 ≤ Mx2 + }.

9

10 11 12

(4.1)

where  > 0 is sufficiently small such that mx1 −  > 0 and mx2 −  > 0. Example 4.1. [Theorem 2.2,2.6,3.1] Take r1 (t) = 9.5 + 51 sin t, K1 (t) = 20 + 15 sin t, r2 = 4.9, K2 = 10 + 15 cos t, d(t) = 2.5(1 + cos t), a(t) = 0.3(1 + cos t), e(t) = 0.2(1 + cos t), h(t) = 0.1(1 + sin t), then System (1.1)-(1.2) becomes ⎧     x1 (t) (0.3(1+cos t))x2 (t) 1 1 ⎪ = x (t) 9.5 + sin t 1 − − ⎨ dx 1 1 dt 5 1+0.03(1+sin t)(1+cos t)x1 (t) 20+ 5 sin t (4.2)    ⎪ 0.06(1+cos t)2 x1 (t) ⎩ dx2 = x2 (t) 4.9 1 − x21(t) − 2.5(1 + cos t) + dt 1+0.03(1+sin t)(1+cos t)x1 (t) 10+ cos t 5

13 14 15

where parameters d(t), a(t), e(t) and h(t) all are nonnegative functions with zero lower bounds, and 24 d(t) ≥ r2 (t) for t ∈ [2kπ−arccos 24 25 , 2kπ+arccos 25 ], k ∈ Z, otherwise, d(t) < r2 (t). By a simple numerical computation, we have Mx1 = K1u = 20.2,

16

Mx2 = 9.8,

mx1 ≈ 7.6119,

mx2 ≈ 2.3113.

Hence, the parameters of System (4.2) satisfy Condition (H1).

17

1.

Due to Theorem 2.2, the set Γ = {(x1 , x2 ) ∈ R+ : 7.6119 ≤ x1 ≤ 20.2, 2.3113 ≤ x2 ≤ 9.8}

18 19

is a positive invariant set of system (4.2), which also can be seen from Figure 1. In Figure 1, the solutions starting from the points (18, 3) and (8, 9) ∈ Γ, respectively, stay in Γ for all t ≥ 0.

20

2.

Choose  = 0.01. Then Theorem 2.2 implies that the set Γ = {(x1 , x2 ) ∈ R+ : 7.6019 ≤ x1 ≤ 20.21, 2.3013 ≤ x2 ≤ 9.81}

21 22

is an ultimately bounded region and System (4.2) is permanent. See Figure 2, in which the solutions of System (4.2) with initial points (23, 1) and (6, 1), respectively, stay in Γ for large t > 0.

23 24 25

3. System (4.2) is a periodic system with ω = 2π. One can check that the condition (2.6) in Theorem 2.6 holds. In fact,   ⎧ h(t)a2 (t)Mx2 r1 (t) e(t)a(t) ⎨ inf K ≈ 0.1952 > 0, − − 2 2 (1+h(t)a(t)mx1 ) (1+h(t)a(t)mx1 ) t∈R  1 (t)  a(t) ⎩ inf r2 (t) − ≈ 0.0185 > 0. K2 (t) 1+h(t)a(t)mx t∈R

1

10

1 2 3

Therefore, according to Theorem 3.1, System (4.2) has a unique positive periodic solution which is globally asymptotically stable. See Figure 3 and Figure 4. In Figure 3, phase diagrams start from initial points (1, 3), (5, 15), (25, 15) and (30, 1), respectively. In Figure 4, the orbit starts from (15, 5).

4

Figure 1: Positive invariant set Γ of System (4.2).

Figure 2: Permanence and ultimately bounded region Γ ( = 0.01) of System (4.2). Here, the orbits start from initial points (23, 1) and (6, 1), respectively, stay in Γ for large t > 0.

Here, the orbits start from initial points (18, 3), (8, 9) ∈ Γ, respectively, stay in Γ for all t ≥ 0.

15

20

x1

18 16 10 14

0

20

40

60

80

100

60

80

100

x2

t 10 5 x2

8 6 4

0

0

5

10

15 x1

20

25

2

30

7

40

Figure 4: Unique periodic solution of System (4.2) with initial value (15, 5).

(4.2). Here, the solutions starting from initial points (1, 3), (5, 15), (25, 15) and (30, 1), respectively, tend to a stable periodic solution.

6

20

t

Figure 3: Globally asymptotical stability of System

5

0

t Example 4.2. [Corollary 2.4] Take r1 (t) = 0.18(1.02+sin t) exp( 1+t ), K1 = 15, r2 (t) = 3+ 15 cos t, K2 (t) = 1 15 exp( t+1 ) + cos t, d(t) = 0.2(1 + cos t), a(t) = 0.2(1 + cos t), e(t) = 0.2(1 + sin t), h(t) = 0.2(1 + cos t), then System (1.1)-(1.2) becomes ⎧     x1 (t) 0.2(1+cos t)x2 (t) t 1 ⎪ = x (t) 0.18(1.02 + sin t) exp( ) 1 − − ⎨ dx 1 dt 1+t 15 1+0.04(1+cos t)2 x1 (t) (4.3)   

 ⎪ x2 (t) 0.04(1+cos t)(1+sin t)x1 (t) ⎩ dx2 = x2 (t) 3 + 1 cos t 1 − − 0.2(1 + cos t) + 1 2 dt 5 1+0.04(1+cos t) x1 (t) 15 exp( )+cos t t+1

11

1 2 3

Here, the death rate d(t) of predator x2 is always smaller than its intrinsic growth rate r2 (t). Parameters d(t), a(t), h(t) and e(t) are nonnegative. Clearly, h(t)a(t)K1 (t) = 1 + 0.6(1 + cos t)2 and sup(h(t)a(t)K1 (t)) = 2.4 > 1 and inf (h(t)a(t)K1 (t)) = 0 < 1. Therefore, the extinct conditions of t≥0

t≥0

4 5

species x1 given in the second part of Theorem 2.3 don’t satisfied. But we can use Corollary 2.4 to show that prey x1 goes extinct for System (4.3). In fact, by a simple numerical computation, we have K2 (t) (r2 (t) − d(t)) = 14, t∈R r2 (t) a(t)m ˜ x2 := r1 (t) − 1+h(t)a(t)K u = 0.18(1.02 1

m ˜ x2 = inf f (t)

7

a(t)m ˜

f(t)

8

2.8(1+cos t) 1+0.6(1+cos t)2 .

x2 From the graph of f (t) = r1 (t) − 1+h(t)a(t)K u , see Figure 5, we can see that its sign constantly changes

1  +∞ a(t)m ˜ x2 with its integral 0 r1 (t) − 1+h(t)a(t)K u dt = −∞. Therefore, by Corollary 2.4, prey x1 will be 1 extinct for System (4.3). It also can be seen from Figure 6, in which the orbits start from initial points (3, 0.5), (2, 0.2) and (1, 1), respectively.

1

3

0.5

2.5

0

2

−0.5

x1

6

t + sin t) exp( 1+t )−

1.5

−1

1

−1.5

0.5

−2

0

20

40

60

80

0

100

t

0

5

10

15

20

25

30

t

Figure 5: The graph of f (t) = r1 (t) − for System (4.3).

a(t)m ˜ x2 1+h(t)a(t)K1u

Figure 6: Orbits of (4.3): Prey x1 will be extinct. The

orbits start from (3, 0.5), (2, 0.1) and (1, 1), respectively.

9

10 11

12 13

Example 4.3. [Theorem 2.3] Take r1 (t) = 1.02 + sin t, K1 = 10, r2 = 0.5, K2 = 15, d(t) = 2.5(1 + t+1) 0.9(sin t+1) cos t), a(t) = 50(sin , h(t) = 0.2(1 + sin t), then System (1.1)-(1.2) becomes (t+1)2 , e(t) = (t+1)2 ⎧ ⎪ ⎪ ⎪ ⎪ ⎨

dx1 dt

⎪ ⎪ ⎪ ⎪ ⎩

dx2 dt



15

1+

(t+1)2

(4.4)

x1 (t)

  One can see that r2 (t) ≥ d(t) for t ∈ (2k + 1)π − arccos 45 , (2k + 1)π + arccos 45 , k ∈ Z+ , otherwise, r2 (t) < d(t). Also we have the follows r2 (t) − d(t) +

14

 50(sin t+1) 

x2 (t) x1 (t) (t+1)2 = x1 (t) (1.02 + sin t) 1 − 10 − 10(sin t+1)2 1+ (t+1)2 x1 (t)   45(sin t+1)2

 x1 (t) x2 (t) (t+1)4 = x2 (t) 0.5 1 − 15 − 2.5(1 + cos t) + 10(sin t+1)2

e(t)a(t)K1u 450(sin t + 1)2 = −2 − 2.5 cos t + u 4 1 + h(t)a(t)K1 (t + 1) + 100(t + 1)2 (sin t + 1)2

 +∞  +∞ being constantly changes its sign in [0, +∞). We can conclude that 0 (r2 (t) − d(t))dt = 0 (−2 −  +∞  +∞ 45(1+sin t)2 2.5 cos t)dt = −∞. Also, we can check that 0 e(t)a(t)dt = 0 dt < 220 (t+1)4 3 . Therefore, (2.4) 12

1 2 3 4 5

6 7 8 9 10

holds and by the first part of Theorem 2.3, predator x2 will be extinct for System (4.4). This also can be seen from Figure 7, in which the orbits start from initial points (16, 1.5), (30, 0.5) and (20, 1), respectively. Example 4.4. [Theorem 3.3] Let r1 (t) = 2 + cos(4πt), r2 (t) = 0.198, K1 (t) = 12 + 6 cos(4πt), K2 (t) = + 34 cos(4πt), a(t) = 0.1(1 + cos(4πt)), e(t) = 0.05(1 + cos(4πt)), h(t) = 0.1(1 + sin(4πt)), d(t) = 0.1(1 + cos(4πt)). Then System (1.1)-(1.2) becomes the following periodic system with ω = 12 : ⎧   

x1 (t) 0.1(1+cos(4πt))x2 (t) 1 ⎪ − = x (t) (2 + cos(4πt)) 1 − ⎨ dx 1 dt 12+6 cos(4πt) 1+0.01(1+cos(4πt))(1+sin(4πt))x1 (t) 

  (4.5) ⎪ 0.005(1+cos(4πt))2 x1 (t) ⎩ dx2 = x2 (t) 0.198 1 − 3 3x2 (t) − 0.1(1 + cos(4πt)) + dt 1+0.01(1+cos(4πt))(1+sin(4πt))x1 (t) 2 + 4 cos(4πt)   1 1 where al = 0, el = 0, dl = 0 and d(t) ≤ r2 (t) for t ∈ k2 − 4π arccos 0.98, k+1 2 + 4π arccos 0.98 , k ∈ Z, otherwise, d(t) > r2 (t). By a simple numerical computation, we have A(r2 − d) = 0.09 > 0 and A(r1 − aS2 exp 2ωA(r2 )) ≥ 1.4 > 0. Then, by Theorem 3.3, System (4.5) has at least one positive 12 periodic solution, which also can be seen from Figure 8. In Figure 8, the trajectory starts from initial point (6, 1). 3 2

1.5

x1

15

10 1 5

0

20

40

60

80

100

120

140

160

180

200

120

140

160

180

200

x2

t 1.4 0.5

x2

1.2 1

0

0.8 0

5

10

15

20

25

30

35

40

45

50

0

20

40

60

Figure 7: Orbits of System (4.4): predator x2 goes extinct. The orbits start from (16, 1.5), (30, 0.5) and (20, 1), respectively.

11 12

13 14

15

80

100

t

t

Figure 8: Positive periodic solution of System (4.5) with initial value (6, 1).

Example 4.5. [Corollary 3.5] In Example 4.2, we take r1 (t) = 0.5(1.02 + sin t), K2 (t) = 15 + cos t, then we have the following 2π-periodic System (1.1)-(1.2) ⎧ 

  x1 (t) 0.2(1+cos t)x2 (t) 1 ⎪ = x (t) 0.5(1.02 + sin t) 1 − − ⎨ dx 2 1 dt 15 1+0.04(1+cos t) x1 (t) (4.6)   

 ⎪ 0.04(1+cos t)(1+sin t)x1 (t) ⎩ dx2 = x2 (t) 3 + 1 cos t 1 − x2 (t) − 0.2(1 + cos t) + 2 dt 5 15+cos t 1+0.04(1+cos t) x1 (t) By (3.2), we have the periodic boundary solution (0, x2 ∗ (t)) = (0, 14) based on the same arguments in Example 4.2, sup(h(t)a(t)K1 (t)) = 2.4 > 1 and inf (h(t)a(t)K1 (t)) = 0 < 1. We can also verify that r1 (t) +

t≥0 a(t)m ˜∗ e(t)a(t)K1u x2 1+h(t)a(t)K1u ≮ 1+h(t)a(t)K1u , t

t≥0

∈ [0, 2π]. In fact, we have  a(t)m ˜ ∗x2 e(t)a(t)K1u |t=0 , r1 (t) + |t=0 = 0.8629 < 1.6471 = u 1 + h(t)a(t)K1 1 + h(t)a(t)K1u   a(t)m ˜ ∗x2 e(t)a(t)K1u |t=π = 0.51 > 0 = |t=π . r1 (t) + u 1 + h(t)a(t)K1 1 + h(t)a(t)K1u



13

1 2 3

Therefore, the globally stable conditions of boundary periodic solution (0, x2 ∗ (t)) given in the second part of Theorem 3.4 aren’t satisfied for System (4.6). Note that K2 (t) r2 (t) (r2 (t) − d(t)) = 14, a(t)m ˜∗ x2 r1 (t) − 1+h(t)a(t)K u = 0.5(1.02 1

m ˜ ∗x2 = inf+ t∈R

g(t) :=

5 6

a(t)m ˜∗

g(t)

7

2.8(1+cos t) 1+0.6(1+cos t)2 .

x2 From the graph of g(t) = r1 (t)− 1+h(t)a(t)K u (see Figure 9), we can see that its sign constantly changes and 1  ∗  2π a(t)m ˜x 2 the integral 0 r1 (t) − 1+h(t)a(t)K dt ≈ −4.7968 < 0. By Corollary 3.5, we know that the periodic u 1 boundary solution (0, x2 ∗ (t)) of System (4.6) is globally asymptotically stable, which also can be seen from Figure 10, in which the orbits start from initial points (6, 1.5), (6, 20), (3, 2) and (5, 30), respectively.

1

30

0.5

25

0

20

−0.5

15

x2

4

+ sin t) −

−1

10

−1.5

5

−2

0

10

20

30

40

0

50

0

1

2

3 x1

t

Figure 9: The graph of g(t) = r1 (t) − for System (4.6).

8 9 10

⎪ ⎩

13

dx1 dt

6

Figure 10: Phase portrait of (4.6): Periodic boundary

solution (0, x∗2 (t)) of System (4.6) is globally asymptotically. Here, the orbits start from initial points (6, 1.5), (6, 20), (3, 2) and (5, 30), respectively.

dx2 dt

  

x1 (t) (1+sin t)x2 (t) − = x1 (t) (1 + 0.2 cos t) 1 − 100+20 cos t 1+0.1(1+sin t)2 x1 (t)    0.115(1+sin t)2 x1 (t) = x2 (t) 0.5 1 − x215(t) − (1.5 + 0.1 cos t) + 1+0.1(1+sin 2 t) x1 (t)

One can verify that 2π  r2 (t) − d(t) + 0

12

2

5

Example 4.6. [The periodic versus constant environment] Take r1 (t) = 1+0.2 cos t, r2 (t) = 0.5, K1 = 100 + 20 cos t, K2 = 15, a(t) = 1 + sin t, d(t) = 1.5 + 0.1 cos t, h(t) = 0.1(1 + sin t), e(t) = 0.115(1 + sin t), then we have the following 2π-periodic System (1.1)-(1.2) ⎧ ⎪ ⎨

11

a(t)m ˜∗ x

1+h(t)a(t)K1u

4

e(t)a(t)K1u 1 + a(t)h(t)K1u





dt

(4.7)

2π

13.8(1 + sin t)2 dt 1 + 12(1 + sin t)2 0 ≈ −2π + 5.1612 < 0. = −2π +

Therefore, according to the first part of Theorem 3.4, we can conclude that the periodic boundary solution (x∗1 (t), 0) of System (4.7) is globally asymptotically stable. It also can be seen from Figure 11, in which

14

1 2

the orbits start form initial points (120, 6), (120, 1), (60, 8) and (100, 4), respectively. Here, x∗1 (t) solves the following Logistic Equation:   dx1 1 = (1 + 0.2 cos t) 1 − x1 (t) x1 (t). dt 100 + 20 cos t

100

x1

x1

100

50

0

0

100

200

300

400

500

600

700

800

900

50

0

1000

0

50

100

150

200

8

8

6

6

4 2 0

5 6 7

8 9 10

350

400

450

500

300

350

400

450

500

2 0

100

200

300

400

500

600

700

800

900

0

1000

Figure 11: Periodic boundary solution (x∗1 (t), 0) of non-autonomous System (4.7) is globally asymptotically. The orbits start from (120, 6), (120, 1), (60, 8) and (100, 4), respectively.

4

300

4

0

50

100

150

200

t

3

250

t

x2

x2

t

250

t

Figure 12: The autonomous system (4.8) with its parameter values being the averages of the corresponding periodic parameters in System (4.7) is permanent. It indicates that non-autonomous model may suppress the permanence of its autonomous version. The initial points are chosen as: (120, 6), (120, 1), (60, 8), (100, 4) and (40, 3).

Although, System (4.7) with periodic coefficients is non-permanent as predator x2 goes extinct, the corresponding autonomous system (4.8) with its parameter values being the average of the corresponding periodic functions in System (4.7) is permanent. Let all parameters in System (4.7) be replaced by their averages, i.e., r1 = 1, r2 = 0.5, K1 = 100, K2 = 15, a = 1, d = 1.5, h = 0.1, e = 0.115, we have the following autonomous system ⎧    x1 (t) x2 (t) 1 ⎪ − = x (t) 1 − ⎨ dx 1 dt 100 1+0.1x1 (t) (4.8)    ⎪ ⎩ dx2 = x2 (t) 0.5 1 − x2 (t) − 1.5 + 0.115x1 (t) dt 15 1+0.1x1 (t) eaK1 14 We can check that r2 − d + 1+ahK = 101 > 0. From Theorem 3.1 in [15], we know that System (4.8) is 1 permanent. See Figure 12, in which the orbits start also from initial points (120, 6), (120, 1), (60, 8) and (100, 4), respectively.

11 12 13 14

15

It is very interesting that Example 4.6 indicates that the non-autonomous model may suppress the permanence of its autonomous version. Now, let us take the following constant coefficients given by Kang and Fewell in [15] K1 = r1 = K2 = 1, r2 = 0.25, h = 4, e = 0.9, a = 5. (4.9) Then we have the following autonomous system ⎧   5x2 (t) 1 ⎪ ⎨ dx dt = x1 (t) (1 − x1 (t)) − 1+20x1 (t) ,   ⎪ ⎩ dx2 = x2 (t) 0.25(1 − x2 (t)) − d + 4.5x1 (t) . dt 1+20x1 (t) 15

(4.10)

1 2 3

From the numerical simulations under this set of parameter values by varying d given by Kang and Fewell [15], System (4.10) has none, one , two, or three coexistence equilibria when r2 > d while the system has none, one, or three coexistence equilibria when r2 < d.

4

If we take periodic parameters

5

r1 (t) ≡ 1, r2 (t) ≡ 0.25, K1 (t) = 1 + 0.5 sin(wπt), K2 (t) = 1 + 0.5 sin(wπt), 6 7

8 9 10

where w > 0, and other parameters a, h, e and d also are chosen as periodic functions such that their averages on a periodic interval are the data in (4.9), then we have the following periodic system ⎧    x1 (t) a(t)x2 (t) 1 ⎪ − = x (t) 1 − ⎨ dx 1 dt 1+0.5 sin(wπt) 1+a(t)h(t)x1 (t) , (4.11)   

⎪ x2 (t) e(t)a(t)x1 (t) ⎩ dx2 = x2 (t) 0.25 1 − − d(t) + . dt 1+0.5 sin(wπt) 1+a(t)h(t)x1 (t) By numerical simulations, we can find that the non-autonomous system (4.11) may promote or suppress the permanence of the corresponding autonomous system (4.10). And we have the following three implications:

11 12 13 14 15 16 17 18

1. Take d = 0.09 < r2 . In this case, system (4.10) has no coexistence equilibrium, and it has three boundary equilibria: E00 = (0, 0) which is a source, E10 = (K1 , 0) = (1, 0) which is a saddle, and E01 = (0, K2 (1 − rd2 )) = (0, 0.64) which is globally asymptotically stable. However, the corresponding periodic system (4.11) with d(t) = 0.09(1 + cos(wπt)) , a(t) = 5(1 + cos(wπt)), h(t) = 4(1 + cos(wπt)), e(t) = 0.9(1 + cos(wπt)) and w = 2.5 is permanent. See Figure 13, in which the orbits start from (1, 0.5), (0.5, 1) and (0.3, 0.7), respectively. It indicates that the non-autonomous model may improve the permanence of its autonomous version.

19

Take d = 0.42 > r2 . In this case, system (4.10) has only one coexistence equilibrium. Since r2 + ≈ 0.4643 > d > r2 , by Theorem 3.1 of Kang and Fewell [15], system (4.10) is permanent. However, the periodic system (4.11) with d(t) = 0.42(1 + sin(wπt)) , a(t) = 5(1 + sin(wπt)), h(t) = 4(1 + sin(wπt)), e(t) = 0.9(1 + sin(wπt)) and w = 2 is non-permanent with its predator x2 goes extinct. See Figure 14, in which the orbits start from (2, 0.5), (0.2, 1) and (0.8, 0.2), respectively. It indicates that the non-autonomous model may suppress the permanence of its autonomous version.

20

2.

21

eaK1 1+ahK1

22 23 24 25 26 27 28 29 30 31 32

eaK1 3. At last, we take d = 0.48 > r2 + 1+ahK ≈ 0.4643. For the larger d, System (4.10) has no coexistence 1 equilibrium, and it has two boundary equilibria E00 = (0, 0), which is a saddle, and E10 = (K1 , 0) = (1, 0) which is globally asymptotically stable. However, the corresponding periodic system (4.11) with d(t) = 0.48(1 + cos(wπt)) , a(t) = 5(1 + sin(wπt)), h(t) = 4(1 + cos(wπt)), e(t) = 0.9(1 + sin(wπt)) and w = 40, is permanent. See Figure 15, in which the orbits start from (1, 0.5), (0.5, 1) and (0.3, 0.7), respectively. It indicates that the non-autonomous model may improve the permanence of its autonomous version.

33 34 35 36

Form the discussion above, we can conclude that for the permanent (or non-permanent) constant coefficient system, if the parameters are replaced by periodic functions, the corresponding periodic system may be non-permanent (or permanent, respectively).

37

16

1.5

2 1.5

x1

x1

1 0.5 0

1 0.5

0

100

200

300

400

500

600

700

800

900

0

1000

0

100

200

300

400

t

500

600

700

800

900

1000

600

700

800

900

1000

t

1.5

1

x2

x2

1 0.5

0.5 0

0

100

200

300

400

500

600

700

800

900

0

1000

0

100

200

t

Figure 13: Take d = 0.09 < r2 , then autonomous sys-

tem (4.10) is globally asymptotically stable at E0,1 = (0, 0.64), while the corresponding periodic system (4.11) is permanent. Here, the solutions starting from (1, 0.5), (0.5, 1) and (0.3, 0.7) tends to a positive periodic solution. It indicates that the non-autonomous model may improve the permanence of its autonomous version.

1

2 3

5

300

400

500

t

Figure 14: Take d = 0.42 > r2 , then System (4.10) is permanent, while the corresponding periodic system (4.11) is non-permanent with predator x2 goes extinct. Here, the initial points are chosen as (2, 0.5), (0.2, 1) and (0.8, 0.2). It indicates that the non-autonomous model may suppress the permanence of its autonomous version.

Conclusions and potential future work

In this paper, we perform a systematic analysis on the dynamics of the non-autonomous predator-prey system (1.1)-(1.2) with Holling Type II functional response and predator being generalist.

4 5 6 7 8 9

In Section 2, we concentrate on the dynamics of System (1.1)-(1.2) with time-dependent parameters. We obtain the sufficient condition (H1) of permanence and positive invariant set Γ of system (1.1)-(1.2) form Theorem 2.2. In Theorem 2.3 and Corollary 2.4, the extinct conditions of species x1 and x2 are given. In Theorem 2.6, the sufficient conditions (H1) combined with (2.6) are given to prove the global attractivity of positive solution by constructing Lyapunov function.

10 11 12 13 14 15 16 17

In Section 3, we devote to the dynamics of System (1.1)-(1.2) with periodic parameters. It is proved in Theorem 3.1 by using Brouwer fixed-point theorem. We show that if Condition (H1) holds then System (1.1)-(1.2) has at least one positive periodic solution and that if in addition Condition (2.6) holds then system (1.1)-(1.2) has unique positive periodic solution that is globally asymptotically stable. Also, we derive the sufficient Condition (H2) of the existence of a positive periodic solution from Theorem 3.3 by continuation theorem. Finally, we get the sufficient conditions for global attractivity of boundary periodic solutions (0, x∗2 (t)) and (x∗1 (t), 0) from Theorem 3.4 and Corollary 3.5.

18 19 20 21 22 23 24 25

In this paper, the parameter a(t), h(t), e(t) and d(t) are allowed to be nonnegative. This is different to some of literature work, for example [11, 13, 14], in which all parameters are bounded below by positive constants. The extinct conditions of integral forms of species x1 and x2 in the first part of Theorem 2.3 and Corollary 2.4, respectively, as well the globally attractive conditions of integral forms for boundary periodic solutions (0, x∗2 (t)) and (x∗1 (t), 0) given in the first part of Theorem 3.4 and Corollary 3.5, respectively. Our results reflect the effects of the long term prey-predator behaviors on the number of species, and are more reasonable than the conditions for the corresponding autonomous system and also

17

x1

1

0.5

0

0

50

100

150

200

250

300

350

400

450

500

300

350

400

450

500

t 1.5

x2

1 0.5 0

0

50

100

150

200

250

t

Figure 15: Take d = 0.48 > r2 , then autonomous system (4.10) is globally asymptotically stable at E10 = (K1 , 0) = (1, 0), while the corresponding periodic system (4.11) is permanent. Here, the solutions starting from (1, 0.5), (0.5, 1) and (0.3, 0.7) tends to a positive periodic solution. It indicates that the non-autonomous model may improve the permanence of its autonomous version.

1 2

the conditions given by supremum and infimum of parameters for non-autonomous system. In Section 4, numerical simulations are performed to verify our theoretical results.

3 4 5

There are some potential study directions to improve our results such that it is more realistic and applicable for real ecosystems.

6 7 8 9 10 11

1. The current results need r1 (t) and r2 (t) to be bounded below by positive constants. However, it is more realistic that not only a(t), h(t), e(t), d(t) but also r1 (t) and r2 (t) are nonnegative functions. To extend our current work, it would be interesting that the more reasonable conditions of integral forms are found for the boundedness, positive invariance, permanence, non-persistence and globally asymptotic stability.

12 13 14 15 16 17

2. From numerical simulations, we find that periodic system may promote or suppress permanence of its autonomous version with parameters being the averages of periodic parameters (see Example 4.6). It would be interesting to explore when and how periodic environment changes the permanence of corresponding autonomous system, especially, how the amplitudes and the periods of periodic parameters influence the permanence of non-autonomous system.

18

19

6

20

Proof of Lemma 2.1

21 22

Proofs

Proof. We only need to show the first part of the lemma. Multiple the both sides of the inequality u (t) ≤ p(t)u(t) − q(t)u2 (t) by u−2 (t), we get   t   t  d −1 u (t) exp p(τ )dτ ≥ q(t) exp p(τ )dτ . dt t0 t0

18

1

Integrating from t0 to t, we have  t  u−1 (t) exp p(τ )dτ ≥

u−1 (t0 ) +



t0

= 2

t→+∞

5

Proof of Theorem 2.2

6

Proof. 1.

8 9 10

−1

u

(t0 ) + inf

t∈R



q(t) p(t) q(t) p(t)

 

12



t t0



15 16 17 18

19

t



p(τ )dτ

ds

 −1 ,

t0

t∈R

we know that both the nonnegative and positive cones of R2 are positively invariant for (1.1)-(1.2). 2. We prove that Γ is a positively invariant set for (1.1)-(1.2). Let (x1 (t), x2 (t)) be the solution of (1.1)(1.2) with (x1 (t0 ), x2 (t0 )) ∈ Γ. From the equation (1.1) and the positivity of the solutions of (1.1)-(1.2), we have t ≥ t0 .

(6.1)

which implies by Lemma 2.1 and 0 < x1 (t0 ) ≤ K1u that x1 (t) ≤ K1u holds for all t ≥ t0 . Then, the equation (1.2) yields

≤ r2 (t) − d(t) +

e(t)a(t)K1u 1+h(t)a(t)K1u



x2 (t) −

r2 (t) 2 K2 (t) x2 (t),

t ≥ t0 ,

(6.2)

and hence, by Lemma 2.1 and 0 < x2 (t0 ) ≤ Mx2 , x2 (t) ≤ Mx2 for t ≥ t0 . It follows dx1 r1 (t) 2 ≥ (r1 (t) − a(t)Mx2 )x1 (t) − x (t) dt K1 (t) 1

14

p(τ )dτ t0

exp



s

From the following expressions     t   a(s)x2 (s) x1 (s) r1 (s) 1 − − ds , x1 (t) = x1 (t0 ) exp K1 (s) 1 + h(s)a(s)x1 (s) t0     t   e(s)a(s)x1 (s) x2 (s) − d(s) + ds , r2 (s) 1 − x2 (t) = x2 (t0 ) exp K2 (s) 1 + h(s)a(s)x1 (s) t0

dx2 dt 13

ds

t∈R

dx1 r1 (t) 2 ≤ r1 (t)x1 (t) − x (t), dt K1 (t) 1 11

t0

p(s) exp

t∈R

4

7

t∈R



p(τ )dτ

p(t) p(t) which implies lim sup u(t) ≤ sup p(t) q(t) . Moreover, if 0 < u(t0 ) ≤ sup q(t) , we have u(t) ≤ sup q(t) , t ∈

[t0 , +∞). The rest proof of lemma is committed.

3

u−1 (t0 ) + inf



s

q(s) exp

t0

≥



t

(6.3)

and therefore for all t ≥ t0 , x1 (t) ≥ mx1 by Lemma 2.1 and x1 (t0 ) ≥ mx1 > 0. Then, furthermore,   dx2 e(t)a(t)mx1 r2 (t) 2 x2 (t) − ≥ r2 (t) − d(t) + x (t) (6.4) dt 1 + h(t)a(t)mx1 K2 (t) 2 which implies, from Lemma 2.1 and x2 (t0 ) ≥ mx2 > 0, that for all t ≥ t0 , x2 (t) ≥ mx2 . Therefore, Γ is positively invariant with respect to (1.1)-(1.2). 3. Assume that Condition (H1) holds. We prove that system (1.1)-(1.2) is permanent. More precisely, we want to show that for any solution (x1 (t), x2 (t)) of (1.1)-(1.2) with xi (t0 ) > 0, i = 1, 2, (i) lim sup x1 (t) ≤ K1u ; t→+∞

19

 (ii) lim sup x2 (t) ≤ Mx2 if inf r2 (t) − d(t) + t∈R

e(t)a(t)K1u 1+h(t)a(t)K1u

2

 (iii) lim inf x1 (t) ≥ mx1 if inf r2 (t) − d(t) + t∈R

e(t)a(t)K1u 1+h(t)a(t)K1u

3

(iv) lim inf x2 (t) ≥ mx2 if (H1) holds.

1

t→+∞

4

5 6

t→+∞

 > 0; 

> 0 and inf (r1 (t) − a(t)Mx2 ) > 0; t∈R

t→+∞

By Condition (H1), we can choose  > 0 sufficiently small such that mx1 −  > 0 and   e(t)a(t)(mx1 − ) > 0. inf (r1 (t) − a(t)(Mx2 + )) > 0, inf r2 (t) − d(t) + t∈R t∈R 1 + h(t)a(t)(mx1 − ) Let (x1 (t), x2 (t)) be any solution of (1.1)-(1.2) with xi (t0 ) > 0, i = 1, 2. By Lemma 2.1, the inequality (6.1) implies that lim sup x1 (t) ≤ K1u . It follows that there exists T0 > t0 such that for t > T0 , x1 (t) ≤ K1u + . t→+∞

7

8

Then, the equation (1.2) implies that   dx2 e(t)a(t)(K1u + ) r2 (t) 2 ≤ r2 (t) − d(t) + x2 (t) − x (t), t ≥ t0 . u dt 1 + h(t)a(t)(K1 + ) K2 (t) 2   e(t)a(t)K1u Since inf r2 (t) − d(t) + 1+h(t)a(t)K > 0, we have by Lemma 2.1 that u 1 t∈R

 e(t)a(t)(K1u +) 2 (t) lim sup x2 (t) ≤ sup K r2 (t) r2 (t) − d(t) + 1+h(t)a(t)(K1u +) t→+∞ t∈R

   e(t)a(t)K1u K2 (t) e(t)a(t) 2 (t) ≤ sup K r + sup . (t) − d(t) + · u u 2 r2 (t) 1+h(t)a(t)K r2 (t) 1+h(t)a(t)K 1

t∈R

9 10

Thus, by the boundedness of

K2 (t) r2 (t)

·

e(t)a(t) 1+h(t)a(t)K1u

1

t∈R

and the arbitrariness of , we have lim sup x2 (t) ≤ Mx2 .

It follows that there exists T1 > t0 such that for t > T1 , x2 (t) ≤ Mx2 +  and

t→+∞

dx1 r1 (t) 2 ≥ (r1 (t) − a(t)(Mx2 + ))x1 (t) − x (t). dt K1 (t) 1 11

Since inf (r1 (t) − a(t)(Mx2 + )) > 0, we have again by Lemma 2.1 that t∈R

lim inf x1 (t) t→+∞

12 13

14

K1 (t) (r1 (t) t∈R r1 (t) K1 (t) (r1 (t) inf t∈R r1 (t)

≥ inf

− a(t)(Mx2 + ))

≥

1 (t) − a(t)Mx2 ) −  sup K r1 (t)

which implies, by the arbitrariness of , that lim inf x1 (t) ≥ mx1 . Then, there exists T2 > t0 such that t→+∞

for t > T2 , x1 (t) ≤ mx1 −  and   dx2 e(t)a(t)(mx1 − ) r2 (t) 2 ≥ r2 (t) − d(t) + x2 (t) − x (t) dt 1 + h(t)a(t)(mx1 − ) K2 (t) 2

 e(t)a(t)(mx1 −) > 0 and Lemma 2.1, we have Thus, from inf r2 (t) − d(t) + 1+h(t)a(t)(m x1 −) t∈R

 e(t)a(t)(mx1 −) 2 (t) r lim inf x2 (t) ≥ inf K (t) − d(t) + 2 r (t) 1+h(t)a(t)(m −) x1 t→+∞ t∈R 2

   e(t)a(t)mx1 e(t)a(t) 2 (t) 2 (t) r −  sup K (t) − d(t) + ≥ inf K 2 r2 (t) 1+h(t)a(t)mx r2 (t) · 1+h(t)a(t) . 1

t∈R

15 16

t∈R

t∈R

By the arbitrariness of , lim inf x2 (t) ≥ mx2 . t→+∞

4.

From the arguments above, we can conclude that the set Γ defined by Γ := {(x1 , x2 ) ∈ R2 | mx1 −  ≤ x1 ≤ K1u + , mx2 −  ≤ x2 ≤ Mx2 + },

17 18

is an ultimately bounded region of (1.1)-(1.2). Here,  > 0 is sufficiently small such that mx1 −  > 0 and mx2 −  > 0. 20

1

Proof of Theorem 2.3 Proof. Let (x1 (t), x2 (t)) be any solution of (1.1)-(1.2) with xi (t0 ) > 0, i = 1, 2. 1. By Lemma 2.1 and equation (1.1), lim sup x1 (t) ≤ K1u . It follows that there exists T0 > t0 such that t→+∞

for t ≥ T0 , x1 (t) ≤ K1u + , where  > 0 sufficiently small. Then, form the equation (1.2), we have   e(t)a(t)(K1u + ) , t ≥ T0 . x2 (t) ≤ x2 (t) r2 (t) − d(t) + 1 + h(t)a(t)(K1u + ) 2 3

 +∞ By our assumptions we know that T0 p(t)dt = −∞. Thus,

  t from the inequality above, we get x2 (t) ≤ x2 (T0 ) exp T0 p(τ )dτ and lim x2 (t) = 0.

Let p(t) = r2 (t) − d(t) +

e(t)a(t)(K1u +) 1+h(t)a(t)(K1u +) .

t→+∞

4

Assume a(t) and h(t) are bounded below by positive constants, and inf (r2 (t) − d(t)) > 0. We prove

5

2.

6

1 (t) that prey x1 will be extinct for system (1.1)-(1.2). By m ˜ x2 > sup ra(t) , there exists  : 0 <  < m ˜ x2 , such

7

that sup(r1 (t) − a(t)(m ˜ x2 − )) ≤ 0. From the equation (1.2), we have that for t ≥ t0 ,

t≥0

t≥0

t≥0

x2 (t) ≥ (r2 (t) − d(t))x2 (t) −

r2 (t) 2 x (t), K2 (t) 2

which implies by inf (r2 (t) − d(t)) > 0 and Lemma 2.1 that t≥0

lim inf x2 (t) ≥ inf t→+∞

8 9

x1 (t) 1+h(t)a(t)x1 (t)

[(r1 (t) − a(t)(m ˜ x2 − )) +

(i) (ii)

r1 (t) K1 (t) (h(t)a(t)K1 (t)

− 1)x1 (t) −

r1 (t)h(t)a(t) 2 x1 (t) K1 (t)

1 (t) If sup(h(t)a(t)K1 (t)) < 1 and m ˜ x2 > sup ra(t) , it is clear that

t≥0

11

K2 (t) (r2 (t) − d(t)) = m ˜ x2 . r2 (t)

It follows that there exists T1 > t0 such that for t > T1 , x2 (t) ≥ m ˜ x2 − . Then, by the equation (1.1), we have

 a(t)(m ˜ x2 −) r1 (t) dx1 ≤ x (t) r (t) − x (t) − 1 1 dt K1 (t) 1 1+h(t)a(t)x1 (t) =

10

t≥0

t≥0

If inf (h(t)a(t)K1 (t)) > 1 and m ˜ x2 > sup t≥0

t≥0

 Δ=

r1 (t) (h(t)a(t)K1 (t) − 1) K1 (t)

dx1 dt

r1 (t)[h(t)a(t)K1 (t)+1]2 , 4h(t)a2 (t)K1 (t)

dx1 dt

 , t > T1 .

< 0 for t > T1 and lim x1 (t) = 0. t→∞

then for each t > T1 ,

2 + 4(r1 (t) − a(t)(m ˜ x2 − ))

r1 (t)h(t)a(t) < 0. K1 (t)

< 0 for t > T1 and lim x1 (t) = 0.

12

It follows that

13

Proof of Corollary 2.4

14

Proof. From the proof of Theorem 2.3, we have lim sup x1 (t) ≤ K1u and lim inf x2 (t) ≥ m ˜ x2 . It follows

t→∞

t→+∞

t→+∞

15 16

that there exists T > t0 such that for t > T , x1 (t) ≤ K1u + , x2 (t) ≥ m ˜ x2 − , where  > 0 is sufficiently small. Thus, by the equation (1.1), we have

 a(t)(m ˜ x2 −) x1 (t) ≤ x1 (t) r1 (t) − 1+h(t)a(t)(K , t > T. u +) 1

21

1

Let q(t) = r1 (t) −

a(t)(m ˜ x2 −) 1+h(t)a(t)(K1u +) .

By (2.5), we have 

x1 (t) ≤ x1 (T ) exp

 +∞ T

q(t)dt = −∞. Therefore,



t

q(τ )dτ

→0

as

t→ + ∞.

T 2

The proof is complete.

3

Proof of Theorem 2.6

4 5

Proof. Let (x1 (t), x2 (t)) be any solution of (1.1)-(1.2) with x1 (t0 ) > 0, x2 (t0 ) > 0. By (2.6), we can choose  : 0 <  < min{mx1 , mx2 } sufficiently small such that ⎧ h(t)a2 (t)Mx2 +) e(t)a(t) ⎨ r1 (t) − K1 (t) (1+h(t)a(t)(mx1 −))2 − (1+h(t)a(t)(mx1 −))2 > 0, (6.5) a(t) ⎩ r2 (t) − > 0. K2 (t)

6 7

1+h(t)a(t)(mx1 −)

Since Γ is an ultimately bounded region of system (1.1)-(1.2), there exists T1 > 0 such that (x1 (t), x2 (t)) ∈ Γ and (x∗1 (t), x∗2 (t)) ∈ Γ for all t ≥ t0 + T1 . Define V (t) = | ln(x1 (t)) − ln(x∗1 (t))| + | ln(x2 (t)) − ln(x∗2 (t))|,

8

(6.6)

By a direct calculation of the right derivative D+ V (t) of V (t) along system (1.1)-(1.2), we have D+ V (t)

=

≤

9

t ≥ t0 + T0 .

r1 (t) a(t)(x2 (t) − x∗2 (t)) |x1 (t) − x∗1 (t)| − sgn{x1 (t) − x∗1 (t)} K1 (t) 1 + h(t)a(t)x∗1 (t) 2 ∗ h(t)a (t)x2 (t)|x1 (t) − x1 (t)| r2 (t) + − |x2 (t) − x∗2 (t)| ∗ (1 + h(t)a(t)x1 (t))(1 + h(t)a(t)x1 (t)) K2 (t) e(t)a(t)(x1 (t) − x∗1 (t)) +sgn{x2 (t) − x∗2 (t)} (1 + h(t)a(t)x∗1 (t))(1 + h(t)a(t)x1 (t))   r1 (t) h(t)a2 (t)(Mx2 + ) e(t)a(t) − − |x1 (t) − x∗1 (t)| − K1 (t) (1 + h(t)a(t)(mx1 − ))2 (1 + h(t)a(t)(mx1 − ))2   r2 (t) a(t) − |x2 (t) − x∗2 (t)|. − K2 (t) 1 + h(t)a(t)(mx1 − ) −

Then from (6.5) there exists a positive constant ν > 0 such that D+ V (t) ≤ −ν(|x1 (t) − x∗1 (t)| + |x2 (t) − x∗2 (t)|)

t ≥ t0 + T0 .

10

Thus, one can use Lemma 2.5 to verify that lim (|x1 (t) − x∗1 (t)| + |x2 (t) − x∗2 (t)|) = 0 (see the proof of

11

Theorem 2.4 in [11]). The proof is complete.

12

Proof of Theorem 3.1

13

Proof. Define a Poincare mapping σ : R2 →R2 as follows

t→+∞

σ((x01 , x02 )) = (x1 (t0 + ω, t0 , (x01 , x02 )), x2 (t0 + ω, t0 , (x01 , x02 )), 14 15 16 17 18

(x01 , x02 ) ∈ R2 ,

where (x1 (t, t0 , (x01 , x02 )), x2 (t, t0 , (x01 , x02 )) denotes the solution of (1.1)-(1.2) through (t0 , (x01 , x02 )). Then σ(Γ) ⊂ Γ by the positive invariant property of Γ for (1.1)-(1.2). The continuity of σ can be guaranteed by the continuity of solution of (1.1)-(1.2) with respect to initial values. Note that Γ is a bounded, closed, convex set in R2 . Therefore, by Brouwer fixed point theorem, the operator σ has at least one fixed point (x∗1 (t), x∗2 (t)) in Γ, which is a positive ω-periodic solution of system (1.1)-(1.2). The proof is complete. 22

1

Proof of Theorem 3.3

2

Proof. Let xi (t) = exp{yi (t)}, i = 1, 2, and rewrite system (1.1)-(1.2) as follows dy1 r1 (t) a(t) exp{y2 (t)} = r1 (t) − exp{y1 (t)} − , dt K1 (t) 1 + h(t)a(t) exp{y1 (t)}

(6.7)

dy2 r2 (t) e(t)a(t) exp{y1 (t)} = r2 (t) − exp{y2 (t)} − d(t) + . dt K2 (t) 1 + h(t)a(t) exp{y1 (t)}

(6.8)

Thus, if (y1 (t), y2 (t)) is a ω-periodic solution of system (6.7)-(6.8), then (x1 (t), x2 (t)) = (exp{y1 (t)}, exp{y2 (t)}) 3 4

is a positive ω-periodic solution of (1.1)-(1.2). Take X = Y = {(y1 , y2 )T ∈ C(R, R2 )| yi (t + ω) = yi (t), i = 1, 2, t ∈ R},

5

endowed with the norm ||(y1 , y2 )|| = max |y1 (t)| + max |y2 (t)|. Then (X, || · ||) and (Y, || · ||) are both t∈[0,ω]

6 7

Banach spaces. Let      y1 y1 L = , y2 y2  N

8 9 10

y1 y2



⎡ =⎣

 P

y1 y2

t∈[0,ω]



 =Q

y1 y2



 =

r1 (t) −

r1 (t) K1 (t)

exp{y1 (t)} −

r2 (t) −

r2 (t) K2 (t)

exp{y2 (t)} − d(t) +

1 ω

ω

0  1 ω ω 0

y1 (t)dt y2 (t)dt

 ,

a(t) exp{y2 (t)} 1+h(t)a(t) exp{y1 (t)} e(t)a(t) exp{y1 (t)} 1+h(t)a(t) exp{y1 (t)}

⎤ ⎦,



y1 y2

 ∈ X.

Then system (6.7)-(6.8) is equivalent to the operator equation Ly = N y, y = (y1 , y2 )T ∈ X. Thus, we only need to verify the conditions in Lemma 3.2 for the operator equation and obtain the existence of at least one ω-periodic solution of system (6.7)-(6.8). It is easy to see that KerL = {(y1 , y2 ) ∈ X| (y1 , y2 ) = (c1 , c2 ) ∈ R2 }, ω ImL = {(y1 , y2 ) ∈ X| yi (t)dt, i = 1, 2} 0

11 12 13

and dim KerL = 2 = codim ImL. Since ImL is closed in Y , L is a Fredholm mapping of index zero. Obviously, P (= Q) is a continuous projection such that ImP = KerL, ImL = KerQ = Im(I − Q). KP−1 : ImL→DomL ∩ KerP exists and is given by  ωt    t y (s)ds − ω1 0 0 y1 (s)dsdt y1 0 1 −1 KP = t . ωt y2 y (s)ds − 1 y (s)dsdt 0

14 15 16

2

ω

0

0

2

Since QN and KP−1 (I − Q)N are continuous, one can verify that N is L-compact on an any closed bounded set in X. Let y = (y1 , y2 ) ∈ X be a solution of Ly = λN y for a certain λ ∈ (0, 1), then   dy1 r1 (t) a(t) exp{y2 (t)} = λ r1 (t) − exp{y1 (t)} − , (6.9) dt K1 (t) 1 + h(t)a(t) exp{y1 (t)}   dy2 r2 (t) e(t)a(t) exp{y1 (t)} = λ r2 (t) − exp{y2 (t)} − d(t) + . dt K2 (t) 1 + h(t)a(t) exp{y1 (t)} 23

(6.10)

1

By integrating over [0, ω], we have  ω r1 (t) a(t) exp{y2 (t)} exp{y1 (t)} + dt, ωA(r1 ) = K1 (t) 1 + h(t)a(t) exp{y1 (t)} 0 ωA(r2 ) =

2



0

 r2 (t) e(t)a(t) exp{y1 (t)} exp{y2 (t)} + d(t) − dt. K2 (t) 1 + h(t)a(t) exp{y1 (t)}

(6.12)

Then we have that by (6.9) and (6.11), ω 0

3

ω

(6.11)

|y1 (t)|dt

ω  ω  r1 (t) ≤ λ 0 r1 (t)dt + λ 0 K exp{y1 (t)} + 1 (t) < 2ωA(r1 ),

a(t) exp{y2 (t)} 1+h(t)a(t) exp{y1 (t)}

 dt

(6.13)

and that by (6.10) and (6.12), ω 0

|y2 (t)|dt

≤ <

ω  ω  r2 (t) λ 0 r2 (t)dt + λ 0 K exp{y2 (t)} + d(t) − 2 (t) 2ωA(r2 ).

e(t)a(t) exp{y1 (t)} 1+h(t)a(t) exp{y1 (t)}

 dt

(6.14)

Denote yi (ξi ) = min yi (t), t∈[0,ω]

4



ω

0

r1 (t) exp{y1 (ξ1 )}dt = ωA(r1 /K1 ) exp{y1 (ξ1 )}, K1 (t)

which implies that A(r1 ) = ln S1 := l11 . A(r1 /K1 )

(6.15)

|y1 (t)|dt ≤ ln S1 + 2ωA(r1 ) := H11 .

(6.16)

y1 (ξ1 ) ≤ ln 6

t∈[0,ω]

From (6.11), we have ωA(r1 ) ≥

5

yi (ηi ) = max yi (t), i = 1, 2.

Thus by (6.13) and (6.15) we have y1 (t) ≤ y1 (ξ1 ) +

ω 0

7

Then from (6.12) and (6.16), we have  ω r2 (t) e(t)a(t)S1 exp{2ωA(r1 )} exp{y2 (ξ2 )} + d(t) − dt ωA(r2 ) ≥ K2 (t) 1 + h(t)a(t)S1 exp{2ωA(r1 )} 0   eaS1 exp{2ωA(r1 )} , = ωA(r2 /K2 ) exp{y2 (ξ2 )} + ωA(d) − ωA 1 + haS1 exp{2ωA(r1 )}

8

which implies that y2 (ξ2 ) ≤ ln

9

A r2 − d +

eaS1 exp{2ωA(r1 )} 1+haS1 exp{2ωA(r1 )}

A(r2 /K2 )

 = ln S2 := l21 .

(6.17)

Thus by (6.14) and (6.17) we get y2 (t) ≤ y2 (ξ2 ) +

ω 0

|y2 (t)|dt ≤ ln S2 + 2ωA(r2 ) := H21 .

24

(6.18)

1

Further, from (6.12) together with (6.18), we have  ω r1 (t) exp{y1 (t)} + a(t) exp{y2 (t)} dt ωA(r1 ) ≤ K1 (t) 0 ≤ ωA(r1 /K1 ) exp{y1 (η1 )} + ωA(a)S2 exp{2ωA(r2 )},

2

which implies that y1 (η1 ) ≥ ln

3

A(r1 ) − A(a)S2 exp {2ωA(r2 )} = ln S3 := l12 . A(r1 /K1 )

Thus, (6.13) and (6.19) yield that y1 (t) ≥ y1 (η1 ) −

4

(6.19)

ω 0

|y1 (t)|dt ≥ ln S3 − 2ωA(r1 ) := H12 ,

(6.20)

which, together with (6.16), reduces to max |y1 (t)| ≤ max{|H11 |, |H12 |} := B1 . t∈[0,ω]

5

6

Moreover, (6.12) and (6.20) yield  ω r2 (t) e(t)a(t)S3 exp{−2ωA(r1 )} exp{y2 (η2 )} + d(t) − dt ωA(r2 ) ≤ K2 (t) 1 + h(t)a(t)S3 exp{−2ωA(r1 )} 0   eaS3 exp{−2ωA(r1 )} = ωA(r2 /K2 ) exp{y2 (η2 )} + ωA(d) − A 1 + haS3 exp{−2ωA(r1 )} which implies that y2 (η2 ) ≥ ln

7

A r2 − d +

eaS3 exp{−2ωA(r1 )} 1+haS3 exp{−2ωA(r1 )}



A(r2 /K2 )

= ln S4 := l22 .

Thus, (6.14) together with (6.21) yields that ω y2 (t) ≥ y2 (η2 ) − |y2 (t)|dt ≥ ln S4 − 2ωA(r2 ) := H22 ,

(6.21)

(6.22)

0

8 9 10 11 12 13 14

15 16

which, together with (6.18), reduces to max |y2 (t)| ≤ max{|H21 |, |H22 |} := B2 . t∈[0,ω]

Now, take B = B1 +B2 +B3 , where B3 > 0 is sufficiently large such that B3 > |l11 |+|l12 |+|l21 |+|l22 |. Define Ω = {(y1 , y2 )T ∈ X| ||(y1 , y2 )|| < B}. Then for each λ ∈ (0, 1), every solution y = (y1 , y2 )T of Ly = λN y satisfies y ∈ / ∂Ω. In order to verify the second condition of Lemma 3.2, we need to consider the following algebraic equations a(t) exp{y2 } μ ω dt = 0, (6.23) A(r1 ) − A(r1 /K1 ) exp{y1 } − ω 0 1 + h(t)a(t) exp{y1 } μ ω e(t)a(t) exp{y1 } A(r2 ) − A(r2 /K2 ) exp{y2 } − A(d) + dt = 0, (6.24) ω 0 1 + h(t)a(t) exp{y1 } where (y1 , y2 ) ∈ R2 and μ ∈ [0, 1] is a parameter. Similar to the arguments above, one can easily check that any solution (y1∗ , y2∗ ) with μ ∈ [0, 1] satisfies l12 ≤ y1∗ ≤ l11 ,

l22 ≤ y2∗ ≤ l21 .

25

(6.25)

1

2 3 4

Then for each y = (y1 , y2 )T ∈ ∂Ω ∩ KerL = ∂Ω ∩ R2 , QN y = 0, i.e., ⎤ ⎡ ω a(t) exp{y2 }     A(r1 ) − A(r1 /K1 ) exp{y1 } − ω1 0 1+h(t)a(t) exp{y1 } dt 0 y1 ⎦ ⎣ = . = QN  ω e(t)a(t) exp{y1 } 0 y2 A(r2 − d) − A(r2 /K2 ) exp{y2 } + ω1 0 1+h(t)a(t) exp{y1 } dt The first part of condition (2) of Lemma 3.2 is valid. Now, we verify that the Brouwer degree deg(JQN, Ω ∩ KerL, 0) = 0. For this purpose, define a homotopy Fμ y = μQN y + (1 − μ)Gy,

5

where

 Gy =

6 7 8

μ ∈ [0, 1], y = (y1 , y2 )T ∈ X,

A(r1 ) − A(r1 /K1 ) exp{y1 } A(r2 − d) − A(r2 /K2 ) exp{y2 }

 .

Since any solution of (6.23)-(6.24) satisfies (6.25), we know that 0 ∈ / Fμ (∂Ω ∩ KerL) for μ ∈ [0, 1]. That is, the Brouwer degree deg(Fμ , Ω ∩ KerL, 0) is well defined. By the invariance property of homotopy, we have deg(Fμ , Ω ∩ KerL, 0) = deg(QN, Ω ∩ KerL, 0) = deg(G, Ω ∩ KerL, 0).

9 10

Notice that the algebraic equation Gy = 0 has a unique solution in R2 and J = I since ImQ = KerL. Thus, deg(JQN, Ω ∩ KerL, 0) = deg(G, Ω ∩ KerL, 0) = 0.

14

By now, all conditions of Lemma 3.2 are verified. Therefore, Ly = N y has at least one solution in ¯ , i.e., system (6.7)-(6.8) has at least one ω-periodic solution (y ∗ (t), y ∗ (t)) ∈ DomL ∩ Ω. ¯ Hence, DomL ∩ Ω 1 2 ∗ ∗ ∗ ∗ (x1 (t), x2 (t)) = (exp{y1 (t)}, exp{y2 (t)}) is a strictly positive ω-periodic solution of system (1.1)-(1.2). The proof is complete.

15

Proof of Theorem 3.4

11 12 13

16 17

Proof. 1. Let (x1 (t), x2 (t)) be any solution of system (1.1)-(1.2) with x1 (t0 ) > 0 and x2 (t0 ) > 0. By the equation (1.1), Lemma 2.1 implies lim sup x1 (t) ≤ K1u . Clearly, (3.3) implies that (2.4) holds. Thus, by t→+∞

18 19

20

Theorem 2.3, we know that lim x2 (t) = 0. Therefore, we only need to prove lim |x1 (t) − x∗1 (t)| = 0. t→+∞

t→+∞

By (3.3), there exists  > 0 such that ω r2 (t) − d(t) +

 e(t)a(t)K1u + (ρ + ρ ) dt = 0. 1 2 1 + h(t)a(t)K1u 0  ω e(t)a(t)K1u dt, where In fact,  = − ω(ρ11+ρ2 ) 0 r2 (t) − d(t) + 1+h(t)a(t)K u 1

   r2 (t)  , ρ1 = sup a(t) − K2 (t)  t∈R+

ρ2 = sup

t∈R+

e(t)a(t) . 1 + h(t)a(t)K1u

Define V (t) = | ln(x1 (t)) − ln(x∗1 (t))| + | ln x2 (t)|, 26

t > 0.

(6.26)

1

Since lim sup x1 (t) ≤ K1u and lim x2 (t) = 0, there exists T0 > t0 such that for all t > T0 , 0 < x1 (t) ≤ t→+∞

t→+∞

2 3

K1u +  and x2 (t) < . Then we have by calculating the right derivative of V (t) along system (1.1)-(1.2) that D+ V (t)

≤

≤

≤ 4

r1 (t) a(t)x2 (t) |x1 (t) − x∗1 (t)| + K1 (t) 1 + h(t)a(t)x1 (t) r2 (t) e(t)a(t)x1 (t) x2 (t) − d(t) + +r2 (t) − K2 (t) 1 + h(t)a(t)x1 (t)    r1l r2 (t)  ∗  − u |x1 (t) − x1 (t)| + a(t) − x2 (t) K1 K2 (t)  e(t)a(t)(K1u + ) +r2 (t) − d(t) + 1 + h(t)a(t)(K1u + )

−

−

r1l e(t)a(t)K1u |x1 (t) − x∗1 (t)| + r2 (t) − d(t) + + (ρ1 + ρ2 ), u K1 1 + h(t)a(t)K1u

t > T0 .

Integrating on both sides from T0 to t, we get by (6.26) that r1l K1u

t

|x1 (s) − x∗1 (s)|ds T   t0 e(s)a(s)K1u ds ≤ V (T0 ) + T0 r2 (s) − d(s) + 1+h(s)a(s)K u + (ρ1 + ρ2 ) 1 < +∞, t > T0 .

V (t) +

5

which implies that |x1 (t) − x∗1 (t)| ∈ L1 ([T0 , +∞)). The boundedness of x∗1 (t) and the ultimate boundedness of x1 (t) imply that both x∗1 (t) and x1 (t) have bounded derivatives for t > T0 (from the equations satisfied by them). Then |x1 (t) − x∗1 (t)| is uniformly continuous on [T0 + ∞). Thus, by Lemma 2.5, we have lim |x1 (t) − x∗1 (t)| = 0. t→+∞

2. Let (x1 (t), x2 (t)) be any solution of system (1.1)-(1.2) with x1 (t0 ) > 0 and x2 (t0 ) > 0. By Lemma 2.1 , we have lim sup x1 (t) ≤ K1u by the equation (1.1) and lim inf x2 (t) ≥ m ˜ ∗x2 > 0 by the equation (1.2). t→+∞

t→+∞

˜ ∗x2 − , where  > 0 is Then there exists T0 > t0 such that for all t > T0 , x1 (t) < K1u +  and x2 (t) ≥ m sufficiently small. Define V (t) = | ln x1 (t)| + | ln(x2 (t)) − ln(x∗2 (t))|. 6

By calculating the right derivative of V (t) along system (1.1)-(1.2), we have D+ V (t)

≤

≤

r1 (t) a(t)x2 (t) x1 (t) − K1 (t) 1 + h(t)a(t)x1 (t) e(t)a(t)x1 (t) r2 (t) |x2 (t) − x∗2 (t)| + − K2 (t) 1 + h(t)a(t)x1 (t) r2 (t) r1 (t) − |x2 (t) − x∗2 (t)| − x1 (t) K2 (t) K1 (t) a(t)(m ˜ ∗x2 − ) e(t)a(t)x1 (t) +r1 (t) − + , 1 + h(t)a(t)x1 (t) 1 + h(t)a(t)x1 (t)

r1 (t) −

27

T > T0 .

1

(1)

If r1 (t) +

e(t)a(t)K1u 1+h(t)a(t)K1u

 Here, ν1 = min

3

(2)

r1l r2l K1u , K2u



a(t)m ˜∗ x2 1+h(t)a(t)K1u , t

∈ [0, ω], then

r2 (t) r1 (t) |x2 (t) − x∗2 (t)| − x1 (t) K2 (t) K1 (t) a(t)(m ˜ ∗x2 − ) e(t)a(t)(K1u + ) +r1 (t) − + u 1 + h(t)a(t)(K1 + ) 1 + h(t)a(t)(K1u + ) r2 (t) r1 (t) < − |x2 (t) − x∗2 (t)| − x1 (t) K2 (t) K1 (t) < −ν1 (|x2 (t) − x∗2 (t)| + x1 (t)), T > T0 .

D+ V (t)

2

< ≤

−

. Thus, lim |(x1 (t), x2 (t)) − (0, x∗2 (t))| = 0. t→+∞

Note a(t)(m ˜ ∗x2 − ) r1 (t) x1 (t) − K1 (t) 1 + h(t)a(t)x1 (t)     r1 (t) r1 (t)h(t)a(t) 2 1 ∗ (r1 (t) − a(t)(m x1 (t) − x1 (t) . ˜ x2 − )) + r1 (t)h(t)a(t) − = 1 + h(t)a(t)x1 (t) K1 (t) K1 (t) r1 (t) −

4 5 6 7

Thus, if either (i) r1 (t) − a(t)m ˜ x2 < 0, h(t)a(t)K1 (t) < 1, t ∈ [0, ω]; or 2 1 (t)+1] (ii) h(t)a(t)K1 (t) > 1, m ˜ x2 > r1 (t)[h(t)a(t)K , t ∈ [0, ω], 4h(t)a2 (t)K1 (t) then r1 (t) −

8

a(t)(m ˜ ∗x2 − ) r1 (t) x1 (t) − < 0, t > T0 , K1 (t) 1 + h(t)a(t)x1 (t)

and lim x1 (t) = 0 by Theorem 2.3. Thus, there exists T1 > T0 such that t→+∞

D+ V (t) ≤ − 9

r2l |x2 (t) − x∗2 (t)|, t > T1 . K2u

It follows lim x2 (t) = x∗2 (t), and hence lim |(x1 (t), x2 (t)) − (0, x∗2 (t))| = 0. The proof is complete. t→+∞

t→+∞

10

Proof of Corollary 3.5

11

From the proof of the second part of Theorem 3.4, we have r2 (t) r1 (t) |x2 (t) − x∗2 (t)| − x1 (t) K2 (t) K1 (t) a(t)(m ˜ ∗x2 − ) e(t)a(t)x1 (t) +r1 (t) − + , T > T0 . 1 + h(t)a(t)(K1u + ) 1 + h(t)a(t)x1 (t)  ω a(t)m ˜∗ x2 dt < 0, we know lim x1 (t) = 0 by Theorem 2.3. Thus, Since 0 r1 (t) − 1+h(t)a(t)K u D+ V (t)

12

≤

−

1

t→+∞

D+ V (t) ≤ − 13

(6.27)

r2l |x2 (t) − x∗2 (t)| K2u

for t→ + ∞, hence we have lim |(x1 (t), x2 (t)) − (0, x∗2 (t))| = 0. The proof is complete. t→+∞

14

28

1

Acknowledgments

6

This research of YK is partially supported by NSF-DMS (Award Number 1313312 & 1716802); NSFIOS/DMS (Award Number 1558127); DARPA -SBIR 2016.2 SB162-005 Phase II; and The James S. McDonnell Foundation 21st Century Science Initiative in Studying Complex Systems Scholar Award (UHC Scholar Award 220020472). The research of D. Bai is partially supported by NSF of China (11771104) . The research of J. Yu is supported by NSF of China (11631005) and PCSIRT of China (IRT1226).

7

References

2 3 4 5

8 9

10 11

12 13

[1] A.D. Bazykin, Aleksandr Iosifovich Khibnik, Bernd Krauskopf, Nonlinear dynamics of interacting populations(Vol. 11). World Scientific, 1998. [2] H.I. Freedman, R.M. Mathsen, Persistence in predator-prey systems with ratio-dependent predator influence. Bull. Math. Biol., 55, 817-827 (1993). [3] G. F. Gause, N. P. Smaragdova and A. A. Witt, Further studies of interaction between predator and prey. J. Animal Ecol., 5,1-18 (1936).

15

[4] M.P. Hassell, R.M. May, Generalist and specialist natural enemies in insect predator-prey interactions. Br. Ecol. Soc., 55, 923-940 (1986).

16

[5] L. Chen, Mathematical Models and Methods in Ecology. Science Press, Beijing, 1998 (in Chinese).

17

[6] J. M. Smith, J. Model in Ecology. Cambridge University Press, London, 1975.

18 19

[7] A. D. Bazykin, Structural and dynamic stability of model predator-prey systems. Institute of Resource Ecology, 1975.

20

[8] A. D. Bazykin, Problems in Mathematical Genetics U.S.S.R. Adac. Sci., Novosibirsk, 1974.

21

[9] J.M. Cushing, Periodic time-dependent predator-prey system. SIAM J. Appl. Math., 32, 82-95 (1977).

14

22 23

24 25

26 27 28

29 30

31 32

[10] G. P. Samanta, Analysis of nonautonomous two species system in a polluted environment. Math. Slovaca, 62(3), 567-586 (2012). [11] M. Fan, Y. Kuang, Dynamics of a nonautonomous predator-prey system with the Beddington-DeAngelis functional response. J. Math. Anal. Appl., 295, 15-39 (2004). [12] M. Marv´ aa, J.G. Alc´ azara, J.-C. Poggialeb, R.Bravo de la Parraa, A simple geometrical condition for the existence of periodic solutions of planar periodic systems. Applications to some biological models. J. Math. Anal. Appl., 423, 1469-1479 (2015). [13] M. Fan, Q. Wang, X. Zou, Dynamics of a non-autonomous ratio-dependent predator-prey system. Proc. Royal Soc. Edinburgh, 133A, 97-118 (2003). [14] H. Li, Y. Takeuchi, Dynamics of the density dependent and nonautonomous predator-prey system with Beddington-Deangelis functional response. Conti. Dynam. Syst. Seri B., 20(4), 1117-1134 (2015).

34

[15] Y. Kang, J. Fewell, Coevolutionary dynamics of a host-parasite interaction model: obligatory v.s. facultative parasitism. Natural Resource Modeling, 28(4), 398-455 (2015).

35

[16] M. P. Hassell, The Dynamics of Arthropod Predator-Prey Systems. Princeton, U.P., Princeton, 1978.

36

[17] Y. Kang, L. Wedekin, Dynamics of a intraguild predation model with generalist or specialist predator. J. Math. Biol., 67(5), 1227-1259 (2013).

33

37

38 39

[18] I. Hanski, L. Hansson, H. Henttonen, Specialist predators, generalist predators, and the microtine rodent cycle. J. Anim. Ecol., 60, 353-367 (1991).

29

1 2

3 4

5 6

7 8

9 10

11 12

[19] W.E. Snyder, A.R. Ives, Interactions between specialist and generalist natural enemies: parasitoids, predators, and pea aphid biocontrol. Ecol. Soc. Am., 84, 91-107 (2003). [20] S. Madec, J. Casas, G. Barles, C. Suppo, Bistability induced by generalist natural enemies can reverse pest invasions. J. Math. Biol. 75(3), 543-575 (2017). [21] A. Erbach, F. Lutscher, G. Seo, Bistability and limit cycles in generalist predator-prey dynamics. Ecol. Compl., 14, 48-55 (2014). [22] S. Chakraborty, The infuence of generalist predators in spatially extended predator-prey systems. Ecol. Compl., 23 :50-60 (2015). [23] Y. Du, J. Shi, Allee effect and bistability in a spatially heterogeneous predator-prey model. Tran. Amer. Math. Soc., 359(9), 4557-4593 (2007). [24] W.F. Fagan, M.A. Lewis, M.G. Neurbert, P. van den Driessche, Invasion theory and biological control. Ecol. Lett., 5(1), 148-157 (2002).

14

[25] C. Magal, C. Cosner, S. Ruan, J. Casas, Control of invasive hosts by generalist parasitoids. Math. Medi. Biol., 25, 1-20 (2008).

15

[26] V. Hutson, A theorem on average Liapunov functions. Monatshefte für Mathematik, 98, 267-275, 1984.

13

16 17

18 19

20 21

[27] I. Barb˘ alat. Systems d’equations differential d’oscillations nonlinearies. Rev. Roumaine Math. Pure Appl., 4, 267-270 (1959). [28] R.E. Gaines, R.M. Mawhin, Coincidence Degree and Nonlinear Differential Equations. Springer-Verlag, Berlin, 1977. [29] M. Fan, K. Wang, Optimal harvesting policy for single population with periodic coefficients. Math. Biosci., 152, 165-177 (1998).

30