Dynamical effects of biocontrol on the ecosystem: Benefits or harm?

Applied Mathematical Modelling 51 (2017) 361–385

Contents lists available at ScienceDirect

Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm

Dynamical effects of biocontrol on the ecosystem: Benefits or harm? Yun Kang a,∗, Dingyong Bai b, Lorenzo Tapia c, Heather Bateman d a

Sciences and Mathematics Faculty, College of Integrative Sciences and Arts, Arizona State University, Mesa, AZ 85212, USA School of Mathematics and Information Science, Guangzhou University, Guangzhou 510006, China c Applied Mathematics for Life and Social Sciences, School of Human Evolution and Social Change, Arizona State University, AZ 85287, USA d Sciences and Mathematics Faculty, College of Letters and Sciences, Arizona State University, Mesa, AZ 85212, USA b

a r t i c l e

i n f o

Article history: Received 11 June 2016 Revised 24 June 2017 Accepted 3 July 2017 Available online 12 July 2017 Keywords: Predation Native vegetation Saltcedar Invasive species Native consumer

a b s t r a c t Motivated by recent field studies on the effects of biocontrol beetle on invasive plants (Tamarix) on riparian ecosystems, we develop a simple three-species (the native consumer, the biocontrol agent, and their predator) model to explore potential dynamical benefits and harm of introducing a non-indigenous species into varied ecosystems with different habitat structures. Our proposed model assumes that (1) habitat consists of both native plants that are main food resources of the native consumer and invasive plants which are targeted food resources of the biocontrol agents; (2) the habitat structure is stable which has a fixed ratio of invasive plants; (3) biocontrol agents can defoliate invasive plants whose defoliations provide as additional food resources for the native consumer that could potentially increase the native consumer’s abundance; and (4) the predators feed on both the native consumer and the biocontrol agents. Our study shows that introducing biocontrol agents into the ecosystem can generate complicated dynamics that could be beneficial or harmful to the ecosystem depending on the environments: biocontrol agents could be beneficial by promoting its biodiversity such that all species could coexist; on the other hand, biocontrol agents could also be harmful by eliminating the native consumer or the predator. In addition, biocontrol agents could stabilize or destabilize the ecosystem depending on the habit structure and the local ecological environments. These theoretical results provide us potential decision-making tools that allow managers to better predict both safety and efficacy of candidate biological control agents are needed. © 2017 Elsevier Inc. All rights reserved.

1. Introduction In the last half century, invasive species have become a global threat as they dramatically change thriving and diverse ecosystems into uniformly similar landscapes, affecting many native plants and animals in the processes [6]. The increasing abundance of invasive species among ecosystems causes community disassembly and a decline in biodiversity, including species extinction [18,38]. To combat this problem, biological control has been considered as one of important strategies in addition to manual and mechanical removal, prescribed burns, chemical control [11,12,36].

∗

Corresponding author. E-mail address: [email protected] (Y. Kang).

http://dx.doi.org/10.1016/j.apm.2017.07.006 0307-904X/© 2017 Elsevier Inc. All rights reserved.

362

Y. Kang et al. / Applied Mathematical Modelling 51 (2017) 361–385

Biological control aims to control or limit the spread or abundance of a invasive species by natural enemies and typically involves an active human role such as releasing biocontrol agents that are specialists, or predators feeding exclusively on a species of concern [15]. Biocontrol agents are selected and introduced to the ecosystem with expectation that [14]: once introduced, biocontrol agents will feast on the invasive species, thus suppress its reproductive capabilities. Over many years the invasive species will die out or reduce to a manageable size. After this, biocontrol agents should also naturally reduce in size due to lack of resources. The successful biological control is expected to control the effects of invasive species over broad spatial and temporal scales while potentially minimizing risk of adverse consequences [23,32]. However, no biological control agent release is without risk, and agents should be selected by simultaneously considering both risks and potential benefits [23,24,31]. There are many examples of biological control agents safely and effectively controlling invasive species populations [8,19,41]; but there are also many cases where biological control agents have caused harm to the native ecosystem [2,9,23]. The after effects of introducing a biocontrol species to the ecosystem has not be fully understood [12,36], and debate over how best to identify and value potential beneficial and harmful effects of non-indigenous species introduced for biological control purposes is an ongoing conversation in society [34]. As an example, the biological control agent (tamarisk leaf beetle-Dihorabda spp.) is actively being used to defoliate exotic saltcedar (also known as tamarisk) in riparian ecosystems in western USA [4]. Non-native saltcedar is the third most dominant woody genus in riparian areas of the arid western USA, after native willow and cottonwood, which can form dense, monotypic thickets supporting little native vegetation [10] and negatively affect some wildlife species [3,6]. The northern tamarisk beetle, a specialist herbivore, has been used as the biological control agent to control saltcedar [25]. Despite the widespread application of the biocontrol beetles, and the success of saltcedar defoliation, however, only limited research has quantified the consequences on biotic communities, abiotic components, and ecosystem services. Hultine et al. [21] pointed to considerable promise for controlling saltcedar over large areas of invaded habitats; however, there may also be possible negative effects of the biocontrol beetles releases on riparian ecology. Bateman et al. [4] suggest that effects of biocontrol defoliation on ecosystem are likely to be site-specific and depend upon the proportion of native riparian trees established prior to biocontrol introduction and defoliation. Because adverse non-target effects of biocontrol agents have been observed, identifying potential hazards to native species and ecosystems is important for resource managers to balance the benefits and risks of biocontrol [8]. Motivated by this, we aim to develop and study a simple three-species (the native consumer, the biocontrol agent, and their predator) model to explore potential dynamical benefits and harm of introducing a non-indigenous species into varied ecosystems with different habitat structures. Mathematical models have been powerful tools that could provide general statements about the important aspects of a biological control agent should possess in nature [33]. Many researchers have used different modeling approaches to address the impacts of biocontrol agents. Andersen et al. [1] used an agent based model to perform a risk assessment on the introduction of a weed biocontrol agent to an ecosystem with both invasive and native plant species. Shea and Kelly [40] used matrix models to study population growth of weeds and plants with biocontrol agents. Frantzen and Müller– Schärer [13] modeled biocontrol agents in a crop-weed ecosystem, concluding that if density of crops and weeds are low, there is a significant benefit to using the agent. Rao and Rao [37] presented a general microbial model to test three biocontrol methods, insuring the current ecosystem is not disturbed due to the biocontrol. Ghosh et al. [16] use a population model to examine the effects of biocontrol agents might have on the infection rate of malaria in mosquitoes and humans. Gourrley and Lou [17] examined the spread of an invasive insect and its biocontrol agent spatially using a population model. Based on simple mathematical models to make qualitative predictions about transient risk [28], and risk at equilibrium [20] of biological control agents, recent work of [23] used a discrete-time Nicholson–Bailey form consumer-resource model to evaluate conditions under which the population of a non-target host would affect a target host population indirectly via a shared parasitoid, and also to examine what conditions are associated with impacts to non-target hosts upon release of a parasitoid. These published mathematical models indeed provides us useful insights on potential dynamical outcomes of biological control. In this article, motivated by recent field studies of the effects of biocontrol beetles on invasive plants (Tamarix) in riparian ecosystems, we develop a simple native consumer-biocontrol-predator model that incorporates both direct and indirect interactions between the biological control agent and resource populations to assess the potential interplay between both harm and benefit of introducing biological agents to ecosystems. Our proposed model is different from the previous modeling work with an assumption that the established invasive species has a fixed ratio of densities in the ecosystem in a short term. This assumption is based on the fact that the biocontrol agent kills the invasive species slowly but does cause some damage such as defoliation, which is suitable for studying the ecological effects of biocontrol beetles of Tamarix in the riparian ecosystem. The resent studies [3–6] support this assumption as they show that effects of biocontrol defoliation are likely to be site-specific and depend upon the proportion of native riparian trees established prior to biocontrol introduction and defoliation. More specifically, our model assumes that (1) habitat consists of both native plants that are main food resources of the native consumer and invasive plants which are targeted food resources of the biocontrol agents; (2) the current ecosystem has a fixed ratio of invasive plants in a short term; (3) biocontrol agents can defoliate the invasive plants whose defoliations provide as additional food resources for the native consumer that could increase the native consumer’s abundance; and (4) the predators feed on both the native consumer and the biocontrol agents. The proposed model is expected to address the following ecological questions:

Y. Kang et al. / Applied Mathematical Modelling 51 (2017) 361–385

363

1. The habitat structure is measured by the ratio of invasive species in the ecosystem. Since this structure is relatively stable because biocontrol kills the invasive species slowly, how does the ratio of invasive species affect the population dynamics of native consumer and predator? 2. Although biocontrol agents can not eradicate invasive species, they do cause damage such as defoliation that can be beneficial as providing additional food resource for native consumer. How does the degree of such benefits of the biocontrol defoliation affect the ecosystem? 3. What are the synergistic effects of the habitat structure and the biocontrol defoliation that could result in harm or benefits to the ecosystem? The remaining article is organized as follows: in Section 2, we provide detailed derivations of our proposed model that could apply to the scenario of the biological control agent (tamarisk leaf beetle) being actively used to defoliate exotic saltcedar in riparian ecosystems in western USA. In Section 3, we provide completed local and global mathematical analysis of the proposed model to investigate varied scenarios of the effects of biocontrol on the ecosystem. In Section 4, we provide bifurcation diagrams and numerical simulations to explore the potential harm and benefits of introducing biocontrol to ecosystems. In Section 5, we conclude our results and the related biological implications. In addition, we provide potential future work. The detailed proofs of theoretical results of our model are given in the Appendix. 2. Model derivations Motivated by recent field studies of the impact of saltcedar biocontrol on small wildlife, we develop a simple threespecies (the native consumer A, the biocontrol agent B, and their predator L) model to explore potential dynamical benefits and harm of introducing a non-indigenous species B into varied ecosystems with different habitat structures. Let A be population of the native consumer such as arthropods; B be the population of the released biocontrol agents such as beetles; and L be the population of predator of arthropods and the biocontrol beetles which could be several species of insectivorous lizards [6]. We propose the following set of nonlinear equations (1)–(3) to describe the interactions among arthropods, biocontrol beetles and their predator in a ecosystem whose habitat structure is measured by  ∈ (0, 1] which is the fixed ratio of invasive plants in the ecosystem:



A = ra A 1 −



B = rb B 1 −

 L = L

A ( 1 −  )K + db B



−

αa AL 1 + αa hA + αb hB

(1)



B αb BL − K 1 + αa hA + αb hB

ca αa A + cb αb B − dl 1 + αa hA + αb hB

(2)

 (3)

where ri is the intrinsic growth rate of species i; α i is the attacking rate of species L to species i; ci is the energy conversion rate of species i to species L; h is the handling time of species L when it consumes species A or B; dl is the natural mortality of species L; and K is the carrying capacity of the ecosystem without biocontrol B. Model (1)-(3) incorporates two important parameters  and db that model both direct and indirect interactions between the biological control agent and resource populations such that we are able to assess the potential interplay between both harm and benefit of introducing biological agents to the ecosystem:  ∈ (0, 1] is the ratio of invasive plant species that is not edible for species A but it is a targeted food resource for the biocontrol B; and db measures the strength of the positive effects of biocontrol B in the ecosystem, e.g., the larger value of db could potentially increase the abundance of the native consumer A. All parameters, ra , rb ,  , db , dl , α a , α b , h, ca , cb and K are positive. See Fig. 1 for the schematic modeling diagram. The detailed ecological assumptions of the proposed model (1)-(3) can be summarized as follows: 1. K is the carrying capacity of the environment, which measures the combined carrying capacity of the native consumer A and biocontrol B. We assume that the established invasive species has a fixed ratio of densities in the ecosystem in a short term. The habitat has a stable mixed vegetation cover where (1 −  )K is edible for A and the remaining  K is the vegetation covered by the invasive plant which is not edible for A but is a preferred food resource for biocontrol B (i.e., biocontrol B only consumes the invasive plant). 2. In the absence of the biocontrol B and the predator L, the native consumer A follows the traditional logistic growth model



A = ra A 1 −

A ( 1 −  )K



where (1 −  )K is the carrying capacity of A from the native vegetation that is edible for A. 3. In the absence of the predator L, biocontrol B also follows the traditional logistic growth model



B = rb B 1 −

B K



364

Y. Kang et al. / Applied Mathematical Modelling 51 (2017) 361–385

Fig. 1. The schematic modeling concept of Model (1)-(3).

where  K is the carrying capacity of B from the invasive plant species that is a targeted food resource for biocontrol B. Thus, in the absence of the predator L, the ecosystem only consists of native consumer A and biocontrol B whose dynamics follow the two equations below:



A = ra A 1 −



B = rb B 1 −



A , ( 1 −  )K + db B

(4)



B , K

(5)

with an assumption that biocontrol B intensively feeds on the invasive plant  K which leads to the defoliation of the invasive plant. The defoliation provides additional food resources and space for the native consumer species A. As a consequence, the carrying capacity of the native consumer species A becomes (1 −  )K + db B where db measures the positive effect of biocontrol defoliations. 4. Species L preys on both A and biocontrol B with Holling type II functional responses since L spends certain time to search for and handle the consumer A or B or both [7,26,27]. It is possible that predator L may have different handling time for native consumer A and biocontrol B, respectively. However, there is no actual data on the handling time. As a starting point, our model assumes that predator L has the same handling time for both native consumer A and biocontrol B. In the case that there is no biocontrol B in the ecosystem, the dynamics of A and L can be described by the following two equations:



A = ra A 1 −

 L = L

A ( 1 −  )K



−

αa AL , 1 + αa hA

(6)



ca αa A − dl . 1 + αa hA

(7)

If there is no native consumer A in the ecosystem, the dynamics of B and L can be described by the following two equations:



B = rb B 1 − L = L



B αb BL − , K 1 + αb hB

 cαB  b b − dl . 1 + αb hB

(8)

(9)

We could like to address that the parameter  , the ratio of invasive plant species that is not edible for species A but it is a targeted food resource for the biocontrol B, is restricted on the interval (0, 1]. In the case that  = 0, the ecosystem has no invasive plant species, in consequence, there is no biocontrol due to no editable plant, then our propose model (1)-(3) is reduced to the A − L subsystem (6)-(7). In the remaining article, we first provide completed analysis of the three subsystems (i.e., the A − B subsystem (4)-(5); the A − L subsystem (6)-(7), and the A − L subsystem (6)-(7)) and Model (1)-(3). In particular, we compare the dynamics of the A − L subsystem (6)-(7)) and Model (1)-(3) to address: 1. How does the ratio of invasive species  affect the population dynamics of native consumer and predator after introducing the biocontrol B in the ecosystem?

Y. Kang et al. / Applied Mathematical Modelling 51 (2017) 361–385

365

Table 1 l l , Db = α (c d−hd . G.A.S means globally asymptotically stable and L.A.S means The local stability of equilibria for subsystems, where Da = αa (cad−hd l) b b l) locally asymptotically stable. Subsystems

Equilibria

Existence conditions

Stability conditions

A − B subsystem

E00∗ E0B∗ EA0∗ EAB∗ E0∗ 0 EA∗ 0 EA∗ L E∗ 00 E∗ B0 E∗ BL

Always Always Always Always Always Always  < 1 − DKa < 1 Always Always  > DKb > 0

Always source Always saddle Always saddle Always G.A.S. Always saddle Da Da L.A.S. if either  1 >  > 1 − K > 0 or Da < 0; Saddle if  < 1 − K < 1.  L.A.S. if 1 − 2KDa + αa1hK <  < 1 − DKa < 1; Source if  < 1 − 2KDa + αa1hK and Da > 0 Always saddle L.A.S. if either  < DKb or Db < 0; Saddle if  > DKb > 0 L.A.S. if 0 < DKb <  < 2KDb + α 1hK ; Source if  > 2KDb + α 1hK and Db > 0

A − L subsystem

B − L subsystem

b

b

2. How does the degree of positive effects db of the biocontrol defoliation affect the ecosystem after introducing the biocontrol B? 3. What are the synergistic effects of the habitat structure and the biocontrol defoliation that could result in harm or benefits to ecosystem? 3. Mathematical analysis In this section, we perform local and global analysis of our proposed model (1)-(3) and the related subsystems to address our research questions. The state space for Model (1)-(3) is R+ . We first show that the model (1)-(3) is positive invariant 3 and bounded in R+ as the following theorem: 3 Theorem 3.1. (Positively invariant and bounded) Model (1)-(3) is positively invariant and bounded in R+ . In addition, we have 3

lim sup A(t ) ≤ (1 −  )K + db  K, t→∞

lim sup B(t ) ≤  K, t→∞

and

lim sup L(t ) ≤ t→∞

M dl

where M = ca (ra + dl )[(1 −  )K + db  K ] + cb (rb + dl ) K. Biological implications: Theorem 3.1 implies that our proposed model (1)-(3) is biologically well-defined, i.e., the population of species in the ecosystem is nonnegative and bounded. Moreover, the abundance of the native consumer A may be above (1 −  )K due to the positive effects from the biocontrol defoliations. 3.1. Dynamics of subsystems In this subsection, we focus on the dynamics of the three subsystems of Model (1)-(3). First, we introduce the following notations

Da :=

dl

αa (ca − hdl )

,

Db :=

dl

αb (cb − hdl )

where Da is the population of species A at the coexistence equilibrium for the A − L subsystem (6)-(7) and Db is the popd

l ulation of biocontrol B at the coexistence equilibrium for the B − L subsystem (8)-(9). The value Di = α (c −hd is positive if i i l) ci > hdl where i = a, b. In addition, we have the following statements:

1. The A − B subsystem (4)-(5) has three boundary equilibria E00∗ := (0, 0), E0B∗ := (0,  K) and EA0∗ := ((1 −  )K, 0 ), and an interior equilibrium EAB∗ := ((1 −  )K + db  K,  K ). 2. The A − L subsystem (6)   -(7) has two boundary equilibria E0∗ 0 := (0, 0) and EA∗0 := ((1 −  )K, 0 ), and an interior equilibrium EA∗L := Da ,

ra Da αa 1 − (1− )K

(1 + αa hDa ) if (1 −  )K > Da > 0 where Da > 0 requires ca > hdl .

3. B (8)  − L subsystem  -(9) has two boundary equilibria E∗ 00 := (0, 0) and E∗ B0 := ( K, 0), and an interior equilibrium E∗BL := rb Db Db , α 1 −  K (1 + αb hDb ) if  K > Db > 0 where Db > 0 requires cb > hdl . b

Theorem 3.2. (Dynamics of subsystems) The existence and local stability of equilibria for Subsystems (4)-(5), (6)-(7) and (8)-(9) are listed in Table 1. In addition, we have the following statements: 1. The A − B subsystem (4)-(5) is globally stable at EAB∗ . 2. The A − L subsystem (6)-(7) is globally stable at EA∗ 0 if either 1 >  > 1 − DKa > 0 or Da < 0 while it is globally stable at EA∗ L     if 1 − 2KDa + α 1hK <  < 1 − DKa < 1. And the A − L subsystem has a unique stable limit cycle if  < 1 − 2KDa + α 1hK and a a Da > 0.

366

Y. Kang et al. / Applied Mathematical Modelling 51 (2017) 361–385

Table 2 l l , Db = α (c d−hd The boundary equilibria and their local stability for Model (1)-(3),where Da = αa (cad−hd l) b b l) Boundary Equilibria

Existence Conditions

Stability Conditions

E000 , E0B0 , EA00 EAB0

Always Always

Always Saddle  L.A.S. if either (1) Da > 0,  K 1 −

E0BL

 K > Db > 0

L.A.S. if



if either (1) Da > 0,  K 1 − Db α ra 1 − αb r

a b

Da Db





− db > K − Da , or (2) Da < 0,  K 1 −





Da Db



− db < K − Da ; Saddle



Da − db < K − Da , or (2) Da < 0,  K 1 − DDa − db > K − Da . Db b <  K < 2Db + α1h . Source if one of the following conditions holds: (1) ααba rra ≥ 1 and b b

αb r a Db 1 1 αb h ; or (2) αa rb < 1 and 2Db + αb h <  K < 1− αb ra . Saddle if one of the following αa r b αb r a 1 conditions holds: (1) αa r ≥ 1 and Db <  K < 2Db + α h ; or (2) ααba rra < 1 and b b

b

 

 K > 2Db +

Db <  K < min

( 1 −  )K > Da > 0

EA0L

Db α r 1− αb ra

a b

L.A.S. if

Da αa r 1− α r b b a

, 2Db + α1h ; or (3) ααba rra < 1 and  K > max b b

Db α r 1− αb ra

a b

, 2Db + α1h . b

< (1 −  )K < 2Da + αa h . Source if one of the following conditions holds: (1) ααa rrab ≥ 1 b 1

and (1 −  )K > 2Da + α1a h ; or (2) ααa rrab < 1 and 2Da + α1a h < (1 −  )K < b

Da α r 1 − αa r b

. Saddle if one of

b a

the following conditions holds: (1) ααa rrab ≥ 1 and Da < (1 −  )K < 2Da + α1a h ; or (2) ααa rrab < 1 and b b

 Da < (1 −  )K < min

Da αa r 1− α r b



(1 −  )K > max

b a

Da αa r 1− α r b b a

, 2Da + α1a h ; or (3) ααa rrab < 1 and b

3. The B − L subsystem (8)-(9) is globally stable at E∗ B0 if either  <

<

2Db K

+

1 αb hK .



, 2Da + α1a h .

Db K

or Db < 0 while it is globally stable at E∗ BL if 0 <

And the B − L subsystem has a unique stable limit cycle if  >

2Db K

+

1 αb hK

Db K

<

and Db > 0.

Biological implications: the results of Theorem 3.2 imply the follows: 1. In the absence of the predator L, the native species consumer A and the biocontrol B can coexist at the equilibrium EAB∗ = ( (1 −  )K + db  K,  K ) in the A − B subsystem (4)-(5). 2. The large handling time h or the large mortality rate dl of L, i.e., hdl > max {ca , cb }, the predator L goes extinct in the subsystems (6)-(7) and (8)-(9). 3. The dynamics of all the three subsystems does not depend on the positive effects of biocontrol B in the ecosystem db . While the dynamics of the subsystems (6)-(7) and (8)-(9) depend on  which is the ratio of invasive plant species that is not edible for species A but it is

preferable for the biocontrol B. More specifically, we can conclude that (a) If 0 < 1 −

Da K



<  < min 1,

 2Da Db K ,1 − K

Db K

, then the predator L goes extinct in the subsystems (6)-(7) and (8)-(9).





2D

+ α 1hK <  < min K b + α 1hK , 1 − DKa < 1, then both the subsystems (6)-(7) and (8)-(9) a b have global stability at their unique interior equilibrium.   2D (c) If K b + α 1hK <  < 1 − 2KDa + α 1hK with Di > 0, i = a, b, then both the subsystems (6)-(7) and (8)-(9) have an unique

(b) If 0 < max

a

b

stable limit cycle. 3.2. Equilibria and stability of Model (1)-(3) In this subsection, we focus on the equilibria and their stability of Model (1)-(3). First, we could conclude that the system (1)-(3) always has the following three axial equilibria

E000 := (0, 0, 0 ),

E0B0 := (0,  K, 0 ),

EA00 := ((1 −  )K, 0, 0 )

and the planner equilibrium EAB0 := ((1 −  )K + db  K,  K, 0 ). In addition, Model (1)-(3) has additional planner equilibria:



E0BL = 0, and



EA0L = Da ,

Db ,

0,

rb

αb ra

αa



1−



1−





Db D (1 + αb hDb ) , if 0 < b <  K K





Da (1 + αa hDa ) , if ( 1 −  )K

 <1−

Da < 1. K

We have the following results regarding the boundary equilibria of Model (1)-(3). Theorem 3.3. (Boundary equilibria) The existence and local stability of the boundary equilibria

E000 , E0B0 , EA00 , EAB0 , E0BL , EA0L of Model (1)-(3) are summarized in Table 2. Biological implications: the results of Theorem 3.3 provide us a baseline information for studying the persistence of species in Model (1)-(3). According to Theorem 3.3, we have the following implications:

Y. Kang et al. / Applied Mathematical Modelling 51 (2017) 361–385

367

D

1. If 0 < Kb <  < 1 − DKa < 1, then the two planner equilibria E0BL and EA0L exist simultaneously, however, they can not be both locally stable. 2. The existence of the planner EAB0 does not depend on  and db but its stability does. Both the existence and stability of E0BL and EA0L depend on  but not db . 3. Biocontrol effects:from Theorem 3.2, we know that EA∗ L in the A − L subsystem (6)-(7) is globally stable if the following inequalities hold

 2D

1−

a

K

+

1 αa hK



< <1−

Da < 1. K

After introducing the biocontrol B into the A − L subsystem, the corresponding equilibrium of EA∗ L in the full A − B − L α r system is EA0L which is locally stable if αa rab < 1 and

 2D

1−

K

a

+

1

αa hK

b

< <1−

D  a  < 1. K 1 − ααa rrab b

The full A − B − L system has the planner equilibria E0BL and EAB0 which could be locally stable. Based on Theorem 3.3, E0BL is locally stable if the following inequalities hold

0<



Db

K 1 − ααba rra b

 <<

2 Db 1 + . K K αb h

EAB0 is locally stable if one of the following inequalities hold either Da > 0 and

db <

Da + db < 1, Db

> 

K − Da

K 1−

Da Db

− db



or Da < 0 and

Da + db < min{1, db }, Db

< 

K − Da

K 1−

Da Db

− db

.

This implies that introducing the biocontrol B to the A − L subsystem could potentially lead to the extinction of the native consumer A locally, i.e., the planner equilibrium E0BL being locally stable; or lead to the extinction of the predator L, i.e., the planner equilibrium EAB0 being locally stable. 4. Multiple boundary attractors: (a) The two planner equilibria EAB0 and EA0L can coexist and both be locally stable if the following conditions hold:



2D 1 Da K − Da  ,1 − a − + db < 1, max Db K K αa h K 1 − DDa − db

< <1−

b

D  a  <1 K 1 − ααa rrab b

d

l which could be satisfied if Db is a small negative value and h is a small positive value such that Da = α (c −hd < a a l)  

1

αa h

1 αa r 1− α r b

−2

holds.

b a

(b) The two planner equilibria EAB0 and E0BL can also coexist and both be locally stable if one of the following conditions hold: either

db <



Da K − Da D   ,  bαb ra  + db < 1, max Da Db K 1 − D − db K 1 − αa r b

<<

b

or



Da Db 2 Db  <  < min + db < min{1, db }, 0 <  + Db K K 1 − ααba rra b

2 Db 1 + K K αb h

1 1 − Da /K , . αb hK 1 − Da /Db − db

This implies that introducing the biocontrol B to the A − L subsystem could generate complicated dynamics depending on initial conditions due to the existence of multiple boundary attractors. Theorem 3.4. (Interior equilibria) Model (1)-(3) can have at most two interior equilibria when Da Db < 0; and it has no interior equilibrium if Da < 0 and Db < 0. The existence conditions of the unique interior equilibrium E∗ (A∗ , B∗ , L∗ ) are summarized as the following Table 3. Biological implications: according to Theorem 3.4, Model (1)-(3) has no interior equilibrium if Da < 0 and Db < 0 (i.e., max {ca , cb } < dl h), which implies that the native consumer A, the biocontrol B, and the predator L can not coexist globally (i.e., no permanence). Model (1)-(3) can have two interior equilibria when Da < 0 and Db > 0 (i.e., ca < dl h < cb ). For example, when

ra = 1, rb = 1.2,

αa = 0.5, αb = .7, ca = 0.025, cb = 0.08, dl = 0.025, h = 2, K = 100, db = .745,  = 0.45,

368

Y. Kang et al. / Applied Mathematical Modelling 51 (2017) 361–385 Table 3 l l , Db = α (c d−hd . The existence conditions of the interior equilibrium E∗ , where Da = αa (cad−hd l) b b l) Unique Interior Equilibrium ∗

∗

∗

∗

E (A , B , L )

Existence Conditions One of the following conditions holds: 1. Da > 0, Db > 0, αb ra = αa rb ; 2. Da > 0, Db < 0, α b ra ≤ α a rb ; 3. Da > 0, αb ra > αa rb , (1 −  )K < Dαaa rb ; 1− α

b ra

4. Da > 0, Db > 0, αb ra < αa rb ,  K <

Db α r 1− αb ra

;

5. Da < 0, Db > 0, αb ra < αa rb ,  K >

Db α r 1− αb ra

.

a b

a b

Table 4 l l , Db = α (c d−hd The conjunct existence conditions of E0BL and E∗ , and the stability of E0BL and EAB0 , where Da = αa (cad−hd l) b b l) Existence Conditions of E0BL and E∗ Da > 0, Db > 0

Da < 0, Db > 0

αb r a αa rb = 1 and

>

αb r a αa rb > 1 and

 > max

Db K



 Db K

,1 −

Stability of E0BL when E∗ exists

Stability of EAB0 when E0BL exists

Unstable

Saddle when Da > 0 and  >

Db K

>0

Da /K αa r 1− α r b b a

αb r a Db Db /K αa rb < 1 and K <  < 1− αb ra αa r b αb r a Db /K αa rb < 1 and  > 1− αb ra αa r

L.A.S. if

b

Db /K α r 1− αb ra

if  > max

a b

<<

Db /K α ra 1 − αb r

a b

,

2Db K

2Db K

+ α 1hK ; Saddle b



+ α 1hK . b

L.A.S. if Da < 0, and 1−Da /K ; Saddle if 0 < DKb <  < 1−D a /D −d b

Da K 1−Da /K . K−Da /Db −db

Da < 0, db <

>

−

Da Db

b

and

Model (1)-(3) has two interior equilibria (A∗ , B∗ , L∗ ) = (69.12, 35.48, 39.28 ) and (A∗ , B∗ , L∗ ) = (84.65, 43.18, 9.15 ) where (A∗ , B∗ , L∗ ) = (69.12, 35.48, 39.28 ) is locally stable. Under this set of parameters’ values, Model (1)-(3) has an additional locally stable boundary equilibrium EAB0 = (88.53, 45, 0 ). Moreover, based on the results of Theorem 3.3 and 3.4, we have the following implications: D

1. If 0 < Kb <  < 1 − DKa < 1, then the two planner equilibria E0BL and EA0L can coexist simultaneously with only one of them being locally stable. However, if the interior equilibrium E∗ (A∗ , B∗ , L∗ ) also exists, then both E0BL and EA0L are unstable where EAB0 is always a saddle point. 2. The planner equilibria E0BL and EAB0 can be locally stable and coexist with the unique interior equilibrium E∗ under the conditions that EA0L does not exist (e.g., Da < 0) and the following inequalities hold:



Da Db 2 Db  + db < min{1, db }, 0 <  + αb ra <  < min Db K K 1 − αa r b

1 1 − Da /K , . αb hK 1 − Da /Db − db

In this case, we should expect that E∗ is locally unstable and there is no coexistence of all A, B, L three species. The conjunct conditions of the existence for E0BL and the interior equilibrium E∗ , and the stability of E0BL and EAB0 are summarized in Table 4. 3. The planner equilibria EA0L and EAB0 can coexist and being locally stable with the unique interior equilibrium E∗ when α r Da > 0, Db < 0, αa rab < 1 and  is small. However, if 1 − Dαaa rb <  , then EA0L and EAB0 can coexist with the unique interior 1− α r b a

b

equilibrium E∗ but EA0L is always unstable. This implies that introducing the biocontrol B to the A − L subsystem could destabilize the system under certain conditions. 3.3. Persistence and extinction of species in Model (1)-(3) In this subsection, we focus on the persistence and extinction of species in Model (1)-(3). By comparing these results to the persistence and extinction of species in the A − L subsystem (6)-(7), we explore the effects of biocontrol B on the ecosystem. Theorem 3.5. (Persistence of single species) The sufficient conditions for the persistence of single species for Model (1)-(3) are summarized in Table 5. Biological implications: according to Theorem 3.2, we can conclude that the native consumer A is always persistent and the predator L is persistent if the inequalities  < 1 − DKa < 1 hold. After introducing the biocontrol B into the A − L system, α r according to Theorem 3.5, the persistence of the native consumer A heavily depends on the values of Db ,  , and αba ra ; and the b

persistence of the predator L depends on both  and db . More specifically, we have the following two conclusions regarding the effects of biocontrol B on ecosystem: 1. If  < 1 −

Da K

< 1, then the persistence of the predator L in Model (1)-(3) requires either

conditions to be persistent:

Da + db < 1 and Db

< 

K − Da

K 1−

Da Db

− db

.

Da Db

+ db > 1 or the following

Y. Kang et al. / Applied Mathematical Modelling 51 (2017) 361–385

369

Table 5 l l , Db = α (c d−hd . Persistent results for Systems (1)–(3), where Da = αa (cad−hd l) b b l) Persistent species

Sufficient conditions

Species A

One of the following conditions is satisfied: (1) Db < 0; (2)  K < Db ; (3) 0 < Db <  K < 2Db + α1h and ααba rra ≥ 1;





b

b

(4) 0 < Db <  K < min 1−(α rDa )b/(αa r ) , 2Db + α1h and ααba rra < 1. b b b b One of the following conditions is satisfied: (1) Da < 0; (2) (1 −  )K < Da ; αa r b 1 (3) 0 < Da < (1 −  )K < 2Da + αa h and α ra ≥ 1;

Species B





b

(4) 0 < Da < (1 −  )K < min 1−(αa rD)a/(α ra ) , 2Da + α1a h and ααa rrab < 1. b b b One of the following conditions is satisfied: (1) (1 − DDa − db ) K < K − Da and b Da > 0; (2) (1 − DDa − db ) K > K − Da and Da < 0.

Species L

b

These additional conditions indicates that introducing the biocontrol B in to the A − L system would potentially destroy the persistence of the predator L, thus lead to its local extinction when the value of db is very large or the values of both db and  are small. 2. If Da < 0, then we could conclude that the predator L goes extinct in the A − L system. However, if the predator L can still be persistent in Model (1)-(3) if the follow inequalities hold

Da + db < 1 and Db

K − Da

> 

Da Db

K 1−

− db

.

This implies that the small values of db , Db and the large value of  can promote the persistence of the predator L after introducing the biocontrol B into the A − L system where the predator L cannot survive in the absence of the biocontrol L. If Da < 0 and Db < 0, then according to Theorem 3.4, we can conclude that Model (1)-(3) has no interior, thus it can not be permanent. From Theorem 3.5, we can obtain the sufficient conditions of the permanence (i.e., the persistence of all three species) of the full model (1)-(3) in the following two scenarios: Da > 0, Db > 0 versus Da Db < 0. Corollary 3.1. (Permanence conditions for the case Da > 0, Db > 0) In the case that Da > 0, Db > 0, the following statements on the permanence of system (1)-(3) in R3+ hold: 1. Both E0BL and EA0L do not exist: system (1)-(3) is permanent in R3+ if one of the following conditions holds: (1-i) 1 −

Da K



<  < min



Db 1−Da /K K , 1−Da /Db −db

1−Da /K Da K , 1−Da /Db −db

(1-ii) max 1 −

<<

Db K

Da Db ,

and db < 1 −

and db > 1 −

or

Da Db .

2. E0BL exists but EA0L does not exist: system (1)-(3) is permanent in R3+ if one of the following conditions holds:



(2-i) max 1 −



Da Db K , K

(2-ii) max 1 −





Da Db K , K

2Db K

<  < min <<

2Db K

+ α 1hK , b

Db /K α ra 1− αb r a b



α r

and αba ra < 1, or b

α r

+ α 1hK and αba ra ≥ 1. b b

3. EA0L exists but E0BL does not exist: system (1)-(3) is permanent in R3+ if one of the following conditions holds: (3-i) 1 −

 2Da K

(3-ii) 1 − min





+ α 1hK <  < min 1 − a





2Da K

+ α 1hK , a

Da /K αa r 1− α r b b a

Da Db K , K

α r

and αba ra ≤ 1, or b



<  < min 1 −

Da Db K , K

α r

and αba ra > 1. b

4. Both E0BL and EA0L exist: system (1)-(3) is permanent in R3+ if one of the following conditions holds:



(4-i) max

Db K ,1

(4-ii) max

−

Db K ,1

 2Da

−

K

 2Da K





+ α 1hK , 1 − a + α 1hK a



Da /K αa r 1− α r b b a



<  < min 1 −



<  < min 1 −

Da 2Db K , K

Da 2Db K , K

+ α 1hK , b

+ α 1hK b

Db /K α ra 1− αb r a b



α r

and αba ra > 1, or b α r

and αba ra ≤ 1. b

Biological implications: based on Theorem 3.2, we know that the A − L system (6)-(7) has global stability at EA∗ 0 (i.e., the predator L goes extinct) if 1 >  > 1 − DKa > 0 holds; while it is permanent if  < 1 − DKa < 1 holds. However, after introducing the biocontrol B in the A − L system, according to Corollary 3.1, Model (1)-(3) can be permanent even if 1 >  > 1 − DKa > 0. This implies that introducing the biocontrol B could benefit the ecosystem by promoting its diversity and coexistence of all species. On the other hand, if both E0BL and EA0L exist, then the permanence of Model (1)-(3) requires additional conditions comparing to the permanence of the A − L system. In this case, introducing the biocontrol B could harm the ecosystem.

370

Y. Kang et al. / Applied Mathematical Modelling 51 (2017) 361–385

In addition, we can conclude that, if Da > 0, Db > 0, then the permanence of Model (1)-(3) requires intermediate values of  (which is the ratio of invasive plant species in the ecosystem that is not edible for species A but it is preferable for the biocontrol B), and it does not depend on the value of db which describes positive effects of biocontrol B in the ecosystem. This is different from the cases when at least one of Da , Db being negative (see Corollary 3.2 below). If one of Da , Db is negative, then we have the following permanence results for Model (1)-(3): Corollary 3.2. (Permanence conditions for the case Da Db < 0) Assume that db < 1 − nence of system (1)-(3) in

R3+

Da Db ,

the following statements on the perma-

hold:

1. Both E0BL and EA0L do not exist: system (1)-(3) is permanent in R3+ if the following conditions holds:

Da > 0, Db < 0, and 0 < 1 −

Da 1 − Da /K << . K 1 − Da /Db − db

2. E0BL exists but EA0L does not exist: system (1)-(3) is permanent in R3+ if Da < 0, Db > 0 and one of the following conditions holds:



(2-i) max

Db 1−Da /K K , 1−Da /Db −db

(2-ii) max

Db 1−Da /K K , 1−Da /Db −db

2Db K

<<

α r

+ α 1hK and αba ra ≥ 1, or b b



2Db K

<  < min

+ α 1hK , b

Db /K α ra 1− αb r a b



α r

and αba ra < 1. b

3. EA0L exists but E0BL does not exist: system (1)-(3) is permanent in R3+ if Da > 0, Db < 0 and one of the following conditions holds: (3-i) 1 −

 2Da K





+ α 1hK <  < min 1 − a





2Da K

(3-ii) 1 − min

+ α 1hK , a

Da /K αa r 1− α r b b a

1−Da /K Da K , 1−Da /Db −db



<  < min 1 −

α r

and αba ra ≤ 1, or b

1−Da /K Da K , 1−Da /Db −db

α r

and αba ra > 1. b

Biological implications: based on Theorem 3.2, we know that the A − L system (6)-(7) has global stability at EA∗ 0 (i.e., the predator L goes extinct) if Da < 0 holds. However, after introducing the biocontrol B in the A − L system, according to Corollary 3.2, Model (1)-(3) can be permanent even if Da < 0. This implies that introducing the biocontrol B could benefit the ecosystem by promoting its diversity and coexistence of all species. Theorem 3.6. (Extinction of predator L) The predator L goes extinct in Model (1)-(3) if one of the following conditions hold 1. 0 <

D 1− Ka D

1− Da −db

< <1−

Da K

<1

b

2. Da < 0 and Db < 0 (i.e., max {ca , cb } < dl h) where Model (1)-(3) has global stability at EAB0 . Biological implications: based on Theorem 3.2, we know that the A − L system (6)-(7) is permanent (i.e., the predator L is persistent) if  < 1 − DKa < 1 holds. However, after introducing the biocontrol B in the A − L system, according to Theorem 3.6, the predator L goes extinct in Model (1)-(3). This implies that introducing the biocontrol B could potentially harm the dl ecosystem by eliminating the predator L when Db = α (c −hd ) is a small negative number (e.g., when cb is very close but b

b

l

smaller than hdl ). On the other hand, if Da < 0 and Db < 0 (i.e., max {ca , cb } < dl h) holds, then the A − L system (6)-(7) has global stability at EA∗ 0 . However, after introducing the biocontrol B in the A − L system, Model (1)-(3) has global stability at EAB0 . 4. Effects of biocontrol on the ecosystem We first summarize the effects of introducing the biological control B into the A − L system based on theoretical results from Section 3 as follows: •

Potential benefits: the biocontrol B could benefit the ecosystem by promoting its biodiversity and coexistence of all species. dl ∗ 1. If Da = α (c −hd ) < 0 (i.e., ca < hdl ), then the A − L system (6)-(7) has global stability at EA 0 (i.e., the predator L goes a

•

a

l

extinct). After introducing the biocontrol B in the A − L system, according to Theorem 3.3 and Corollary 3.2, Model (1)-(3) can be permanent or the predator L is persistent under certain conditions even if Da < 0. 2. If the inequalities 1 >  > 1 − DKa > 0 hold, the A − L system (6)-(7) has global stability at EA∗ 0 (i.e., the predator L goes extinct). However, after introducing the biocontrol B in the A − L system, according to Corollary 3.1, Model (1)-(3) can be permanent under certain conditions. Potential harm: the biocontrol B could harm the ecosystem by eliminating the predator L or making the native species A extinct locally.

Y. Kang et al. / Applied Mathematical Modelling 51 (2017) 361–385

371

Fig. 2. The two parameter bifurcation diagram on the number of interior equilibria of Model (1)-(3) when ra = 1, rb = 1.2, αa = 5/12, cb = 0.08, dl = 0.025, h = 2, K = 100, ca = 0.025,  ∈ (0, 1] and db ∈ [0, 2]: no interior equilibrium in white regions; one interior equilibrium in black regions; and two interior equilibria in red regions. In Fig. 2(a), αb = .4; and in Fig. 2(b), αb = .7. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article).

1. According to Theorem 3.3, introducing the biocontrol B to the A − L subsystem could complicate the dynamics by generating multiple boundary attractors, which potentially lead to the extinction of the native consumer A locally, i.e., the planner equilibrium E0BL being locally stable; or lead to the extinction of the predator L, i.e., the planner equilibrium EAB0 being locally stable. 2. If the inequalities 0 < 1 − DKa < 1 hold, then the A − L system (6)-(7) is permanent but after introducing the biocontrol B in the A − L system, according to Theorem 3.3 and Corollary 3.2, the predator L goes extinct locally (or globally when D 1− Ka D

1− Da −db

<  also holds) in Model (1)-(3).

b

In the rest of this section, we use numerical simulations and bifurcation diagrams to explore the impacts of biocontrol B on the ecosystem. We focus on the following three scenarios to study specific impacts of  ∈ [0, 1] which is the ratio of invasive plant species that is not edible for species A but it is preferable for the biocontrol B, and db which is the positive effects of biocontrol B, after introducing biocontrol B to the A − L ecosystem. For convenience, we fix the values of the following parameters that could provide typical dynamics:

ra = 1, rb = 1.2,

αa = 5/12, cb = 0.08, dl = 0.025, h = 2, K = 100.

We consider the following typical scenarios according to the sign of Da and the ratio of

ra αb rb αa :

1. If Da < 0 (e.g., ca = 0.025), the A − L system (6)-(7) has global stability at EA∗ 0 . Based on the ratio of following two cases of

ra αb rb αa

> 1 and

ra αb rb αa

ra αb rb αa ,

we have the

< 1.

Da < 0 and ra α b < rb α a . By setting αb = .4, ca = 0.025, the two dimensional bifurcation (see Fig. 2(a)) provides the insights on the effects of  and db on the existence of interior equilibria. It suggests that the small value of  drives the predator L extinct (e.g., when  < 0.1, the white regions in the left of Fig. 2(a)); the large values of  and the small values of db (e.g., the white regions in the right side of Fig. 2(a)) can generate oscillating dynamics that all species coexist without an interior equilibrium (the larger  and the smaller db give the smaller period and amplitude). The system (1)-(3) converges to EAB0 in the black regions of Fig. 2(a). • Da < 0 and ra α > r α a . By setting α = .7, ca = 0.025, the two dimensional bifurcation suggests that large values of b b b db and  can generate bistability between the interior attractor (all species coexist) and the boundary attractor EAB0 (see the blue regions in Figs. 3 and 4 where the system has two interior equilibria). In addition, the large value of  can destabilize the system while the large value of db more likely stabilizes the system (however, if db is too large, the system can drive L extinct). In the white regions of Fig. 2(b) (i.e., no interior equilibrium), Model (1)-(3) converges to the boundary attractor EAB0 ; in the black regions of Fig. 2(b) (i.e., the unique interior equilibrium), Model (1)-(3) converges to the boundary attractor that consists only B and L with oscillating dynamics (i.e., A goes extinct). Biological implications: when Da < 0, the A − L system (6)-(7) can only have the existence of the native consumer A. Our bifurcation diagrams and numerical simulations suggest that introducing the biocontrol B into the A − L system (6)-(7) can be beneficial to the ecosystem by promoting the existence of the predator L (see the blue regions of Figs. 3 and •

372

Y. Kang et al. / Applied Mathematical Modelling 51 (2017) 361–385

Fig. 3. The bifurcation diagram on the number of interior equilibria of Model (1)-(3) when ra = 1, rb = 1.2, αa = 5/12, αb = .7, cb = 0.08, dl = 0.025, h = 2, K = 100, ca = 0.025, db = 0.745, and  ∈ (0, 1]: the interior equilibrium in green is saddle; in blue is locally stable; in red is source. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article).

Fig. 4. The bifurcation diagram of db of Model (1)-(3) when ra = 1, rb = 1.2, αa = 5/12, αb = .7, cb = 0.08, dl = 0.025, h = 2, K = 100, ca = 0.025,  = 0.45, and db ∈ [0, 2]: the interior equilibrium in green is saddle; in blue is locally stable; in red is source. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.).

4 where all species coexist with stable dynamics; and the white regions in the right of Fig. 2(a) where all species coexist with oscillating dynamics). The downside is that the native consumer A could go extinct and is replaced by the coexistence of B and L (e.g., see the black regions of Fig. 2(b)).   2. If Da > 0 (e.g., when ca = 0.055, we have 1 − DKa = 0.88 and 1 − 2KDa + α 1hK = 0.748), the A − L system (6)-(7) has global a

stability at EA∗ 0 when  > 1 − DKa > 0; it has global stability at the unique interior equilibrium EA∗ L under the conditions   that 1 − 2KDa + α 1hK <  < 1 − DKa < 1; it has a unique stable limit cycle around EA∗ L under the conditions that  <

 2Da



a

r α

r α

r α

+ α 1hK < 1. Based on the ratio of ra αab , we have the following two cases of ra αab > 1 and ra αab < 1. a b b b • Da > 0 and ra α < r α a by setting α = .4, ca = 0.055 (see Fig. 5(a)). The system (1)-(3) has oscillating dynamics with b b b all species coexistence even in the white regions of Fig. 5(a) (i.e.,  < 0.1). Moreover, the larger value of  give the smaller period and amplitude of the oscillations. • Da > 0 and ra α > r α a by setting α = .7, ca = 0.055 (see Fig. 5(b)). In some white regions of Fig. 5(b) when  is small b b b (e.g.,  < 0.1), the system (1)-(3) has oscillating dynamics with all species coexistence (see the interesting oscillation patterns shown in Fig. 6(a)). In black regions of Fig. 5(b) (i.e.,  > 0.6), the system (1)-(3) converges to oscillating dynamics with only species B and L. Biological implications: for the chosen parameter values, we should expect that the A − L system (6)-(7) can only have the existence of the native consumer A when  > 0.88 and the coexistence of A and the predator L when  < 0.88. Our 1−

K

Y. Kang et al. / Applied Mathematical Modelling 51 (2017) 361–385

373

Fig. 5. The two parameter bifurcation diagram on the number of interior equilibria of Model (1)-(3) when ra = 1, rb = 1.2, αa = 5/12, cb = 0.08, dl = 0.025, h = 2, K = 100, ca = 0.055,  ∈ (0, 1] and db ∈ [0, 2]: no interior equilibrium in white regions and one interior equilibrium in black regions. In Fig. 5(a), αb = .4; and in Fig. 5(b), αb = .7.

Fig. 6. The population dynamics of Model (1)-(3) when ra = 1, rb = 1.2, αa = 5/12, cb = 0.08, dl = 0.025, h = 2, K = 100. In Fig. 6(a), αb = .7, ca = 0.055, db = 0.15,  = 0.05; and in Fig. 6(b), αb = .4, ca = 0.025, db = 0.0745,  = 0.75.

bifurcation diagrams and numerical simulations suggest that introducing the biocontrol B into the A − L system (6)-(7) can be beneficial to the ecosystem by promoting the existence of the predator L (see the white regions Fig. 5(a) when  > 0.88 where all species coexist with oscillating dynamics). However, introducing the biocontrol B into the A − L system (6)-(7) can be potentially harmful to the ecosystem by replacing the native species A with the biocontrol B (see the black regions of Fig. 5(b) when 0.6 <  < 0.88 where only the biocontrol B and the predator L coexist, however, the original A − L system (6)-(7) has the coexistence of A and the predator L). According to bifurcation diagrams and numerical simulations, we can conclude that Model (1)-(3) has rich and complicated dynamics. Introducing the biocontrol B into the A − L system (6)-(7) could be beneficial or harmful to the ecosystem depending on the specific ecological environment defined by the parameter values. Moreover, we have the following observations: 1. If Da > 0, it seems that the value of db has no effects on the existence of interior equilibrium. The larger value of  is r α more likely to generate interior equilibrium when ra αab < 1 while the smaller value of  is more likely to generate interior equilibrium when

ra αb rb αa

b

> 1.

2. If Da < 0, then the large values of db are more likely to generate multiple interior equilibrium. More specifically, when Da < 0, the large values of db and  can generate two interior equilibria, and thus create a bistability between the interior coexistence attractor and the boundary attractor where the predator L goes extinct.

374

Y. Kang et al. / Applied Mathematical Modelling 51 (2017) 361–385

3. The large values of  could destabilize the dynamics of Model (1)-(3); while the large values of db could stabilize the dynamics. 4. Model (1)-(3) is not permanent when it has no interior equilibrium. However, it could have oscillating dynamics with all species coexistence for most initial conditions (see Fig. 6). In this case, the smaller values of  generate oscillating dynamics with smaller periods. 5. Discussion Nonnative invasive plants are considered a major conservation threat to native ecosystems as they can impact both composition and habitat physiognomy affecting ecosystem function [6]. For example, throughout riparian areas of the southwestern United States, non-native saltcedar (also known as tamarisk) can form dense, monotypic stands and have detrimental effects on native plants and habitat quality [39]. To control saltcedar, the biocontrol agent, northern tamarisk beetle, has been used to defoliate nonnative saltcedar in USA western riparian systems since 2001 [5]. Saltcedar biocontrol, like other control methods, has the potential to affect non-target species, thus impact biotic communities and climatic conditions in affected riparian areas [4]. Since adverse non-target effects of biocontrol agents have been observed, identifying potential hazards to native species and ecosystems is important for resource managers to balance the benefits and risks of biocontrol [5,8]. Motivated by this, we propose and study a simple native consumer-biocontrol-predator model that incorporates both direct and indirect interactions between the biological control agent and resource populations to investigate how the synergistic effects of the habitat structure (i.e., the proportion of invasive plant species established prior to biocontrol introduction) and the biocontrol caused defoliation that could result in harm or benefits to ecosystem after introducing biocontrol. Our model assumes that the proportion of invasive plant species established in the ecosystem remains the same as prior to biocontrol introduction. This bases on the fact that biocontrol agents do not eradicate the invasive species but does cause the plant to defoliate and can slowly cause plant mortality. In addition, the biocontrol defoliation provides an additional resource for native consumer thus could potentially increase the abundance of native consumer. The proportion of the invasive species and the positive effect of biocontrol defoliation on native consumer’s abundance are two important ecological factors that allow us to access the potential interplay between both harm and benefits of introducing biological agents to the ecosystem. Our theoretical study shows that introducing biocontrol of invasive plant species into ecosystems could generate complicated dynamics that may possess both risks and benefits depending ecological environments. Our work shows that (1) the larger proportion of invasive plant species established in the ecosystem is more likely to destabilize the system by generating oscillating dynamics (e.g., see Fig. 3); (2) the stronger positive effect of biocontrol defoliation on native consumer’s abundance could stabilize the system (e.g., see Fig. 4); and (3) if the predator can survive before introducing biocontrol, the ecosystem with the smaller proportion of invasive plant species is more likely to have coexistence of all species (see discussion on the case of Da > 0 and Fig. 5). The synergistic interactions of the habitat structure (i.e., the proportion of the invasive species) and biocontrol defoliation (i.e., the positive effect on native consumer’s abundance) could create complicated dynamical outcomes including the following two types of multiple attractors. 1. Two boundary attractors: according to Theorem 3.3, the full A − B − L system could have two types of multiple boundary attractors after introducing the biocontrol B into the A − L system: • The coexistence of the native consumer and the biocontrol (i.e., EAB0 ) versus the coexistence of the biocontrol and the predator (i.e. E0BL ): this scenario is more likely to occur when the ecosystem has the smaller proportion of invasive plant species and not too strong positive effect of biocontrol defoliation on native consumer’s abundance. For 5 example, when ra = 1, rb = 1.2,  = .125, αa = 12 , αb = .25, db = .25, ca = 0.025, cb = 0.08, K = 50, h = 1.75, dl = 0.025, the full system has EAB0 = (45.31, 6.25, 0 ) and E0BL = (0, 2.76, 5.92 ) being locally stable. • The coexistence of the native consumer and the biocontrol (i.e., EAB0 ) versus the coexistence of the native consumer and the predator (i.e. EA0L ): this scenario is more likely to occur when the ecosystem has the larger proportion of invasive plant species and not too strong positive effect of biocontrol defoliation on native consumer’s 5 abundance. For example, when ra = 1, rb = 1.2,  = .55, αa = 12 , αb = 1.55, db = .25, ca = 0.055, cb = 0.025, K = 50, h = 2., dl = 0.025, the full system has EAB0 = (29.375, 27.5, 0 ) and EA0L = (12, 0, 12.32 ) being locally stable. 2. Bistability between the boundary attractor where A goes extinct and the interior coexistence attractor: this scenario is more likely to occur when (a) the predator cannot survive before introducing biocontrol (i.e., Da < 0); (b) the relative growth a of native consumer is greater than its relative attacking rate by the predator (i.e, rra > α α ); and (c) the ecosystem has b

b

the intermediate proportion of invasive plant species and the stronger (but not too strong) positive effect of biocontrol defoliation on native consumer’s abundance (see red regions of Fig. 2(b)). The full A − B − L system could have coexistence of all species through oscillating dynamics even without an interior equilibrium when (a) the predator cannot survive before introducing biocontrol; (2) the relative growth of native consumer is a smaller than its relative attacking rate by the predator (i.e, rra < α α ); and (c) the ecosystem has the large proportion of inb

b

vasive plant species and the weak positive effect of biocontrol defoliation on native consumer’s abundance (see Fig. 6(b)). In the case that the predator can survive before introducing biocontrol and the relative growth of native consumer is larger than its relative attacking rate by the predator, then the full A − B − L system could have coexistence of all species through oscillating dynamics even without an interior equilibrium when ecosystem has a very small proportion of invasive plant species (see Fig. 6(a)).

Y. Kang et al. / Applied Mathematical Modelling 51 (2017) 361–385

375

Fig. 7. Comparison between data (Fig. 7(a)) and simulations (Fig. 7(b)) of Model (1)-(3) regarding population dynamics of A; B; and L before versus after introducing biocontrol agents-B. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article).

Our study suggests that the biocontrol B could harm the ecosystem by eliminating the predator L or making the native species A extinct locally. For example, introducing the biocontrol B into the A − L system (6)-(7) could drive the extinction of the native species A which is replaced by the biocontrol B (see the black regions of Fig. 5(b) when 0.6 <  < 0.88 where only the biocontrol B and the predator L coexist, however, the original A − L system (6)-(7) has the coexistence of A and the predator L). According to Theorem 3.3 and Corollary 3.2, introducing the biocontrol B in the A − L system could result in the extinction of predator L when the proportion of invasive plant species is large even though the original A − L system has the coexistence of both native consumer A and predator L globally. On the other hand, the biocontrol B could also benefit the ecosystem by promoting its biodiversity such that all species could coexist. According to Theorem 3.3 and Corollary 3.2, introducing the biocontrol B in the A − L system could lead to the coexistence of all species even when the predator cannot survive before introducing biocontrol (i.e., Da < 0). The numerical examples of the cases when Da < 0 also support this (see Fig. 6(b) and the red regions of Fig. 2(b)). Future study direction: our current research is a theoretical study motivated by the biocontrol beetle of invasive plant salcedar, which is expected to apply for a general setting of biocontrol issues in ecological communities. Our results indeed inform us the rich dynamics generated by the biocontrol agents, and the potential interplay between both risks and benefits of introducing biological agents to ecosystems. As an illustration, we provide an example of data collected by Bateman and her lab in [30] regarding the effects of introducing biocontrol agents into ecosystem versus our simulations of the similar situation to address potential applications of our work in biocontrol. As shown in Fig. 7: Fig. 7(a) shows time series of data collected on capture rates of arthropods (orange color) and lizards (cyan color) before (from 2009 to 2010) and after (from 2011 to 2013) introducing biocontrol agents. The data suggests that arthropod taxa (that lizards prefer as food) increased in their study sites following biocontrol defoliation. This observation has been confirmed quantitatively in our simulations shown

376

Y. Kang et al. / Applied Mathematical Modelling 51 (2017) 361–385

in Fig. 7(b) of Model (1)-(3) by taking ra = 1, rb = 1.2, αa = 5/12, αb = 0.15, ca = 0.075, cb = 0.08, db = 0.98, dl = 0.035, h = 2, K = 50, and  = 0.25. The case of before introducing biocontrol agents B is given by A(0 ) = 50; B(0 ) = 0; L(0 ) = 60; while the case of after introducing biocontrol agents B is given by A(0 ) = 50; B(0 ) = 40; L(0 ) = 60. Fig. 7 only provides a glimpse of the potential application of our proposed model (1)-(3). Much more work is needed in the connection of our model and results to the field data. The important extension of our current study is to apply the more detailed data collected in the field [3–6] to potentially modify and/or validate our model, and perform parameter estimations of life history parameters that are not easy to obtain in the filed work. This is our ongoing project. There are many potential paths to improve our model to make it more realistic and applicable for real ecosystems. For examples, 1. Dynamics of invasive species: the current model assumes that invasive plant species established in the ecosystem has a fixed proportion due to biocontrol killing invasive species slowly. However, in a long run, the ratio of invasive species should be expected to decrease due to the biocontrol defoliations. It would be more realistic to incorporate multiple time scales in the model to include the dynamics of the ratio of invasive plant species in the ecosystem occurring in the larger time scale. 2. Spatial effects: our current work does not include space which is an important contribute to the dynamics of native species, invasive species, their consumers and the related biocontrol agents. An individual plant species growing in different or changing environments or during different time periods can vary greatly in its life history, which in return, could affect population dynamics of its consumer or biocontrol agents. One of the most recent work in the biocontrol modeling that includes space is by Parshad et al. [35]. They propose a spatiotemporal population model for a threespecies food chain, where an invasive species functions as the top predator, that depredates on a middle predator, which in turn depredates on a prey species. Their proposed model provides a predictive tool, that can be tuned depending on whatever a conservation biologist or manager conceives as a high level, and to then model the effect of our proposed controls in driving down the invasive population, before that level is reached. To extend our current work by including space, it would be interesting to model the vegetation according to spatial location, i.e.,  (x), the ratio of the invasive species is a function of its spatial location x, to explore how landscape driven spatial variations affect biological control processes and outcomes. The large population of biocontrol B feeding on the invasive plant causes defoliation, thus decreases the shady area needed for the predator which may create less comfortable environment for L. Thus, it would be interesting to include the negative effects of biocontrol defoliations on the predator in addition to its positive effects on the abundance of the native consumer. Also, some studies have found that ants can act as predators on the biocontrol and cause it to not become fully established [5]. Moreover, seasonal effects are important that should be considered for the future modeling study. Acknowledgments This research of YK is partially supported by NSF-DMS (Award Number 1313312); NSF- IOS/DMS (Award Number 1558127) and The James S. McDonnell Foundation 21st Century Science Initiative in Studying Complex Systems Scholar Award (UHC Scholar Award 220020472). D-Y. B.’s research is supported by the Natural Science Fund of China (11371107). Appendix: Proofs Proof of Theorem 3.1 Proof. For A ≥ 0, B ≥ 0 and L ≥ 0, we have

dA dB dL |A=0 = 0, |B=0 = 0, |L=0 = 0. dt dt dt which implies that A = 0, B = 0 and L = 0 are invariant manifolds, respectively. Due to the continuity of the system, we can conclude that Model (1)-(3) is positively invariant in R+ . 3 By the Eq. (2), we have B (t ) ≤ rb B(1 − BK ), which implies that lim sup B(t ) ≤  K. This indicates that for any δ > 0, there t→∞

exists T1 > 0 large enough, such that B(t ) ≤  K + δ for all t > T1 . Then, by the Eq. (1),

A (t ) ≤ ra A(1 −



≤ ra A 1 − =

A

( 1 −  )K + db B

)

A ( 1 −  )K + db (  K + δ ) ra A

( 1 −  )K + db (  K + δ )



[(1 −  )K + db ( K + δ ) − A],

t > T1 .

Y. Kang et al. / Applied Mathematical Modelling 51 (2017) 361–385

377

It follows lim sup A(t ) ≤ (1 −  )K + db  K. Hence, there exists T2 > T1 , such that A(t ) ≤ (1 −  )K + db  K + δ for all t ≥ T2 . Now, t→∞

Let N = ca A + cb B + L. Then

dN dA dB dL = ca + cb + dt dt dt dt ≤ ca ra A + cb rb B − dl L = ca ( ra + dl )A + cb ( rb + dl )B − dl N ≤ c a ( r a + d l )[ ( 1 −  ) K + d b  K + δ ] + c b ( r b + d l )[  K + δ ] − d l N ≤ M + δ [ca (ra + dl ) + cb (rb + dl )] − dl N,

t ≥ T2 .

Here, M = ca (ra + dl )((1 −  )K + db  K ) + cb (rb + dl ) K. Thus, we can conclude that lim sup N (t ) ≤ (3) is bounded in

R+ . 3

t→∞



M . dl

Therefore, Model (1)-

Proof of Theorem 3.2 Proof. 1. Consider A − B Subsystem (4)-(5). It is easy to see that (4)-(5) has three boundary equilibria E00∗ := (0, 0), E0B∗ := (0,  K) and EA0∗ := ((1 −  )K, 0 ), and an interior equilibrium EAB∗ := ((1 −  )K + db  K,  K ). The stability of the these equilibria is determined by the eigenvalues λi (i = 1, 2 ) of the following Jacobian matrix J associated to Subsystem (4)-(5), evaluated at these equilibria:

⎛

ra −

⎜

J |E = ⎝

2ra A ( 1 −  )K + db B 0

⎞

ra db A2 ((1 −  )K + db B )2 ⎟ ⎠ 2r B rb − b K

For the boundary equilibrium E00∗ = (0, 0 ), since



J|E00∗ =

ra 0

0 rb



which gives λ1 = ra , λ2 = rb . Thus, E00∗ = (0, 0 ) is always a source. For the boundary equilibrium E0B∗ = (0,  K ), since



ra 0

J |E0B∗ =

0 −rb



which gives λ1 = ra , λ2 = −rb . Thus, E0B∗ = (0,  K ) is always a saddle point. For the boundary equilibrium EA0∗ = ((1 −  )K, 0 ), since



−ra 0

J |EA0∗ =

ra db rb



which gives λ1 = −ra , λ2 = rb . Thus, EA00 = ((1 −  )K, 0, 0 ) also is always a saddle point. For the interior equilibrium EAB∗ = ((1 −  )K +  K db ,  K ), since



−ra 0

J|EAB∗ =

ra db −rb



which gives λ1 = −ra , λ2 = −rb . Thus, EAB∗ is always locally asymptotically stable. 2. Consider A − L Subsystem (6)-(7). It is easy  to see that  (6)-(7) has two boundary equilibria E0∗ 0 := (0, 0) and EA∗0 := ((1 −  )K, 0 ), and an interior equilibrium EA∗L := Da , αraa 1 − (1−Da )K (1 + αa hDa ) if (1 −  )K > Da > 0. The stability of these equilibria is determined by the eigenvalues λi (i = 1, 2 ) of the following Jacobian matrix J associated to Subsystem (6)-(7), evaluated at these equilibria:

⎛

J |E

αa L ra − (12−raA)K − (1 + αa hA )2 ⎜ =⎝ αa ca L (1 + αa hA )2

⎞ αa A 1 + αa hA ⎟ ⎠ ca αa A − dl 1 + αa hA −

For the boundary equilibrium E0∗0 = (0, 0 ), since



J|E0∗0 =

ra 0

0 −dl



which gives λ1 = ra , λ2 = −dl . Thus, E0∗0 = (0, 0 ) is always a saddle point.

378

Y. Kang et al. / Applied Mathematical Modelling 51 (2017) 361–385

For the boundary equilibrium EA∗0 = ((1 −  )K, 0 ), since

⎛

−ra

⎜

J|EA∗0 = ⎝

0

⎞ (1 −  )K αa 1 + (1 −  )K αa h ⎟ ⎠ ca αa (1 −  )K − dl 1 + αa h(1 −  )K −

which gives λ1 = −ra , λ2 =

λ2 =

ca αa (1− )K 1+αa h (1− )K

− dl . Notice

ca αa (1 −  )K − dl < 0 ⇔ (1 −  )K αa (ca − hdl ) < dl . 1 + αa h(1 −  )K d

l Thus, if either ca − hdl < 0, i.e., Da = α (c −hd < 0, or (1 −  )K < Da , then EA∗ 0 is locally asymptotically stable; if (1 −  )K > a a l) ∗ Da > 0, then EA 0 is a saddle point.  

For the interior equilibrium EA∗L = Da ,

ra Da αa 1 − (1− )K

istence of EA∗ L . Since

⎛

ra −

⎜

J|EA∗L = ⎝

2 r a Da ra (1 − Da /((1 −  )K )) − ( 1 −  )K 1 + αa hDa ca ra (1 − Da /((1 −  )K )) 1 + αa hDa

−

(1 + αa hDa ) . Let (1 −  )K > Da > 0, which guarantees the ex-

⎞ αa Da 1 + αa hDa ⎟ ⎠ 0

the eigenvalues λ1 and λ2 are the roots of the following equation

  2 r a Da ra (1 − Da /((1 −  )K )) ca αa ra Da (1 − Da /((1 −  )K )) λ2 − λ ra − − + = 0. ( 1 −  )K 1 + αa hDa (1 + αa hDa )2

Since (1 −  )K > Da > 0, the real parts of λ1 and λ2 have the same signs. It is easy to check that

ra (1 − Da /((1 −  )K )) <0 1 + αa hDa 1 ⇔ ( 1 −  )K < 2 Da + . αa h

Reλ1 < 0, Reλ2 < 0 ⇔ ra −

2 r a Da

( 1 −  )K

−

Therefore, if 0 < Da < (1 −  )K < 2Da + α1h , then EA∗ L is locally asymptotically stable; if (1 −  )K > 2Da + α1h and Da > 0, a a then EA∗ L is a source. 3. Similar to the arguments above, one can prove the local stability of equilibria for B − L Subsystem (8)-(9). To prove the results on the global dynamics, we first introduce a result about the existence of limit cycles for Subsystem (6)-(7) and (8)-(9), which can be obtained directly from the Theorem 1.1 in [29]. Lemma A.1. (1) Let Da > 0. Subsystem (6)-(7) has one limit cycle if and only if (1 −  )K > 2Da + α1h . a (2) Let Db > 0. Subsystem (8)-(9) has one limit cycle if and only if  K > 2Db + α1h . b

Now we are presenting the proof of the global dynamics as follows: 1. From the results on the local dynamics, we know that E00∗ is a source, E0B∗ and EA0∗ both are saddle points, and EAB∗ is locally asymptotically stable for the A − B subsystem (4)-(5). In order to prove the global stability of EAB∗ , we only need to show that the subsystem (4)-(5) has no limit cycle in the interior of R2+ . Define a function φ (A, B ) = A−2 B−2 . Then it is easy to check that

∂ ( φ A ) ∂ ( φ B ) + = −A−2 B−2 (ra + rb ) < 0 ∂A ∂B holds in the interior of R2+ . It follows by Dulac’s criterion (Guckenheimer and Holmes [22]) that the subsystem (4)-(5) has no limit cycle in the interior R+ . (Note: the proof of the global stability of EAB∗ also can be completed by showing that for any 2 initial condition taken in R+ with A(0)B(0) > 0, limt→∞ B(t ) =  K and limt→∞ A(t ) = (1 −  )K + db  K.) 2 2. If either (1 −  )K < Da or Da < 0, then the A − L subsystem (6)-(7) has only two boundary equilibria E0∗ 0 , which is saddle, and EA∗ 0 , which is local asymptotically stable. Thus, from the Poincare´ –Bendixson theorem (Guckenheimer and Holmes [22]), we can conclude that the subsystem (6)-(7) is globally stable at EA∗ 0 . If (1 −  )K > Da > 0, Lemma A.1 indicates that (6)-(7) is globally stable at EA∗ L . 3. Similar to the arguments above for the A − L subsystem (6)-(7), one can get the global stability of B − L subsystem (8)-(9).  Proof of Theorem 3.3 Proof. The existence of the three axial equilibria E000 , E0B0 , EA00 and the three planner equilibria EAB0 , E0BL , EA0L is easy to check. The stability of these boundary equilibria is determined by the eigenvalues λi (i = 1, 2, 3 ) of the following Jacobian

Y. Kang et al. / Applied Mathematical Modelling 51 (2017) 361–385

379

matrix J associated to the system (1)-(3), evaluated at these equilibria:

⎛ ⎜ ⎜ ⎜ J |E = ⎜ ⎜ ⎝

αa L(1 + αb hB ) 2ra A − (1 −  )K + db B (1 + αa hA + αb hB )2 αa αb hBL (1 + αa hA + αb hB )2 αa [ca + (ca − cb )αb hB]L (1 + αa hA + αb hB )2

ra −

ra db A2 αa αb hAL + ((1 −  )K + db B )2 (1 + αa hA + αb hB )2 2r B αb L(1 + αa hA ) rb − b − K (1 + αa hA + αb hB )2 αb [cb + (cb − ca )αa hA]L (1 + αa hA + αb hB )2

−

αa A

1 + αa hA + αb hB ⎟ ⎟ αb B ⎟ − ⎟ 1 + αa hA + αb hB ⎟ ⎠ c a αa A + c b αb B − dl 1 + αa hA + αb hB

1. For the boundary equilibrium E000 = (0, 0, 0 ), since



J |E0 =

ra 0 0

0 rb 0

0 0 −dl



which gives λ1 = ra , λ2 = rb , λ3 = −dl . Thus, E000 = (0, 0, 0 ) is always a saddle point. 2. For the boundary equilibrium E0B0 = (0,  K, 0 ), since

⎛

ra

⎞

0

⎜ ⎜ 0 −rb ⎝

J |E0B0 = ⎜

0

0

0 ⎟  K αb ⎟ − 1 +  K αb h ⎟ ⎠ cb αb  K − dl 1 + αb h K

which gives λ1 = ra , λ2 = −rb , λ3 =

cb αb  K 1+αb h K

− dl . Thus, E0B0 = (0,  K, 0 ) is always a saddle point.

3. For the boundary equilibrium EA00 = ((1 −  )K, 0, 0 ), since

⎛

−ra

⎜ ⎜ ⎝

ra db

J|EA00 = ⎜ 0

rb

0

0

⎞ (1 −  )K αa 1 + (1 −  )K αa h ⎟ ⎟ 0 ⎟ ⎠ ca αa (1 −  )K − dl 1 + αa h(1 −  )K −

ca αa (1− )K which gives λ1 = −ra , λ2 = rb , λ3 = 1+ αa h(1− )K − dl . Thus, EA00 = ((1 −  )K, 0, 0 ) also is always a saddle point. 4. For the planner equilibrium EAB0 = ((1 −  )K +  Kdb ,  K, 0 ), since

⎛

−ra

J|EAB0

⎜ ⎜ ⎜ =⎜ 0 ⎜ ⎝ 0

ra db −rb 0

⎞ αa ((1 −  )K + db  K ) 1 + αa h((1 −  )K + db  K ) + αb h K ⎟ ⎟ rb  K ⎟ − ⎟ 1 + αa h((1 −  )K + db  K ) + αb h K ⎟ ⎠ ca αa ((1 −  )K + db  K ) + cb αb  K − dl 1 + αa h((1 −  )K + db  K ) + αb h K −

which gives λ1 = −ra , λ2 = −rb , λ3 =

ca αa ((1− )K+db  K )+cb αb  K 1+αa h ((1− )K+db  K )+αb h K

ca αa ((1 −  )K + db  K ) + cb αb  K 1 + αa h((1 −  )K + db  K ) + αb h K ⎧ D ⎪ ⎨(1 −  )K + db  K < Da − a  K if Db ⇔ Da ⎪ ⎩(1 −  )K + db  K > Da −  K if Db ⎧  ⎨ K 1 − Da − db > K − Da if Da Db  ⇔ ⎩ K 1 − Da − db < K − Da if Da Db dl >

Da > 0 Da < 0 >0 <0

and a saddle point if

ca αa ((1 −  )K + db  K ) + cb αb  K 1 + αa h((1 −  )K + db  K ) + αb h K ⎧  D ⎪ ⎨ K 1 − a − db < K − Da if Da > 0 Db ⇔  D ⎪ K 1 − a − d > K − D if D < 0. ⎩ a a b Db dl <

⎞

− dl . Thus, EAB0 is locally asymptotically stable if

380

Y. Kang et al. / Applied Mathematical Modelling 51 (2017) 361–385



r



D





5. Consider the stability of the planner equilibrium E0BL = 0, Db , αb 1 −  Kb (1 + αb hDb ) . Let  > b antees the existence of E0BL . Since

⎛

ra −

J|E0BL

⎜ ⎜ ⎜ =⎜ ⎜ ⎝

⎞ D αa rb  1− b 0 0 αb K ⎟ ⎟ 2 r b Db rb (1 − Db /( K )) ∂ F2 αb Db ⎟ |E0BL rb − − − ⎟ ∂A K 1 + αb hDb 1 + αb hDb ⎟ ⎠ ∂ F3 cb rb (1 − Db /( K )) |E0BL 0 ∂A 1 + αb hDb

we get the characteristic equation as follows

   D αr  λ − ra − a b 1 − b · αb K     2r D r (1 − Db /( K )) c r (1 − Db /( K ))αb Db + b b = 0. λ2 − λ rb − b b − b K 1 + αb hDb (1 + αb hDb )2 

α r

D



Thus, the eigenvalue λ1 = ra − αa b 1 −  Kb , and λ2 and λ3 are the roots of the following equation b

  2 r b Db rb (1 − Db /( K )) c r (1 − Db /( K ))αb Db λ − λ rb − − + b b = 0. K 1 + αb hDb (1 + αb hDb )2 2

In order to determine the stability, we consider the signs of λi (i = 1, 2, 3 ). First,

λ1 = ra − ⇔

D αa rb  1− b <0 αb K

αb ra Db < 1 and  K > ( > Db ) . αa rb 1 − ααba rra b

Since  K > Db > 0, the real parts of λ2 and λ3 have the same signs.

2 r b Db r (1 − Db /( K )) − b <0 K 1 + αb hDb 1 ⇔  K < 2 Db + . αb h

Reλ2 < 0, Reλ3 < 0 ⇔ rb −

Therefore, the planner equilibrium E0BL is locally asymptotically stable if λi < 0(i = 1, 2, 3 ). That is,

Db 1 <  K < 2 Db + . αb h 1 − ααba rra b

It is a source if λi > 0(i = 1, 2, 3 ). That is,

⎧α r b a ⎪ ⎨α r ≥1 a b

or

⎪ ⎩  K > 2 Db + 1 αb h

⎧ αb ra ⎪ ⎨ αa rb < 1 1 Db ⎪ < K < . ⎩ 2 Db + αb h 1 − ααba rra b

It is a saddle point if either (λ1 > 0, λ2 < 0, λ3 < 0) or (λ1 < 0, λ2 > 0, λ3 > 0) holds. That is,

⎧α r b a ⎪ ⎨α r ≥1 a b

⎪ ⎩ 0 < Db <  K < 2 Db + 1 αb h α r or αba ra < 1 and b

or

⎧ αb ra ⎪ ⎪ ⎨α r <1 a b

 Db 1 ⎪ ⎪ ⎩ 0 < Db <  K < min 1 − αb ra , 2Db + α h b αa rb

 K > max

 Db α ra 1− αb r a b

1

, 2Db + α h . b

Db K

> 0, which guar-

Y. Kang et al. / Applied Mathematical Modelling 51 (2017) 361–385





6. Consider the stability of the planner equilibrium EA0L = Da , 0, αraa 1 − (1−Da )K which guarantees the existence of EA0L . Since

⎛

ra −

⎜ ⎜ ⎜ =⎜ ⎜ ⎝

J |EA0L

2 r a Da ra (1 − Da /((1 −  )K )) − ( 1 −  )K 1 + αa hDa 0 ca ra (1 − Da /((1 −  )K )) 1 + αa hDa



381

(1 + αa hDa ) . Let (1 −  )K > Da > 0,

⎞ ∂ F1 αa Da |EA0L − ∂B 1 + αa hDa ⎟ ⎟ Da αb ra  ⎟ rb − 1− 0 ⎟ αa ( 1 −  )K ⎟ ⎠ ∂ F3 |EA0L 0 ∂B

we get the characteristic equation as follows



  Da αr  λ − rb − b a 1 − · αa ( 1 −  )K     2 r a Da ra (1 − Da /((1 −  )K )) ca αa ra Da (1 − Da /((1 −  )K )) 2 − + = 0. λ − λ ra − ( 1 −  )K 1 + αa hDa (1 + αa hDa )2 

α r

Thus, the eigenvalue λ2 = rb − αbaa 1 − (1−Da )K



and λ1 and λ3 are the roots of the following equation

  2 r a Da ra (1 − Da /((1 −  )K )) ca αa ra Da (1 − Da /((1 −  )K )) 2 λ − λ ra − − + = 0. ( 1 −  )K 1 + αa hDa (1 + αa hDa )2

Since (1 −  )K > Da > 0, the real parts of λ1 and λ3 have the same signs. It is easy to check that

λ2 = rb − ⇔

Da αb ra  1− <0 αa ( 1 −  )K

αa rb Da < 1 and (1 −  )K > ( > Da ) αb ra 1 − ααa rrab b

and

ra (1 − Da /((1 −  )K )) <0 1 + αa hDa 1 ⇔ ( 1 −  )K < 2 Da + . αa h

Reλ1 < 0, Reλ3 < 0 ⇔ ra −

2 r a Da

( 1 −  )K

−

Therefore, the planner equilibrium EA0L is locally asymptotically stable if

Da 1 αa rb < (1 −  )K < 2Da + α h . 1 − α ra a b

It is a source if

⎧α r a b ⎪ ⎨α r ≥1 b a

or

⎪ ⎩ ( 1 −  )K > 2 Da + 1 αa h

⎧ αa rb ⎪ ⎨ αb ra < 1 1 Da ⎪ < ( 1 −  )K < ⎩ 2 Da + αa h 1 − ααa rrab b

It is a saddle point if

⎧α r a b ⎪ ⎨α r ≥1 b a

⎪ ⎩ Da < ( 1 −  ) K < 2 Da + 1 αa h

α r

or αa rab < 1 and (1 −  )K > max b 

Da αa r 1− α r b b a

or

⎧ αa rb ⎪ ⎪α r <1 ⎨ b a

 Da 1 ⎪ ⎪ ⎩ Da < (1 −  )K < min 1 − αa rb , 2Da + αa h αb ra 

, 2Da + α1h . The proof is complete. a

Proof of Theorem 3.4 Proof. An interior equilibrium E∗ (A∗ , B∗ , L∗ ) of the system (1)-(3) has to satisfy the following three equations



ra



A αa L 1− − =0 ( 1 −  )K + db B 1 + αa hA + αb hB

(A.1)

382

Y. Kang et al. / Applied Mathematical Modelling 51 (2017) 361–385



rb 1 −



B αb L − =0 K 1 + αa hA + αb hB

(A.2)

ca αa A + cb αb B − dl = 0 1 + αa hA + αb hB

(A.3)

From the equations (A.1) and (A.2), we get



A = ((1 −  )K + db B ) 1 −

αa rb αa rb B + αb ra αb ra  K

(A.4)

which is a parabola, opening up, in (B, A)-plane. The parabola has two intersection points



MB1 =

−

( 1 −  )K db



,0 ,

MB2 =



1−

αb ra  K, 0 αa rb

with the abscissa axis A = 0 and an intersection point

 

MA = 0 , 1 −

αa rb ( 1 −  )K αb ra

with the vertical axis B = 0. From the equation (A.3), we get

A = Da −

Da B Db

(A.5)

which is a straight line in (B, A)-plane. It has an intersection point QB = (Db , 0 ) with the abscissa axis A = 0 and an intersection point QA = (0, Da ) with the vertical axis B = 0. Since E∗ (A∗ , B∗ , L∗ ) is a positive equilibrium of the system (1)-(3), (B∗ , A∗ ) must be the intersection point of the parabola (A.4) with the straight line (A.5) in the first quadrant of (B, A)-plane. We have the following four scenarios: (I) Da > 0, Db > 0. In this case, QA and QB lie on the positive half axis of B = 0 and A = 0, respectively. Thus, there exists an intersection point (B∗ , A∗ ) of the parabola (A.4)  with the straight line (A.5) in the first quadrant of (B, A)-plane if and only if α r

α r

QB lies on the right-hand of MB2 , i.e., 1 − αba ra  K < Db and QA lies on the top of the point MA , i.e., 1 − αa rab (1 −  )K < Da . b b Therefore, 1. when αb ra = αa rb , (1)-(3) has a unique interior equilibrium. 2. when α b ra < α a rb , (1)-(3) has a unique interior equilibrium if  K <

Db α ra 1− αb r

a b

3. when α b ra > α a rb , (1)-(3) has a unique interior equilibrium if (1 −  )K <

. Da αa r 1− α r b

.

Da αa r 1− α r b

. In the case that (1 −  )K >

b a

(II) Da > 0, Db < 0. In this case, QA lies on the positive half A-axis and QB lies on the negative half B-axis. Then, there exists an intersection point (B∗ , A∗ ) of the parabola (A.4) with the straight line (A.5) in the first quadrant of (B, A)-plane if α r

1 − αa rab (1 −  )K < Da . Therefore, b (1) when α b ra ≤ α a rb , then (1)-(3) has a unique interior equilibrium; (2) when α b ra > α a rb , then (1)-(3) has an interior equilibrium if (1 −  )K <

and only if QA lies on the top of MA , i.e.,

b a

Da αa r 1− α r b

,

b a

(III) Da < 0, Db > 0. In this case, QA lies on the negative half A-axis and QB lies on the positive half B-axis. Then, there exists an intersection point (B∗ , A∗ ) of the parabola  (A.4) with the straight line (A.5) in the first quadrant of (B, A)-plane if and only if MB2 lies on the right-hand of QB , i.e.,

α r

1 − αba ra  K > Db . Therefore, if α a rb > α b ra and  K > b



Db α ra 1− αb r

, then (1)-(3)

a b α r has a unique interior equilibrium E∗ (A∗ , B∗ , L∗ ). If 1 − αba ra  K < Db holds, then it is possible for the parabola (A.4) and the b

straight line (A.5) either have two positive intercepts or no intercepts. (IV) Da < 0, Db < 0. In this case, QA and QB lie on the negative half axis of B = 0 and A = 0, respectively. Therefore, (1)-(3) has no interior equilibrium. The proof is complete.  Proof of Theorem 3.5 Proof. If L = 0, Model (1)-(3) reduces to the A − B subsystem (4)-(5). According to Theorem 3.2, the omega limit set of of the A − B subsystem (4)-(5) is EAB0 = ((1 −  )K + db  K,  K, 0 ). Applying Theorem 2.5 of Hutson [42], we can conclude that dL species L is persistent in R3+ if Ldt |EAB0 > 0. Notice that

dL |E = Ldt AB0



ca αa A + cb αb B − dl 1 + αa hA + αb hB



EAB0

Y. Kang et al. / Applied Mathematical Modelling 51 (2017) 361–385

383

ca αa ((1 −  )K + db  K ) + cb αb  K − dl > 0 1 + αa h((1 −  )K + db  K ) + αb h K ⎧ D ⎪ ⎨(1 − a − db ) K < K − Da if Da > 0, Db ⇔ Da ⎪ ⎩ (1 − − db ) K > K − Da if Da < 0. Db =

Thus, species L is persistent in the full system if (1 −

(1 −

Da − db ) K > K − Da Db

Da Db

− db ) K < K − D a

when Da > 0 or

when Da < 0.

2. If B = 0, Model (1)-(3) reduces to the A − L subsystem (6)-(7). According to Theorem 3.2, the omega limit set of of the A − L subsystem (6)-(7) is EA00 = ((1 −  )K, 0, 0 ) if (1 −  )K < Da or Da < 0. Under this condition, according to Theorem 2.5 dB of Hutson [42], species B is persistent in R3+ if Bdt |EA00 > 0. Notice that

 

B dB |E = rb 1 − Bdt A00 K



−



αb L

1 + αa hA + αb hB

= rb > 0.

EA00

Therefore, species B is persistent in R3+ if (1 −  )K < Da or Da < 0. If 0 < Da < (1 −  )K < 2Da + α1h , then from Theorem 3.2, the omega limit set of the A − L subsystem of (1)-(3) is EA0L =





ra Da αa 1 − (1− )K

Da , 0,

0. Notice that



a (1 + αa hDa ) . Applying Theorem 2.5 of Hutson [42] again, species B is persistent in R3+ if

 



B αb L − K 1 + αa hA + αb hB  αa rb ⇔ 1− ( 1 −  )K < Da . αb ra

dB |E = Bdt A0L

rb 1 −

Thus, if 0 < Da < (1 −  )K < 2Da + α1h and a



α r







dB | Bdt EA0L

>

>0

EA0L

α r



1 − αa rab (1 −  )K < Da , then species B is persistent in R3+ . The inequality b

1 − αa rab (1 −  )K < Da naturally holds if α b ra ≤ α a rb , while it is equivalent to (1 −  )K < b

Da αa r 1− α r b

if α b ra > α a rb . Thus, the

b a

sufficient conditions of the persistence of species B is proved. 3. Similar to the proof of the persistence of species B above, we would provide the sufficient conditions of the persistence of species A in R3+ . Here we ignore the details.  Proof of Theorem 3.6 Proof. If  < 1 −

Da K

< 1 and 1 −

hold

1− 1−

Da K

Da Db

− db

< <1−

Da Db

− db > 1 (i.e.,

Da Db

+ db < 0) hold, then we can conclude that the following inequalities

Da . K

According to Theorem 3.1, we have

lim sup A(t ) ≤ (1 −  )K + db  K, and t→∞

From the full model (1)-(3), since the function

 L = L

ca αa A + cb αb B − dl 1 + αa hA + αb hB

Notice that if  < 1 −

Da K



< 1 and 1 −



lim sup B(t ) ≤  K. t→∞

ca αa A+cb αb B 1+αa hA+αb hB

is increasing with respect to A, B respectively, thus we have



ca αa ((1 −  )K + db  K ) + cb αb  K ≤L − dl . 1 + αa h((1 −  )K + db  K ) + αb h K Da Db

− db > 1 hold, then we have

1 − DKa Da ca αa ((1 −  )K + db  K ) + cb αb  K < dl ⇔ < 1 − . 1 + αa h((1 −  )K + db  K ) + αb h K K 1 − DDa − db b

Therefore,  < 1 −

dL Ldt

=

Da K

< 1 and 1 −

Da Db

− db > 1 hold, then we have

ca αa ((1 −  )K + db  K ) + cb αb  K ca αa A + cb αb B − dl ≤ − dl < 0 1 + αa hA + αb hB 1 + αa h((1 −  )K + db  K ) + αb h K

384

Y. Kang et al. / Applied Mathematical Modelling 51 (2017) 361–385

which implies that

lim sup L(t ) = 0. t→∞

If Da < 0 and Db < 0, then we have max {ca , cb } < dl h, this implies that

dL Ldt

=

ca αa A + cb αb B − dl = 1 + αa hA + αb hB

αa A(ca − dl h ) + αb B(cb − dl h ) − dl < 0. 1 + αa hA + αb hB

Therefore, we have lim sup L(t ) = 0 when Da < 0 and Db < 0. t→∞

If lim sup L(t ) = 0, then Model (1)-(3) is converging to the A − B subsystem (4)-(5) when time is large enough. Thus, act→∞

cording to Theorem 3.2, we can conclude that Model (1)-(3) has global stability at EAB0 when one of the following conditions hold 1. 0 <

D 1− Ka D

1− Da −db

< <1−

Da K

< 1 or

b

2. Da < 0 and Db < 0 (i.e., max {ca , cb } < dl h) .  References [1] M.C. Andersen, M. Ewald, J. Northcott, Risk analysis and management decisions for weed biological control agents: ecological theory and modeling results, Biol. Control 35 (3) (2005) 330–337. [2] H.L. Bateman, H.L. Snell, A. Chung-MacCoubrey, D.M. Finch, Growth, activity, and survivorship from three sympatric parthenogenic whiptails (family Teiidae), J. Herpetol. 44 (2) (2010) 301–306. [3] H. Bateman, T. Dudley, D. Bean, S. Ostoja, K. Hultine, M. Kuehn, A river system to watch: documenting the effects of saltcedar (Tamarix spp.) biocontrol in the Virgin river valley, Ecol. Restor. 28 (4) (2010) 405–410. [4] H. Bateman, D. Merritt, E. Glenn, P. Nagler, Indirect effects of biocontrol of an invasive riparian plant (Tamarix) alters habitat and reduces herpetofauna abundance, Biol. Invasions 17 (2014) 1–11, doi:10.1007/s10530- 014- 0707- 0. [5] H. Bateman, P. Nagler, E. Glenn, Plot-and landscape-level changes in climate and vegetation following defoliation of exotic saltcedar (Tamarix sp.) from the biocontrol agent Diorhabda carinulata along a stream in the Mojave desert (USA), J. Arid Environ. 89 (2013) 16–20. [6] H. Bateman, S. Ostoja, Invasive woody plants affect the composition of native lizard and small mammal communities in riparian woodlands, Anim. Conserv. 15 (3) (2012) 294–304. [7] B.W. Dale, L.G. Adams, R.T. Bowyer, Functional response of wolves preying on barren-ground caribou in a multiple-prey ecosystem, J. Anim. Ecol. 63 (1994) 644–652. [8] P. De Clercq, P. Mason, D. Babendreier, Benefits and risks of exotic biological control agents, Biocontrol 56 (2011) 681–698. [9] J.S. Denslow, C.M. D’Antonio, After biocontrol: assessing indirect effects of insect releases, Biol. control 35 (3) (2005) 307–318. [10] J. DiTomaso, Impact, biology, and ecology of saltcedar (Tamarix spp.) in the southwestern United States., Weed Technol. 12 (1998) 326–336. [11] J.M. DiTomaso, Invasive weeds in rangelands: species, impacts, and management, Weed Sci. 48 (2) (20 0 0) 255–265. [12] A. El-Sayed, D. Suckling, C. Wearing, J. Byers, Potential of mass trapping for long-term pest management and eradication of invasive species, J. Econ. Entomol. 99 (5) (2006) 1550–1564. [13] J. Frantzen, H. Müller-Schärer, Modeling the impact of a biocontrol agent, Puccinia lagenophorae, on interactions between a crop, Daucus carota, and a weed, Senecio vulgaris, Biol. Control 37 (3) (2006) 301–306. [14] D. Fravel, Commercialization and implementation of biocontrol, Annu. Rev. Phytopathol. 43 (2005) 337–359. [15] E. Gerber, L. Hariet, Biology, impact and interactions of potential biocontrol agents on garlic mustard, in: Proceedings of the Symposium on the Biology, Ecology, and Management of Garlic Mustard (Alliaria petiolata) and European Buckthorn (Rhamnus Cathartica), 2005. [16] M. Ghosh, A.A. Lashari, X.-Z. Li, Biological control of malaria: a mathematical model, Appl. Math. Comput. 219 (15) (2013) 7923–7939. [17] S. Gourley, Y. Lou, A mathematical model for the spatial spread and biocontrol of the Asian longhorned beetle, SIAM J. Appl. Math. 74 (3) (2014) 864–884. [18] J. Gurevitch, D.K. Padilla, Are invasive species a major cause of extinctions? Trends Ecol. Evol. 19 (9) (2004) 470–474. [19] M. Hoddle, Restoring balance: using exotic species to control invasive exotic species, Conserv. Biol. 18 (2004) 38–49. [20] R. Holt, M. Hochberg, Indirect interactions, community modules and biological control: a theoretical perspective, in: E. Wajenberg, J. Scott, P. Quimby (Eds.)Evalutating Indirect Ecological Effects of Biological Control. CABI, New York, NY. (2001) 13–37. [21] K. Hultine, J. Belnap, C. van Riper III, J. Ehleringer, P. Dennison, M. Lee, P. Nagler, K. Snyder, S. Uselman, J. West, Tamarisk biocontrol in the western United States: ecological and societal implications., Front. Ecol. Environ. 8 (2010) 467–474. [22] J. Guckenheimer, P. Holmes, Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields, Springer (1983). [23] J. Kaser, G. Heimpel, Linking risk and efficacy in biological control hostparasitoid models, Biol. Control 90 (2015) 49–60. [24] D. Kriticos, M. Watt, T. Withers, A. Leriche, M. Watson, A process based population dynamics model to explore target and non-target impacts of a biological control agent, Ecol. Model. 220 (2009) 2035–2050. [25] P. Lewis, C. DeLoach, A. Knutson, J. Tracy, T. Robbins, Biology of Diorhabda elongata deserticola (Coleoptera: Chrysomelidae), an Asian leaf beetle for biological control of saltcedars (Tamarix spp.) in the United States., Front. Ecol. Environ. 2 (2003) 101–116. [26] B. Liu, Z. Teng, L. Chen, Analysis of a predator–prey model with Holling II functional response concerning impulsive control strategy, J. Comput. Appl. Math. 193 (1) (2006) 347–362. [27] X. Liu, L. Chen, Complex dynamics of Holling type II Lotka–Volterra predator–prey system with impulsive perturbations on the predator, Chaos Solitons Fractals 16 (2) (2003) 311–320. [28] L. Lynch, A. Ives, J. Waage, M. Hochberg, M. Thomas, The risks of biocontrol: transient impacts and minimum non-target densities., Ecol. Appl. 12 (2002) 1872–1882. [29] M. Hesaaraki, S.M. Moghadas, Existence of limit cycles for predator-prey systems with a class of functional responses, Ecol. Model. 142 (2001) 1–9. [30] M. Picciano, The effect of arthropod biomass on lizard abundance, Undergraduate honor thesis, Barrett, The Honors College at Arizona State University (2014). [31] P. McEvoy, E. Coombs, Biological control of plant invaders: regional patterns, field experiments, and structured population models, Ecol. Appl. 9 (1999) 387–401. [32] R. Messing, M. Wright, Biological control of invasive species: solution or pollution? Front. Ecol. Environ. 4 (2006) 132–140. [33] W. Murdoch, J. Chesson, P. Chesson, Biological control in theory and practice, Am. Nat. 125 (1985) 344–366.

Y. Kang et al. / Applied Mathematical Modelling 51 (2017) 361–385

385

[34] K. Nelson, D. Andow, M. Banker, Problem formulation and option assessment (PFOA) linking governance and environmental risk assessment for technologies: a methodology for problem analysis of nanotechnologies and genetically engineered organisms, J. Law Med. Ethics J. Am. Soc. Law Med. Ethics 37 (2009) 732–748. [35] R.D. Parshad, E. Quansah, K. Black, M. Beauregard, Biological control via “ecological” damping: an approach that attenuates non-target effects, Math. Biosci. 273 (2016) 23–44. [36] A.F.L.A. Powell, P. Stouffer, Effects of prescribed burns and bison (Bos bison) grazing on breeding bird abundances in tallgrass prairie, Auk 123 (1) (2006) 183–197. [37] V.S.H. Rao, P.R.S. Rao, Mathematical models and stabilizing bio-control mechanisms for microbial populations in a cultured environment, Chaos Solitons Fractals 28 (5) (2006) 1222–1251. [38] N.J. Sanders, N.J. Gotelli, N.E. Heller, D.M. Gordon, Community disassembly by an invasive species, Proc. Nat. Acad. Sci. 100 (5) (2003) 2474–2477. [39] P. Shafroth, J. Cleverly, T. Dudley, C. van Riper III, J. Stuart, Control of Tamarix in the western United States: implications for water salvage, wildlife use, and riparian restoration., Environ. Manag. 35 (2005) 231–246. [40] K. Shea, D. Kelly, Estimating biocontrol agent impact with matrix models: Carduus nutans in New Zealand, Ecol. Appl. 8 (3) (1998) 824–832. [41] D. Simberloff, P. Stiling, How risky is biological control? Ecology 77 (7) (1996) 1965–1974. [42] V. Hutson, A theorem on average Liapunov functions, Mon. Math. 98 (1984) 267–275.