Contents lists available at ScienceDirect

Journal of Theoretical Biology journal homepage: www.elsevier.com/locate/yjtbi

Competition and allelopathy with resource storage: Two resources$ James P. Grover a,n, Feng-Bin Wang b a b

Department of Biology and Program in Earth and Environmental Science, University of Texas at Arlington, P.O. Box 19498, Arlington, TX 76019, United States Department of Natural Sciences in the Center for General Education, Chang Gung University, Tao-Yuan 333, Taiwan

H I G H L I G H T S

G R A P H I C A L

A B S T R A C T

We study competition for two resources with production of allelopathic toxins. Outcomes include competitive exclusion, coexistence, multiple attractors. At high resource supplies, competition by allelopathy destabilizes coexistence. Multiple attractors could make toxic algal blooms unpredictable.

art ic l e i nf o

a b s t r a c t

Article history: Received 8 October 2013 Received in revised form 6 February 2014 Accepted 10 February 2014 Available online 21 February 2014

Allelopathy is added to a familiar mathematical model of competition between two species for two essential resources in a chemostat environment. Both species store the resources, and each produces a toxin that induces mortality in the other species. The corresponding model without toxins displays outcomes of competitive exclusion independent of initial conditions, competitive exclusion that depends on initial conditions (bistability), and globally stable coexistence, depending on tradeoffs between competitors in growth requirements and consumption of the resources. Introducing toxins that act only between, and not within species, can destabilize coexistence leading to bistability or other multiple attractors. Invasibility of the missing species into a resident's semitrivial equilibrium is related to competitive outcomes. Mutual invasibility is necessary and sufﬁcient for a globally stable coexistence equilibrium, but is not necessary for coexistence at a locally stable equilibrium. Invasibility of one semitrivial equilibrium but not the other is necessary but not sufﬁcient for competitive exclusion independent of initial conditions. Mutual non-invasibility is necessary but not sufﬁcient for bistability. Numerical analysis suggests that when competitors display bistability in the absence of toxin production, increases in the overall magnitude of resource supply cause bistability to arise over a larger range of supply ratios between the two resources. When competitors display coexistence in the absence of toxin production, increases in overall resource supply destabilize coexistence and produce bistability or other conﬁgurations of multiple attractors over large ranges of supply ratios. The emergence of multiple attractors at high resource supplies suggests that blooms of harmful algae producing allelopathic toxins could be difﬁcult to predict under such rich conditions. & 2014 Elsevier Ltd. All rights reserved.

Keywords: Droop's model Harmful algae Coexistence Multiple attractors Global stability

1. Introduction ☆

Research partially supported by the National Center for Theoretical Science, National Tsing Hua University and National Council of Science, Taiwan. n Corresponding author. Tel.: þ 1 817 272 2405; fax: þ1 817 272 2855. E-mail addresses: [email protected] (J.P. Grover), [email protected], [email protected] (F.-B. Wang). http://dx.doi.org/10.1016/j.jtbi.2014.02.013 0022-5193 & 2014 Elsevier Ltd. All rights reserved.

Among the many environmental changes raising concern during recent decades is the apparent increase of harmful algal blooms in coastal and inland waters (Fu et al., 2012; Hallegraeff, 1993).

10

J.P. Grover, F.-B. Wang / Journal of Theoretical Biology 351 (2014) 9–24

Many factors are identiﬁed as contributing to such blooms, including both the absolute supply of essential nutrients such as nitrogen and phosphorus, and the relative supplies (or supply ratios) of such nutrients (Heisler et al., 2008; Smayda, 1990). A large literature documents how algal communities depend on nutrient resource conditions (Smith and Bennett, 1999; Sommer, 1989, 2002), and a long tradition of ecological theory addresses competition for nutrient resources in this context (Litchman et al., 2006; Tilman, 1981; Tilman et al., 1982). Algae that produce harmful blooms are notable for their production of toxins, which can have inhibitory, or allelopathic effects on their algal competitors (Fistarol et al., 2003; HattenrathLehman and Gobler, 2011; Zhai et al., 2013), on grazers of algae (Remmel et al., 2011; Sopanen et al., 2006; Tillman, 2003), and on organisms further removed in the food web, such as ﬁsh (Heil et al., 2001; Southard et al., 2010). Some researchers have suggested that allelopathy contributes to the formation of harmful blooms (Granéli et al., 2008; Legrand et al., 2003; but see Jonsson et al., 2009). Theory addressing competition for nutrient resources has a long history in ecology and is a mature area of research. However, theory for allelopathic interactions is currently in a stage of development. This study extends the previous theoretical work on competition by allelopathy and resource exploitation by considering competition for two resources which are stored within individuals. To our knowledge, models with allelopathy and competition for two resources have only been analyzed in the case where the competitors are plasmid-bearing and plasmid-free strains of the same species (Bhattacharyya et al., 2012). Outcomes in this case are constrained by the fact that the plasmid-bearing strain cannot persist without the plasmid-free strain also being present, because a proportion of plasmid-bearing individuals spontaneously lose the plasmid. Many models address competition between two independent species with allelopathy and competition for one resource. There are usually only two ultimate outcomes in such models, either competitive exclusion that is independent of initial conditions, or competitive exclusion that depends on initial conditions, i.e. bistability (e.g. Grover and Wang, 2013; Hsu and Waltman, 2004; Martines et al., 2009). Coexistence does not arise without additional complications, such as density-dependent regulation of toxin production by competitor populations (Braselton and Waltman, 2001), through quorum sensing for example. In the model studied here, the two resources are essential, sometimes also referred to as complementary, meaning that a minimum amount of each resource is required for population growth. Without competition by allelopathy, exploitative competition for two such resources is well understood (Li and Smith, 2007; Tilman, 1982): outcomes of competitive exclusion independent of initial conditions, coexistence, and bistability can all occur depending on parameters. In the model we studied, storage of resources is included because variations in storage provide a biological basis for the consumption patterns leading either to coexistence or bistability (Li and Smith, 2007; Turpin, 1988). Without competition by allelopathy, coexistence is associated with stable equilibria in which each species consumes relatively more of the resource mostly limiting to its own growth (León and Tumpson, 1975; Tilman, 1982). Consumption at equilibrium is directly related to storage, and thus parameters governing storage dynamics also govern stability, and whether coexistence or bistability arises. Given that both coexistence and bistability are possible when populations compete for two resources, an important question is whether and how such outcomes are modiﬁed when allelopathy also occurs. We address that question here, using a model of two competing species, each of which produces dissolved toxins that

induce mortality in the other species. As represented here, allelopathy contributes only to interspeciﬁc competition, and toxins have no intraspeciﬁc effects. Strengthening of interspeciﬁc competition is expected to destabilize coexistence (Chesson, 2000), and that was indeed found in the analyses that follow. The organization of the paper is as follows. The model system is described in next section. The model analysis is presented in Sections 3 and 4. The simulation results and biological interpretations are presented in Sections 5 and 6, respectively. Some technical proofs are deferred to the Appendix.

2. The model In this current paper, we consider the following system, which represents competition between two species in a chemostat, or continuous ﬂow environment where populations, resources, and toxins all experience the same turnover rate due to dilution by ﬂow: 8 dS > > ¼ ðSð0Þ SÞD f S1 ðS; Q S1 ÞN 1 f S2 ðS; Q S2 ÞN 2 > > > dt > > > > þ m12 ðP 2 ÞN1 Q S1 þ m21 ðP 1 ÞN 2 Q S2 ; > > > > dR > > > ¼ ðRð0Þ RÞD f R1 ðR; Q R1 ÞN 1 f R2 ðR; Q R2 ÞN2 > > dt > > > > þ m12 ðP 2 ÞN1 Q R1 þ m21 ðP 1 ÞN 2 Q R2 ; > > > > > dN 1 > > ¼ ½ð1 k1 Þ min μS1 ðQ S1 Þ; μR1 ðQ R1 Þ m12 ðP 2 Þ DN 1 ; > > > dt > > > > dQ S1 > > ¼ f S1 ðS; Q S1 Þ ð1 k1 Þ min μS1 ðQ S1 Þ; μR1 ðQ R1 Þ Q S1 ; > > dt > > > > > > dQ R1 ¼ f ðR; Q Þ ð1 k1 Þ minμ ðQ Þ; μ ðQ ÞQ ; > R1 S1 R1 R1 R1 < dt S1 R1 ð1Þ dN 2 > > ¼ ½ð1 k2 Þ min μS2 ðQ S2 Þ; μR2 ðQ R2 Þ m21 ðP 1 Þ DN 2 ; > > > dt > > > > dQ S2 > > ¼ f S2 ðS; Q S2 Þ ð1 k2 Þ min μS2 ðQ S2 Þ; μR2 ðQ R2 Þ Q S2 ; > > dt > > > > dQ R2 > > > ¼ f R2 ðR; Q R2 Þ ð1 k2 Þ min μS2 ðQ S2 Þ; μR2 ðQ R2 Þ Q R2 ; > > dt > > > > dP 1 > > ¼ k1 min μS1 ðQ S1 Þ; μR1 ðQ R1 Þ N 1 DP 1 ; > > > dt > > > > dP 2 > > > > dt ¼ k2 min μS2 ðQ S2 Þ; μR2 ðQ R2 Þ N 2 DP 2 ; > > > > Sð0Þ Z 0; Rð0Þ Z 0; P i ð0Þ Z 0; N i ð0Þ Z 0; > > > > : Q Si ð0Þ Z Q min;Si ; Q Ri ð0Þ Z Q min;Ri ; i ¼ 1; 2: Here S(t) and R(t) denote the concentrations of the limiting resources in the chemostat at time t. N i ðtÞ denotes the concentration of species i at time t. QSi (QRi) represents the amount of cell quota of resource S (R) per individual of species i at time t. μSi ðQ Si Þ and μRi ðQ Ri Þ are the growth rates of species i as a function of cell quota QSi and QRi, respectively. f Si ðS; Q Si Þ (f Ri ðR; Q Ri Þ) is the per capita uptake rate of species i as a function of resource concentration S (R) and cell quota QSi (QRi). D is the dilution rate of the chemostat. Each nutrient is supplied at the rate D, and both input concentrations are Sð0Þ and Rð0Þ . Q min;Hi denotes threshold cell quota below which no growth of species i occurs, where H ¼S, R. Growth rate for species is determined by the minimum of two Droop functions, that is, “Liebigs Law of the Minimum” is used to describe the dependence of species growth on cell quotas. This law reﬂects that the two resources are complementary, not substitutable. The functions mij ðP j Þ describe the mortality effect on species i from the toxin produced by species j. The constant ki is the fraction of production devoted to the production of the inhibitor. Hence, 0 o ki o1. The functions mij ðP j Þ satisfy mij ðP j Þ Z 0;

mij ð0Þ ¼ 0

and

m0ij ðP j Þ 4 0; 8 P j Z 0:

J.P. Grover, F.-B. Wang / Journal of Theoretical Biology 351 (2014) 9–24

Typically, mij ðP j Þ was taken as mij ðP j Þ ¼

with initial values in the domain: X ¼ fðS; R; N; Q S ; Q R ; PÞ A R6þ : Q S Z Q min;S ; Q R Z Q min;R g:

mmax;ij P j : K ij þ P j

According to Droop (1968, 1973, 1974), for H¼ S, R and i¼ 1, 2, the growth rate μHi ðQ Hi Þ takes the forms: Q μHi ðQ Hi Þ ¼ μ1;Hi 1 min;Hi ; Q Hi or

μHi ðQ Hi Þ ¼ μ1;Hi

ðQ Hi Q min;Hi Þ þ ; AHi þðQ Hi Q min;Hi Þ þ

where Q min;Hi is the minimum cell quota necessary to allow cell division and ðQ Hi Q min;Hi Þ þ is the positive part of ðQ Hi Q min;Hi Þ and μ1;Hi is the maximal growth rate at inﬁnite quotas (i.e., as Q Hi -1) of the species i. According to Grover (1992), for H¼S, R and i¼1, 2, the uptake rate f Hi ðH; Q Hi Þ takes the form: f Hi ðH; Q Hi Þ ¼ ρHi ðQ Hi Þ

H ; K Hi þ H

ρ

high high max;Hi ð max;Hi

It is easy to show that X is positively invariant for the system (2). Let Z S ¼ Sð0Þ S Q S N;

Z R ¼ Rð0Þ R Q R N:

Then (2) becomes 8 dN > > ¼ ½ð1 kÞ min μS ðQ S Þ; μR ðQ R Þ DN; > > dt > > > > > dQ S ð0Þ > > > > dt ¼ f S ðS Q S N Z S ; Q S Þ ð1 kÞ minfμS ðQ S Þ; μR ðQ R ÞgQ S ; > > > > dQ R > > > ¼ f R ðRð0Þ Q R N Z R ; Q R Þ ð1 kÞ minfμS ðQ S Þ; μR ðQ R ÞgQ R ; > > < dt dP ¼ k minfμS ðQ S Þ; μR ðQ R ÞgN DP; > > > dt > > > > dZ S > > ¼ DZ S ; > > dt > > > > > dZ R > > ¼ DZ R ; > > dt > > : Nð0Þ Z 0; Q ð0Þ Z Q S min;S ; Q R ð0Þ Z Q min;R ; Pð0Þ Z 0; Z S ð0Þ Z 0; Z R ð0Þ Z 0;

ð3Þ

where ρHi ðQ Hi Þ is deﬁned as follows:

ρHi ðQ Hi Þ ¼ ρ

11

ρ

low max;Hi Þ

with initial values in the domain:

Q Hi Q min;Hi ; Q max;Hi Q min;Hi

X ¼ fðN; Q S ; Q R ; P; Z S ; Z R Þ A R6þ : Q H Z Q min;H ;

here Q min;Hi r Q Hi rQ max;Hi . Cunningham and Nisbet (1980, 1983) took ρHi ðQ Hi Þ to be a constant. Motivated by the above classical examples, the functions μHi ðQ Hi Þ and f Hi ðH; Q Hi Þ satisfy the following assumptions: (i¼1, 2 and H ¼S, R) (C1) μHi ðQ Hi Þ is continuously differentiable for Q Hi ZQ min;Hi , μHi ðQ min;Hi Þ ¼ 0, μHi ðQ Hi Þ Z 0 and μ0Hi ðQ Hi Þ 4 0 for Q Hi ZQ min;Hi . (C2) f Hi ðH; Q Hi Þ is continuously differentiable for H 4 0 and Q Hi ZQ min;Hi , f Hi ð0; Q Hi Þ ¼ 0, f Hi ðH; Q Hi Þ Z 0, ∂f Hi ðH; Q Hi Þ= ∂H 40 and ∂f Hi ðH; Q Hi Þ=∂Q Hi r 0 for H 4 0 and Q Hi Z Q min;Hi . Two models with internal storage are closely related to the system (1). First, putting ki ¼ 0, mmax;ij ¼ 0 and P i ðtÞ 0 in (1), then we arrive at a two-species, two-resource competition model without inhibitor effects (Li and Smith, 2007). Second, in Grover and Wang (2013), the authors investigated a Droop model which involves competition for one nutrient only, with inhibitor production. We shall compare the results of those special models and our model system (1) in the Discussion section.

3. Single population model

Q H N þ Z H r H ð0Þ ; for H ¼ S; Rg: Biologically, SðtÞ ¼ Sð0Þ Q S N Z S should be nonnegative. Indeed, if there exists a t0 such that Sð0Þ Q S ðt 0 ÞNðt 0 Þ Z S ðt 0 Þ ¼ 0 then S0 ðt 0 Þ ¼ ðSð0Þ Q S N Z S Þ0 ðt 0 Þ ¼ f S ðSð0Þ Q S ðt 0 ÞNðt 0 Þ Z S ðt 0 Þ; Q S ðt 0 ÞÞNðt 0 Þ þ D½Q S ðt 0 ÞNðt 0 Þ þ Z S ðt 0 Þ ¼ DSð0Þ Z0; which implies that SðtÞ Z 0 for all t Z0. Similarly, we can show that RðtÞ ¼ Rð0Þ Q R N Z R is also nonnegative. The equations for N, QS, QR and P along with (C1) and (C2) imply that NðtÞ Z 0, Q S Z Q min;S , Q R Z Q min;R and PðtÞ Z 0 for all t Z 0. Obviously, Z S ðtÞ, Z R ðtÞ-0 as t-1. Therefore, solutions of (2) are ultimately bounded on X. By putting Z S ¼ Z R ¼ 0 in (3), we obtain the following limiting system: 8 dN > > > dt ¼ ½ð1 kÞ minfμS ðQ S Þ; μR ðQ R Þg DN; > > > > > dQ S > > > ¼ f S ðSð0Þ Q S N; Q S Þ ð1 kÞ minfμS ðQ S Þ; μR ðQ R ÞgQ S ; > > dt < dQ R ð4Þ ¼ f R ðRð0Þ Q R N; Q R Þ ð1 kÞ minfμS ðQ S Þ; μR ðQ R ÞgQ R ; > > > dt > > > > dP > > ¼ k minfμS ðQ S Þ; μR ðQ R ÞgN DP; > > > > dt : Nð0Þ Z 0; Q S ð0Þ Z Q min;S ; Q R ð0Þ Z Q min;R ; Pð0Þ Z 0; with initial values in the domain:

This section is devoted to the study of the single population model. Mathematically, it simply means that we set that ðN1 ; P 1 Þ ¼ ð0; 0Þ or ðN 2 ; P 2 Þ ¼ ð0; 0Þ in Eqs. (1) and all subscripts are dropped in the remaining equations, that is, we consider 8 dS > > ¼ ðSð0Þ SÞD f S ðS; Q S ÞN; > > dt > > > > dR > > > ¼ ðRð0Þ RÞD f R ðR; Q R ÞN; > > > dt > > > > dN > > ¼ ½ð1 kÞ minfμS ðQ S Þ; μR ðQ R Þg DN; > > dt < dQ S ¼ f S ðS; Q S Þ ð1 kÞ minfμS ðQ S Þ; μR ðQ R ÞgQ S ; > > dt > > > > > dQ R > > ¼ f R ðR; Q R Þ ð1 kÞ minfμS ðQ S Þ; μR ðQ R ÞgQ R ; > > dt > > > > dP > > > ¼ k minfμS ðQ S Þ; μR ðQ R ÞgN DP; > > dt > > : Sð0ÞZ 0; Rð0Þ Z 0; Nð0Þ Z 0; Q ð0Þ Z Q ; Q ð0Þ Z Q S

min;S

R

Y ¼ fðN; Q S ; Q R ; PÞ A R4þ : Q H Z Q min;H ; Q H N r H ð0Þ ; for H ¼ S; Rg: Eq. (4) has at most two equilibrium solutions. One of these, which we label as E 0 , corresponds to the absence of the microorganism. It is given by E 0 ¼ ðN; Q S ; Q R ; PÞ ¼ ð0; Q 0S ; Q 0R ; 0Þ;

ð2Þ

min;R ;

Pð0Þ Z 0;

where ðQ 0S ; Q 0R Þ satisﬁes ( f S ðSð0Þ ; Q 0S Þ ð1 kÞ minfμS ðQ 0S Þ; μR ðQ 0R ÞgQ 0S ¼ 0; f R ðRð0Þ ; Q 0R Þ ð1 kÞ minfμS ðQ 0S Þ; μR ðQ 0R ÞgQ 0R ¼ 0:

ð5Þ

By the same arguments as those in Li and Smith (2007, Lemma 9), the existence and uniqueness of a solution of (5) can be proved. We note that E 0 always exists. Further, it is also shown that (Li and

12

J.P. Grover, F.-B. Wang / Journal of Theoretical Biology 351 (2014) 9–24

Smith, 2007, Lemma 9) for an open and dense set of ðSð0Þ ; Rð0Þ Þ, we have μS ðQ 0S Þ a μR ðQ 0R Þ, that is, the system (4) is differentiable at E 0 . A straightforward computation shows that the Jacobian matrix at E 0 has three negative eigenvalues and the fourth eigenvalue is given by ð1 kÞ minfμS ðQ 0S Þ; μR ðQ 0R Þg D: Thus, E 0 is locally asymptotically stable if ð1 kÞ minfμS ðQ 0S Þ; μR ðQ 0R Þg D o0, and unstable if ð1 kÞ minfμS ðQ 0S Þ; μR ðQ 0R Þg D 4 0. The other possible equilibrium, labeled as E, corresponds to the presence of microorganism: E ¼ ðN^ ; Q^ S ; Q^ R ; P^ Þ;

8 ð1 kÞ minfμS ðQ^ S Þ; μR ðQ^ R Þg D ¼ 0; > > > > ð0Þ > ^ Q^ S Þ ¼ DQ^ S ; > < f S ðS Q^ S N; ð0Þ ^ Q^ R Þ ¼ DQ^ R ; f R ðR Q^ R N; > > > > k > > : P^ ¼ N^ : 1k

ð6Þ

In order to discuss the existence and uniqueness of E, we deﬁne extended real numbers sS , sR , λS , λR such that (see also Li and Smith, 2007, Section 4) D ; 1k

fS S

ð0Þ

! Rð0Þ λR ^ ^ Q S ; Q S DQ^ S ¼ 0:

sR

ð10Þ

The following result describes the asymptotic behavior of (4), whose proof is given in the Appendix. This and the following theorem extend Lemma 3.1 to demonstrate that the single-species equilibrium is globally stable, attracting any initially positive population. Theorem 3.2. The following statements are true:

where

μS ðsS Þ ¼

and

μR ðsR Þ ¼

D ; 1k

f S ðλS ; sS Þ ¼ DsS ;

f R ðλR ; sR Þ ¼ DsR :

ð7Þ

The following Lemma is concerned with the existence and uniqueness of E, whose proof is given in the Appendix. It establishes that a single-species equilibrium is feasible when the supply concentrations of both nutrients exceed the breakeven concentrations that support growth rates balancing losses, after accounting for the cost of toxin production. Feasibility of the single-species equilibrium implies instability of the trivial equilibrium: biologically, the species invades an empty habitat. Moreover, population growth at the single species equilibrium is limited by one nutrient or the other; co-limitation arises only for a biologically unlikely balance of parameters.

(i) If E 0 is locally asymptotically stable (i.e. ð1 kÞ minfμS ðQ 0S Þ; μR ðQ 0R Þg D o 0, or equivalently if and only if λS 4 Sð0Þ or λR 4 Rð0Þ ), then every solution of (4) satisﬁes lim ðNðtÞ; Q S ðtÞ; Q R ðtÞ; PðtÞÞ ¼ E 0 : t-1

(ii) If ðSð0Þ λS Þ=ðRð0Þ λR Þ a ðsS =sR Þ and E 0 is unstable (i.e. ð1 kÞ minfμS ðQ 0S Þ; μR ðQ 0R Þg D 4 0, or equivalently if and only if λS o Sð0Þ and λR o Rð0Þ ), then every solution of (4) with Nð0Þ 4 0 satisﬁes lim ðNðtÞ; Q S ðtÞ; Q R ðtÞ; PðtÞÞ ¼ E;

t-1

where E is deﬁned in (6). Next, we are able to use the theory of chain transitive sets (see, e.g., Zhao, 2003, Section 1.2) to lift the dynamics of the limiting system (4) to the full system (2). By Theorem 3.2 and the arguments similar to those in Grover and Wang (2013, Theorem 3.2), we have the following results: Theorem 3.3. The following statements are true: (i) If ð1 kÞ minfμS ðQ 0S Þ; μR ðQ 0R Þg D o0 (or equivalently if and only if λS 4 Sð0Þ or λR 4 Rð0Þ ), then every solution of (2) satisﬁes lim ðSðtÞ; RðtÞ; NðtÞ; Q S ðtÞ; Q R ðtÞ; PðtÞÞ ¼ ðSð0Þ ; Rð0Þ ; 0; Q 0S ; Q 0R ; 0Þ; t-1

Lemma 3.1. The following statements hold: (i) E 0 is stable if and only if minfμS ðQ 0S Þ; μR ðQ 0R Þg oðD=ð1 kÞÞ, or equivalently if and only if λS 4Sð0Þ or λR 4 Rð0Þ . (ii) E 0 is unstable if and only if minfμS ðQ 0S Þ; μR ðQ 0R Þg 4 ðD=ð1 kÞÞ, or equivalently if and only if λS oSð0Þ and λR o Rð0Þ . (iii) If ðSð0Þ λS Þ=ðRð0Þ λR Þ a sS =sR and E 0 is unstable, then E exists uniquely. More precisely, (a) if Rð0Þ λR Sð0Þ λS

4

sR ; sS

(ii) If ðSð0Þ λS Þ=ðRð0Þ λR Þ a ðsS =sR Þ and ð1 kÞ minfμS ðQ 0S Þ; μR ðQ 0R Þg D 4 0 (or equivalently if and only if λS o Sð0Þ and λR o Rð0Þ ), then every solution of (2) with Nð0Þ 4 0 satisﬁes ^ R; ^ N^ ; Q^ S ; Q^ R ; P^ Þ; lim ðSðtÞ; RðtÞ; NðtÞ; Q S ðtÞ; Q R ðtÞ; PðtÞÞ ¼ ðS;

t-1

where ðN^ ; Q^ S ; Q^ R ; P^ Þ is deﬁned in (6), S^ ¼ Sð0Þ Q^ S N^ ^ R^ ¼ Rð0Þ Q^ R N.

and

ð8Þ

and E 0 is unstable, then E ¼ E S ¼ ððSð0Þ λS Þ=sS ; sS ; Q^ ; ðk=ð1 kÞÞðSð0Þ λS Þ=sS Þ;

4. Two population model

R

(b) if S

ð0Þ

λS

Rð0Þ λR

4

In this section, we investigate the dynamics of the system (1). The following set is the region of interest for the system (1):

sS ; sR

Ω ¼ fðS; R; N1 ; Q S1 ; Q R1 ; N 2 ; Q S2 ; Q R2 ; P 1 ; P 2 Þ A R10 þ:

and E 0 is unstable, then E ¼ E R ¼ ððRð0Þ λR Þ=sR ; Q^ S ; sR ; ðk=ð1 kÞÞðRð0Þ λR Þ=sR Þ. Here Q^ R and Q^ S are respectively determined by ! Sð0Þ λS ^ ^ ð0Þ f R Q ;Q DQ^ ¼ 0; R

sS

R

R

R

Q Hi ZQ min;Hi ; i ¼ 1; 2; H ¼ S; Rg:

ð9Þ

It is easy to show that Ω is positively invariant for (1) and any solution of (1) with initial value in Ω exists globally on ½0; 1Þ. Let W H ðtÞ ¼ H ð0Þ H Q H1 N 1 Q H2 N 2 ;

H ¼ S; R:

ð11Þ

J.P. Grover, F.-B. Wang / Journal of Theoretical Biology 351 (2014) 9–24

Then we can rewrite (1) as follows: 8 dN 1 > > > ¼ ½ð1 k1 Þ minfμS1 ðQ S1 Þ; μR1 ðQ R1 Þg m12 ðP 2 Þ DN 1 ; > > dt > > > > dQ > S1 > ¼ f S1 ðSð0Þ Q S1 N 1 Q S2 N 2 W S ; Q S1 Þ > > > dt > > > > ð1 k1 Þ minfμS1 ðQ S1 Þ; μR1 ðQ R1 ÞgQ S1 ; > > > > > dQ R1 > > ¼ f R1 ðRð0Þ Q R1 N 1 Q R2 N 2 W R ; Q R1 Þ > > dt > > > > ð1 k1 Þ minfμS1 ðQ S1 Þ; μR1 ðQ R1 ÞgQ R1 ; > > > > > > dN 2 > > ¼ ½ð1 k2 Þ minfμS2 ðQ S2 Þ; μR2 ðQ R2 Þg m21 ðP 1 Þ DN 2 ; > > dt > > > > > dQ S2 > > ¼ f S2 ðSð0Þ Q S1 N 1 Q S2 N 2 W S ; Q S2 Þ > > < dt ð1 k2 Þ minfμS2 ðQ S2 Þ; μR2 ðQ R2 ÞgQ S2 ; > > > > dQ R2 ¼ f ðRð0Þ Q N Q N W ; Q Þ > > R R1 1 R2 2 R2 R2 > dt > > > > Þ minf μ ðQ Þ; μ ðQ ÞgQ ð1 k > 2 S2 R2 R2 ; S2 R2 > > > > dP > 1 > > ¼ k1 minfμS1 ðQ S1 Þ; μR1 ðQ R1 ÞgN 1 DP 1 ; > > > > dt > > dP 2 > > > ¼ k2 minfμS2 ðQ S2 Þ; μR2 ðQ R2 ÞgN 2 DP 2 ; > > dt > > > > dW S > > ¼ DW S ; > > > dt > > > > dW R > > ¼ DW R ; > > > dt > > : N i ð0Þ Z 0; Q Hi ð0Þ Z Q min;Hi ; P i ð0Þ Z 0; W H ð0Þ Z 0; i ¼ 1; 2; H ¼ S; R: ð12Þ with initial values in the domain:

Σ~ ¼ fðN 1 ; Q S1 ; Q R1 ; N2 ; Q S2 ; Q R2 ; P 1 ; P 2 ; W S ; W R Þ A R10 þ: Q Hi Z Q min;Hi ; Q H1 N 1 þ Q H2 N 2 þ W H rH ð0Þ ; H ¼ S; Rg: ð0Þ

ð13Þ ð0Þ

Biologically, SðtÞ ¼ S Q S1 N1 Q S2 N 2 W S and RðtÞ ¼ R Q R1 N 1 Q R2 N2 W R should be nonnegative. In fact, it is not hard to prove that SðtÞ; RðtÞ Z0 for all t Z 0. The equations for Ni, QSi, QRi and Pi, along with (C1) and (C2) imply that N i ðtÞ Z 0, Q S;i ðtÞ Z Q min;Si , Q R;i ðtÞ ZQ min;Ri and P i ðtÞ Z 0 for all t Z 0, i¼ 1, 2. Since WH satisﬁes dW H =dt ¼ DW H , H ¼S, R and then limt-1 W H ðtÞ ¼ 0. Therefore, solutions of (1) are ultimately bounded on Ω. Putting W S ¼ W R ¼ 0 in (12), then the reduced system of (1) takes the form 8 dN 1 > > ¼ ½ð1 k1 Þ minfμS1 ðQ S1 Þ; μR1 ðQ R1 Þg m12 ðP 2 Þ DN 1 ; > > > dt > > > > dQ S1 ð0Þ > > > > dt ¼ f S1 ðS Q S1 N 1 Q S2 N 2 ; Q S1 Þ ð1 k1 Þ minfμS1 ðQ S1 Þ; μR1 ðQ R1 ÞgQ S1 ; > > > > > dQ R1 > > ¼ f R1 ðRð0Þ Q R1 N 1 Q R2 N 2 ; Q R1 Þ ð1 k1 Þ minfμS1 ðQ S1 Þ; μR1 ðQ R1 ÞgQ R1 ; > > dt > > > > dN 2 > > > ¼ ½ð1 k2 Þ minfμS2 ðQ S2 Þ; μR2 ðQ R2 Þg m21 ðP 1 Þ DN 2 ; > > dt < dQ S2 ¼ f S2 ðSð0Þ Q S1 N 1 Q S2 N 2 ; Q S2 Þ ð1 k2 Þ minfμS2 ðQ S2 Þ; μR2 ðQ R2 ÞgQ S2 ; > > dt > > > > > dQ R2 > > ¼ f R2 ðRð0Þ Q R1 N 1 Q R2 N 2 ; Q R2 Þ ð1 k2 Þ minfμS2 ðQ S2 Þ; μR2 ðQ R2 ÞgQ R2 ; > > dt > > > > dP 1 > > > ¼ k1 minfμS1 ðQ S1 Þ; μR1 ðQ R1 ÞgN 1 DP 1 ; > > dt > > > > dP > 2 > ¼ k2 minfμS2 ðQ S2 Þ; μR2 ðQ R2 ÞgN 2 DP 2 ; > > > dt > > : N ð0ÞZ 0; Q ð0Þ Z Q i Hi min;Hi ; P i ð0Þ Z 0; i ¼ 1; 2; H ¼ S; R;

13

The washout equilibrium of (14), labeled as E0, corresponds to the absence of both species. It is given by E0 ¼ ðN 1 ; Q S1 ; Q R1 ; N 2 ; Q S2 ; Q R2 ; P 1 ; P 2 Þ ¼ ð0; Q 0S1 ; Q 0R1 ; 0; Q 0S2 ; Q 0R2 ; 0; 0Þ and it always exists. Here ðQ 0Si ; Q 0Ri Þ, i¼1,2, satisﬁes ( f Si ðSð0Þ ; Q 0Si Þ ð1 ki Þ minfμSi ðQ 0Si Þ; μRi ðQ 0Ri ÞgQ 0Si ¼ 0; f Ri ðRð0Þ ; Q 0Ri Þ ð1 ki Þ minfμSi ðQ 0Si Þ; μR ðQ 0Ri ÞgQ 0Ri ¼ 0: By the same arguments as those in Li and Smith (2007, Lemma 9), we can show that for an open and dense set of ðSð0Þ ; Rð0Þ Þ, we have

μSi ðQ 0Si Þ a μRi ðQ 0Ri Þ;

ð15Þ

that is, the system (14) is differentiable at E0. The next results establish the local stability of E0, whose proof is given in the Appendix. Lemma 4.1. The following statements are true: (i) E0 is locally asymptotically stable if minfμ ðQ 0 Þ; μ ðQ 0 Þg S1 R1 S1 R1 oðD=ð1 k1 ÞÞ and minfμS2 ðQ 0S2 Þ; μR2 ðQ 0R2 Þg o ðD= ð1 k2 ÞÞ. (ii) E0 is unstable if minfμ ðQ 0 Þ; μ ðQ 0 Þg 4 ðD=ð1 k1 ÞÞ or S1 R1 S1 R1 minfμS2 ðQ 0S2 Þ; μR2 ðQ 0R2 Þg 4 ðD=ð1 k2 ÞÞ. Remark 4.2. For i¼1, 2, the following statements are true: (i) minfμSi ðQ 0Si Þ; μRi ðQ 0Ri Þg o ðD=ð1 ki ÞÞ if and only if λSi 4 Sð0Þ or λRi 4 Rð0Þ ; (ii) minfμSi ðQ 0Si Þ; μRi ðQ 0Ri Þg 4 ðD=ð1 ki ÞÞ if and only if λSi o Sð0Þ and λRi o Rð0Þ . One of the semi-trivial steady-state solutions of (14), labeled as E1, corresponds to the presence of species 1 and the absence of species 2. It is given by E1 ¼ ðN 1 ; Q S1 ; Q R1 ; N 2 ; Q S2 ; Q R2 ; P 1 ; P 2 Þ ¼ ðN n1 ; Q nS1 ; Q nR1 ; 0; Q nS2 ; Q nR2 ; P n1 ; 0Þ;

where 8 ð1 k1 Þ minfμS ðQ nS1 Þ; μR1 ðQ nR1 Þg D ¼ 0; > > > > ð0Þ n n n n > > < f S1 ðS Q S1 N 1 ; Q S1 Þ ¼ DQ S1 ; f R1 ðRð0Þ Q nR1 N n1 ; Q nR1 Þ ¼ DQ nR1 ; > > > > k1 > n > Nn ; : P1 ¼ 1 k1 1 and ( f S2 ðSð0Þ Q nS1 Nn1 ; Q nS2 Þ ð1 k2 Þ minfμS2 ðQ nS2 Þ; μR2 ðQ nR2 ÞgQ nS2 ¼ 0; f R2 ðRð0Þ Q nR1 Nn1 ; Q nR2 Þ ð1 k2 Þ minfμS2 ðQ nS2 Þ; μR2 ðQ nR2 ÞgQ nR2 ¼ 0: ð16Þ

with initial values in the domain:

In order to discuss the existence, uniqueness and stability of E1, we will denote E1 by ES1 (ER1) when the species one is limited by resource S (R). At ES1, we can use the same arguments as those in the proof of Lemma 3.1 and Li and Smith (2007, Section 5) to show that Q nS1 ¼ sS1 and N n1 ¼ ðSð0Þ λS1 Þ=sS1 , that is, ! Sð0Þ λS1 k1 Sð0Þ λS1 n n n S E 1 ¼ E1 ¼ ; sS1 ; Q R1 ; 0; Q S2 ; Q R2 ; ;0 : ð17Þ sS1 1 k1 sS1

Σ ¼ fðN 1 ; Q S1 ; Q R1 ; N 2 ; Q S2 ; Q R2 ; P 1 ; P 2 Þ A R8þ :

Further, ES1 exists if λS1 o Sð0Þ , λR1 o Rð0Þ and

ð14Þ

Q Hi Z Q min;Hi ; Q H1 N 1 þ Q H2 N2 rH ð0Þ ; H ¼ S; Rg: Next, we deﬁne extended real numbers sSi , sRi , λSi , λRi for i ¼ 1; 2 such that (see also Li and Smith, 2007, Section 5):

μSi ðsSi Þ ¼

D ; 1 ki

μRi ðsRi Þ ¼

D ; 1 ki

f Si ðλSi ; sSi Þ ¼ DsSi ;

f Ri ðλRi ; sRi Þ ¼ DsRi :

Rð0Þ λR1 Sð0Þ λS1

4

sR1 : sS1

At ER1, we can show that Q nR1 ¼ sR1 and Nn1 ¼ ðRð0Þ λR1 Þ=sR1 , that is, ! Rð0Þ λR1 n k1 Rð0Þ λR1 E1 ¼ ER1 ¼ ; Q S1 ; sR1 ; 0; Q nS2 ; Q nR2 ; ;0 : sR1 1 k1 sR1

14

J.P. Grover, F.-B. Wang / Journal of Theoretical Biology 351 (2014) 9–24

Further, ER1 exists if λS1 oSð0Þ , λR1 o Rð0Þ and Sð0Þ λS1

Rð0Þ λR1

4

nn ð0Þ At ER2, we can show that Q nn λR2 Þ=sR2 , that is, R2 ¼ sR2 and N 2 ¼ ðR ! ð0Þ λR2 nn k2 Rð0Þ λR2 nn nn R R ; Q S2 ; sR2 ; 0; E2 ¼ E2 ¼ 0; Q S1 ; Q R1 ; : sR2 1 k2 sR2

sS1 : sR1

The following Lemma is concerned with the existence, uniqueness and stability of E1, whose proof is given in the Appendix. This Lemma and the following one address relationships between invasibility of the single-species equilibria, and ultimate outcomes of competition. Lemma 4.3 establishes that if the nutrient supply concentrations exceed the breakeven concentrations for species 1, then its single-species equilibrium exists and is locally stable if species 2 is unable to invade under the nutrient and toxin concentrations established at this equilibrium. A symmetrical result holds for species 2. Lemma 4.4 and Remark 4.5 demonstrate that if species 2 cannot invade the equilibrium of species 1 alone, while species 1 successfully invades the equilibrium of species 2, then there is competitive exclusion of species 2. A symmetrical result holds in the reverse case. Mutually non-invasible singlespecies equilibria imply bistability (or multiple stable attractors), while mutually invasible equilibria imply coexistence.

Further, ER2 exists if λS2 o Sð0Þ , λR2 o Rð0Þ and Sð0Þ λS2

Rð0Þ λR2

4

sS2 : sR2

By the same arguments as those in Lemma 4.3, we have the results on the existence, uniqueness and stability of E2. Lemma 4.4. The following statements are true: (i) Let

λS2 o Sð0Þ ;

λR2 o Rð0Þ ;

and

Sð0Þ λS2

Rð0Þ λR2

a

sS2 : sR2

ð21Þ

Then E2 exists uniquely. (ii) Let E2 exist. Then E2 is locally asymptotically stable if nn nn ð1 k1 Þ minfμS1 ðQ nn S1 Þ; μR1 ðQ R1 Þg m12 ðP 2 Þ D o 0;

ð22Þ

and unstable if nn nn ð1 k1 Þ minfμS1 ðQ nn S1 Þ; μR1 ðQ R1 Þg m12 ðP 2 Þ D 4 0:

Lemma 4.3. The following statements are true:

ð23Þ

(i) Let

λS1 o Sð0Þ ;

λR1 o Rð0Þ ;

and

Sð0Þ λS1 R

ð0Þ

λR1

a

sS1 : sR1

ð18Þ

Then E1 exists uniquely. (ii) Let E1 exist. Then E1 is locally asymptotically stable if ð1 k2 Þ minfμS2 ðQ nS2 Þ; μR2 ðQ nR2 Þg m21 ðP n1 Þ D o 0;

ð19Þ

and unstable if ð1 k2 Þ minfμS2 ðQ nS2 Þ; μR2 ðQ nR2 Þg m21 ðP n1 Þ D 4 0:

Remark 4.5. Assume that both E1 and E2 exist (i.e. (18) and (21) hold). We have following predictions for the system (14): (i) Competitive exclusion occurs when (19) and (23) are met or (20) and (22) are met. (ii) Bi-stability occurs when (19) and (22) are met. (iii) Coexistence occurs when (20) and (23) are met.

ð20Þ Let

Σ 0 ¼ fðN1 ; Q S1 ; Q R1 ; N2 ; Q S2 ; Q R2 ; P 1 ; P 2 Þ A Σ : N1 4 0; N2 40g; The other semi-trivial steady-state solution of (14), labeled as E2, corresponds to the presence of species 2 and the absence of species 1. It is given by nn nn nn nn nn E2 ¼ ðN 1 ; Q S1 ; Q R1 ; N 2 ; Q S2 ; Q R2 ; P 1 ; P 2 Þ ¼ ð0; Q nn S1 ; Q R1 ; N 2 ; Q S2 ; Q R2 ; 0; P 2 Þ;

where 8 nn ð1 k2 Þ minfμS2 ðQ nn > S2 Þ; μR2 ðQ R2 Þg D ¼ 0; > > > ð0Þ nn nn nn nn > < f S2 ðS Q S2 N 2 ; Q S2 Þ ¼ DQ S2 ; nn nn nn f R2 ðRð0Þ Q nn > R2 N 2 ; Q R2 Þ ¼ DQ R2 ; > > > > k2 : P nn N nn 2 ¼ 2 ;

1 k2

and ( nn nn nn nn nn f S1 ðSð0Þ Q nn S2 N 2 ; Q S1 Þ ð1 k1 Þ minfμS1 ðQ S1 Þ; μR1 ðQ R1 ÞgQ S1 ¼ 0; nn nn nn nn nn f R1 ðRð0Þ Q nn R2 N 2 ; Q R1 Þ ð1 k1 Þ minfμS1 ðQ S1 Þ; μR1 ðQ R1 ÞgQ R1 ¼ 0:

Then the following result is concerned with the possibility of the coexistence for the system (14), whose proof is given in the Appendix. This and the following theorem establish that coexistence under mutual invasibility is robust, in the form of uniform persistence where both species' populations have positive lower bounds to their long term dynamics. Theorem 4.6. Assume that both E1 and E2 exist (i.e. (18) and (21) hold). Let (20) and (23) hold. Then system (14) is uniformly persistent with respect to ðΣ 0 ; ∂Σ 0 Þ in the sense that there is an η 4 0 such that for any ðN 1 ð0Þ; Q S1 ð0Þ; Q R1 ð0Þ; N 2 ð0Þ; Q S2 ð0Þ; Q R2 ð0Þ; P 1 ð0Þ; P 2 ð0ÞÞ A Σ 0 ; the solution ðN 1 ðtÞ; Q S1 ðtÞ; Q R1 ðtÞ; N 2 ðtÞ; Q S2 ðtÞ; Q R2 ðtÞ; P 1 ðtÞ; P 2 ðtÞÞ of (14) satisﬁes lim inf N i ðtÞ Z η; t-1

We will denote E2 by ES2 resource S (R). At ES2,

(ER2)

when the species two is limited by we can show that Q nn S2 ¼ sS2 and ð0Þ N nn λS2 Þ=sS2 , that is, 2 ¼ ðS ! ð0Þ λS2 k2 Sð0Þ λS2 nn S nn ; Q ; ; s ; Q ; 0; : E2 ¼ ES2 ¼ 0; Q nn S2 S1 R1 R2 sS2 1 k2 sS2 Further, ES2 exists if λS2 o Sð0Þ , λR2 o Rð0Þ and Rð0Þ λR2 Sð0Þ λS2

4

sR2 : sS2

∂Σ 0 ≔Σ \Σ 0 :

i ¼ 1; 2:

Further, system (14) admits at least one positive (coexistence) steadystate solution. We can lift the dynamics of the limiting system (14) to the full system (1), whose proof is given in the Appendix. Theorem 4.7. Let (18), (21), (20) and (23) hold. Then system (1) admits at least one positive (coexistence) steady-state solution, and there is an η 4 0 such that for any initial value ðSð0Þ; Rð0Þ; N 1 ð0Þ; Q S1 ð0Þ; Q R1 ð0Þ; N 2 ð0Þ; Q S2 ð0Þ; Q R2 ð0Þ; P 1 ð0Þ; P 2 ð0ÞÞ A Ω

J.P. Grover, F.-B. Wang / Journal of Theoretical Biology 351 (2014) 9–24

with N 1 ð0Þ 4 0 and N 2 ð0Þ 4 0, the corresponding solution of (1) satisﬁes lim inf N i ðtÞ Z η; t-1

i ¼ 1; 2:

ð24Þ

5. Numerical simulations 5.1. Functions used For numerical simulations and biological interpretation we make the following choices for the necessary functions in the model system (1). First, f Hi ðH; Q Hi Þ ¼ ρmax;Hi

H ; K Hi þ H

so that the uptake function of species i ( ¼1,2) for resource H ( ¼ S; R) is independent of quota and has a constant maximal uptake rate (Cunningham and Nisbet, 1980). This choice simpliﬁes parameterization and interpretation; some potential errors involved were discussed by Klausmeier et al. (2004) and are summarized in the discussion section. Second, mij ðP j Þ ¼

mmax;ij P j ; K ij þ P j

as in Grover and Wang (2013). Third, Q μHi ðQ Hi Þ ¼ μ1;Hi 1 min;Hi ; Q Hi as in Droop (1968). Furthermore, we assume that μ1;Si ¼ μ1;Ri ¼ μ1;i , again to simplify parameterization and interpretation (Klausmeier et al., 2004). We note that the true maximal growth rate of a population is neither equal to the parameter μ1;i nor to this parameter discounted by the cost of inhibitor production, ð1 ki Þμ1;i because the latter rate would be obtained only in the limit of an inﬁnite quota. The maximal achievable quota of each resource is ﬁnite, constrained by the species' capability to consume resources. The analysis of the maximal growth rate presented by Klausmeier et al. (2004) applies here, with the parameter μ1;i everywhere multiplied by the factor ð1 ki Þ to account for the cost of toxin production.

15

Competitive exclusion independent of initial conditions arises when one species has the lowest breakeven concentration for both resources: e.g. λR1 o λR2 and λS1 o λS2 . In this case, species 2 goes extinct (is competitively excluded) for all nonzero initial conditions, and species 1 approaches the equilibrium it would reach as a single species growing alone (semitrivial equilibrium). Biologically, species 1 is the superior competitor for both resources, and it dominates for all supply concentrations ðRð0Þ ; Sð0Þ Þ sufﬁcient to exceed its breakeven concentrations. The other two outcomes of coexistence and bistability can arise when there is a tradeoff such that each species has the lowest breakeven concentration for a different resource: e.g. λR1 o λR2 and λS1 4 λS2 . Biologically, species 1 is a superior competitor for resource R while species 2 is a superior competitor for resource S, and a coexistence equilibrium for the two-species system exists, but is not necessarily feasible or stable. It is feasible only if the nutrient supply concentrations sufﬁcient to permit both species to have positive population densities. For extreme imbalances in nutrient supply (very high or low ratios of Rð0Þ =Sð0Þ ), the coexistence equilibrium is infeasible, and competitive exclusion occurs independent of initial conditions, the winner being the species with the lowest breakeven concentration for the resource of lower supply. For an intermediate range of resource supplies, the coexistence equilibrium is feasible and its stability then depends on relative consumption rates at the equilibrium. The coexistence equilibrium is stable (globally) when each species also consumes relatively more of the resource that limits its own growth (León and Tumpson, 1975; Tilman, 1982). If, for example λR1 o λR2 and λS1 4 λS2 , the stability condition is Cn

Cn

f R2 4 f R1

Cn

Cn

f S1 4 f S2 ;

and

where the superscript n indicates that these consumption functions are evaluated at the coexistence equilibrium. When both of the inequalities of this condition are reversed, the two-species equilibrium is unstable and there is bistability. Depending on initial conditions, one of the semitrivial, single-species equilibria is approached asymptotically, and the other species goes extinct. Since dilution and growth rates balance at equilibrium Cn

n

f R1 ¼ DQ CR1 ;

Cn

n

f R2 ¼ DQ CR2 ;

Cn

n

f S1 ¼ DQ CS1 ;

and

Cn

n

f S2 ¼ DQ CS2 ;

so the consumption for stable coexistence can be written in terms of n the equilibrium quotas Q CHi : n

n

n

n

5.2. Scenarios investigated

Q CR2 4 Q CR1

We now summarize some well-known properties of the singlespecies model that can be derived as a special case of the full system, and two-species models that can be derived as additional special cases. These properties aid in understanding the biological predictions of the full model. The one-species model has a non-trivial equilibrium when the true maximal growth rate exceeds the dilution rate, μmax;i 4 D, and when the nutrient supply concentrations are larger than the breakeven nutrient concentrations characterizing the non-trivial equilibrium, ðRð0Þ ; Sð0Þ Þ 4 ðλRi ; λSi Þ. There is a familiar graphical representation of this single-species equilibrium, based on a zero-net-growth isocline and consumption vector (Tilman, 1982). A well-known two-species, two-resource competition model arises for the special case when there are no inhibitor effects ðki ¼ 0 ; mmax;ij ¼ 0; P i ðtÞ 0Þ (e.g. Li and Smith, 2007; Tilman, 1982). Assuming that both species reach semitrivial equilibria when growing alone, then generically the associated breakeven concentrations of resources are distinct, and there are three possible asymptotic outcomes: competitive exclusion of one species independent of initial conditions, stable coexistence, and competitive exclusion of a species whose identity depends on initial conditions (bistability).

Another simpler subsystem derived from the full model system involves competition for one nutrient only, with inhibitor production (Grover and Wang, 2013). Asymptotic outcomes for this subsystem can include competitive exclusion and bistability. The latter case arises when one species is a superior competitor for the nutrient, but is both more sensitive to the inhibitor produced by the other species, and less able to produce an effective inhibitor against its competitor. In the case of such a tradeoff, bistability occurs when nutrient supply is high. When nutrient supply is low, there is competitive exclusion of the more toxic species less able to compete for nutrient. The analysis in Grover and Wang (2013) suggests that stable coexistence might also occur for some choices of functions governing nutrient uptake, growth, or mortality. But for conventional choices of these functions (such as those adopted here), stable coexistence did not appear to be possible in this one-nutrient subsystem.

and

Q CS1 4 Q CS2 :

5.3. Parameterization For parameterization of numerical models based on the full model, it is convenient to measure population densities ðN 1 ; N 2 Þ as carbon biomass concentration (μmol C L 1), and then also to measure toxins by their molar carbon concentrations (μmol C L 1).

16

J.P. Grover, F.-B. Wang / Journal of Theoretical Biology 351 (2014) 9–24

The dimensionless toxin production parameters ki then have a straightforward interpretation as the fraction of total carbon produced by a species (toxin and biomass) that is represented by toxin. Nutrient resources are measured by molar concentrations, and assumed to be nitrogen (R, μmol N L 1) and phosphorus (S, μmol P L 1), and thus quotas are the molar ratios of nutrient to carbon in phytoplankton cells. Parameter values were then chosen to represent nitrogen- and phosphorus-limited phytoplankton (Klausmeier et al., 2004; Morel, 1987; Rhee, 1974, 1978). We assign parameters representing inhibitors that are potent at low concentrations (μg to mg of toxin per liter, or ppb to ppm), have low costs to produce (e.g. contain essentially no limiting nutrients and demand o 0:1 of the cellular carbon budget, Grover et al., 2011; Grover and Wang, 2013), and induce mortality that is less than the maximal growth potential of the target species. Two parameter sets were deﬁned for the numerical studies reported here (Table 1). For the ﬁrst, called the “coexistence parameter set”, stable coexistence arises when there are no inhibitor effects, for biologically reasonable ranges of nutrient supplies. For the second, called the “bistable parameter set”, bistability arises when there are no inhibitor effects, for biologically reasonable ranges of nutrient supplies. To analyze asymptotic outcomes numerically across a wide range of nutrient supply conditions, a supply ratio gradient was examined: 2:25 r Rð0Þ =Sð0Þ r 144. Nutrient supplies were also constrained by the relation Rð0Þ Sð0Þ ¼ 16a2 , where a is a scaling factor that allows adjusting the magnitudes of supply concentrations at a ﬁxed ratio. Nutrient supplies deﬁned by these relationships include the ranges likely to occur for nitrogen and phosphorus supplies in natural aquatic habitats (Guildford and Hecky, 2000).

5.4. Coexistence We ﬁrst summarize results obtained for the coexistence parameter set. When there are no inhibitor effects ðki ¼ 0 ; mmax;ij ¼ 0; P i ðtÞ 0Þ, species 2 competitively excludes species 1 independent of initial conditions, for low supply ratios with Rð0Þ =Sð0Þ less than about 13 (Fig. 1a). Growth of both species is limited by resource R for which species 2 has a lower breakeven concentration. In contrast, species 1 competitively excludes species 2 independent of initial conditions, for high supply ratios with Rð0Þ =Sð0Þ exceeding about 22. Growth of both species is limited by resource S for which species 1 has a lower breakeven concentration. For intermediate supply ratios between about 13 and 22, coexistence is globally stable. Throughout this intermediate range, the species have a tradeoff in breakeven concentrations for different resources, growth of each species is limited by a different resource, and each species consumes resource relatively more of the resource that

limits its own growth. These are the classical outcomes of coexistence or competitive exclusion expected when there is competition for two resources (Li and Smith, 2007; Tilman, 1982). Competitive outcomes are considerably more complex when inhibitor effects are introduced into this scenario (Fig. 2b). For a sufﬁciently low magnitude of resource supply, the classical outcomes observed without inhibitor effects are preserved. In particular, coexistence is observed for intermediate supply ratios based on the tradeoffs of breakeven concentrations, growth limitation, and consumption. But as the magnitude of resource supply increases, biomasses of the competitors increase, which increases their production of toxins and alters competitive outcomes for some supply ratios. With the parameters assigned, species 1 produces more toxin than species 2, and species 2 is also more susceptible to toxin-induced mortality. Thus, for sufﬁciently high resource supply ratios, inhibitor effects reinforce the advantage that species 1 has in competing for resource S. For these high supply ratios, species 1 excludes species 2 independent of initial conditions, a fortiori. For sufﬁciently low resource supply ratios, both species are limited by resource R for which species 2 is the superior competitor based on breakeven concentrations. However, species 1 has the advantage in competition by inhibition, and with sufﬁciently high resource supply magnitude, species 1 competitively excludes species 2 for some initial conditions. For other initial conditions, species 2 competitively excludes species 1, and thus the outcome is bistable. This outcome qualitatively resembles those seen in many one-resource models with inhibition, and a tradeoff between abilities to compete for the resource and by inhibition (Defreitas and Frederickson, 1978; Grover and Wang, 2013; Hsu and Waltman, 2004). The most complex competitive outcomes arise in situations of intermediate resource supply ratio, and sufﬁciently high resource supply magnitude (labeled as “M” in Fig. 1b). Within this parameter space, there are multiple stable attractors, including at least one coexistence equilibrium that attracts some initial conditions. Numerical simulations suggest that in most of this part of parameter space, there are two locally stable attractors: the semitrivial equilibrium of species 1, and a coexistence equilibrium. For a very small part of this parameter space, three locally stable attractors were observed: the two semitrivial equilibria and a coexistence equilibrium.

5.5. Bistability We now summarize the simpler results obtained when the bistable parameter set is used. When there are no inhibitor effects ðki ¼ 0; mmax;ij ¼ 0; P i ðtÞ 0Þ, species 1 competitively excludes species 2 independent of initial conditions, for low supply ratios

Table 1 Parameter values used in numerical simulations, for the coexistence and bistable parameter sets. Parameter

Coexistence value

Bistable value

Sð0Þ

Varies

Varies

Rð0Þ D k1, k2 ρmax;R1 , ρmax;R2 ρmax;S1 , ρmax;S2 KR1, KR2 KS1, KS2 μ1;1 , μ1;2 Q min;R1 , Q min;R2 Q min;S1 , Q min;S2 K12, K21 mmax;12 , mmax;21

Varies

Varies

0.3 d 1 0.008, 0.002 8.2, 7.8 μmol N (μmol C) 1 d 1 0.49, 0.51 μmol P (μmol C) 1 d 1 7.3, 8.7 μmol N L 1 0.57, 0.53 μmol P L 1 1, 1 d 1 0.19, 0.13 μmol N (μmol C) 1 0.008, 0.012 μmol P(μmol C) 1 0.2, 0.1 μmol C L 1 0.02, 0.08 d 1

0.3 d 1 0.002, 0.008 10, 8 μmol N (μmol C) 1 d 1 0.5, 0.7 μmol P (μmol C) 1 d 1 6, 8 μmol N L 1 0.8, 0.5 μmol P L 1 1, 1 d 1 0.16, 0.2 μmol N (μmol C) 1 0.014, 0.01 μmol P (μmol C) 1 0.1, 0.2 μmol C L 1 0.3, 0.02 d 1

J.P. Grover, F.-B. Wang / Journal of Theoretical Biology 351 (2014) 9–24

1

Supply Magnitude

Supply Magnitude

1 0.8 0.6 2

C

1

0.4 0.2

0.8 0.6 1

1

2

10

10

N:P Supply Ratio

1

Supply Magnitude

1

Supply Magnitude

2

1

2

10

N:P Supply Ratio

0.8 B 1 0.4 0.2

B

0.4 0.2

10

0.6

17

M C

2 1

10

0.8 0.6 2 0.4

B 1

0.2 2

10

N:P Supply Ratio

1

10

2

10

N:P Supply Ratio

Fig. 1. Asymptotic competitive outcomes in relation to the magnitude of nutrient supply (parameter a), and the supply ratio (Rð0Þ =Sð0Þ ). Outcomes are coded: 1–species 1 excludes species 2 independent of initial conditions; 2–species 2 excludes species 1 independent of initial conditions; C–coexistence occurs independent of initial conditions; B–bistability occurs and either species 1 or species 2 is excluded depending on initial conditions; M–multiple stable attractors occur, including coexistence for at least some initial conditions. (a) Coexistence parameter set with no inhibitor effects. (b) Coexistence parameter set with inhibitor effects. (c) Bistable parameter set with no inhibitor effects. (d) Bistable parameter set with inhibitor effects.

with Rð0Þ =Sð0Þ less than about 6.5 (Fig. 1c). Growth of both species is limited by resource R for which species 1 has a lower breakeven concentration. In contrast, species 2 competitively excludes species 1 independent of initial conditions for high supply ratios with Rð0Þ =Sð0Þ exceeding about 23. Growth of both species is limited by resource S for which species 2 has a lower breakeven concentration. For intermediate supply ratios between about 6.5 and 23, bistability occurs and either species 1 or species 2 is competitively excluded depending on initial conditions. Throughout this intermediate range, the species have a tradeoff in breakeven concentrations for different resources, growth of each species is limited by a different resource, but each species consumes resource relatively more of the resource that limits the growth of the other species. These are the classical outcomes of bistability or competitive exclusion expected when there is competition for two resources (Li and Smith, 2007; Tilman, 1982). The principal effect of introducing inhibitor effects into this scenario is to expand the range of resource supply ratios under which bistability arises, as the resource supply magnitude increases (Fig. 1d). This expansion reduces the range of supply ratios for which the species inferior at competition by inhibition is able to exclude the other species. For the parameters assigned here, species 1 is inferior in competition by inhibition, so increasing resource supply magnitude reduces its ability to exclude species 2, as toxin production and inhibition effects increase. When stable coexistence is possible for the corresponding model with no inhibitor effects, outcomes of competition can involve two or more locally stable equilibria. We illustrate these results for the coexistence parameter set under a variety of resource supply conditions (Fig. 2). When one species is competitively excluded independent of initial conditions (Fig. 2a, b),

dynamics beginning far from equilibrium approach a manifold whose projection on the N1 N 2 -plane is approximately a curve connecting the two semitrivial equilibria. Dynamics then move along this manifold towards the semitrivial equilibrium that is stable, and away from the one that is unstable. When supply conditions permit stable coexistence, trajectories approach a similar manifold and then move towards the stable interior equilibrium and away from the unstable semitrivial equilibria (Fig. 2c). When supply conditions produce bistability, trajectories approach the manifold connecting the semitrivial equilibria, but then diverge from the unstable interior equilibrium, and approach one of the stable semitrivial equilibria (Fig. 2d). Some cases of multiple stable equilibria involve two stable and two unstable equilibria, with the stable equilibria being a semitrivial equilibrium and an interior equilibrium, and the unstable equilibria being the other semitrivial equilibrium and another interior equilibrium (Fig. 2e). These four equilibria lie on the manifold connecting the semitrivial equilibria. Trajectories from various initial conditions approach the manifold, and then move towards one of the stable equilibria. All the examples of this type found for the parameters explored here involve a stable semitrivial equilibrium for species 1, in addition to the stable interior equilibrium leading to coexistence for some initial conditions. An additional result with ﬁve equilibria on the manifold connecting the semitrivial equilibria arises for a small range of parameters (Fig. 2f). In these cases, the two semitrivial equilibria are stable, as is an interior equilibrium and these stable equilibria are separated by unstable interior equilibria. Once again, trajectories initially approach the manifold, and then move towards one of the stable equilibria, including the possibility of coexistence for some initial conditions.

J.P. Grover, F.-B. Wang / Journal of Theoretical Biology 351 (2014) 9–24

Biomass Species 2

18

25

25

20

20

15

15

10

10

5

5

0

0

Biomass Species 2

0

5

10

15

20

25

25

25

20

20

15

15

10

10

5

5

5

10

15

20

25

0

5

10

15

20

25

0

5 10 15 20 Biomass Species 1

25

0

0 0

Biomass Species 2

0

5

10

15

20

25

25

25

20

20

15

15

10

10

5

5 0

0 0

5 10 15 20 Biomass Species 1

25

Fig. 2. Competitive dynamics of the full model projected onto the N 1 N 2 -plane. Trajectories are shown for 80 initial conditions; solid circles show locally stable equilibria, and open circles show unstable equilibria. For all panels the coexistence parameter set was used (Table 1), with resource supplies as follows: (a) supply magnitude a ¼0.2, supply ratio Rð0Þ =Sð0Þ ¼ 10, Rð0Þ ¼ 2:530, Sð0Þ ¼ 0:2530; (b) supply magnitude a¼ 0.2, supply ratio Rð0Þ =Sð0Þ ¼ 30, Rð0Þ ¼ 4:382, Sð0Þ ¼ 0:1461; (c) supply magnitude a¼ 0.15, supply ratio Rð0Þ =Sð0Þ ¼ 16, Rð0Þ ¼ 2:4, Sð0Þ ¼ 0:15; (d) supply magnitude a¼ 0.3, supply ratio Rð0Þ =Sð0Þ ¼ 10, Rð0Þ ¼ 3:795, Sð0Þ ¼ 0:3975; (e) supply magnitude a¼ 0.2, supply ratio Rð0Þ =Sð0Þ ¼ 16, Rð0Þ ¼ 3:2, Sð0Þ ¼ 0:2; (f) supply magnitude a ¼0.3, supply ratio Rð0Þ =Sð0Þ ¼ 12:93, Rð0Þ ¼ 4:315, Sð0Þ ¼ 0:3337.

6. Discussion This study extended a classical model of competition for two essential resources that can be stored within individuals by introducing additional competition through allelopathy. The analyses address conditions under which one or both species are able to persist, that is, questions of coexistence, persistence, and competitive exclusion. Such questions are intimately linked to breakeven requirements for resources, quantifying the resources needed for population growth to balance losses. When storage occurs, two types of breakeven quantities can be deﬁned, the breakeven quota within individuals that is required for population growth to balance losses, and the corresponding breakeven concentration in the environment that is required for consumption to replenish the breakeven quota as it is apportioned to population growth. Understanding persistence of a single species is a prelude to understanding coexistence. Persistence of a single species is associated with an unstable trivial equilibrium and a stable single species equilibrium that is semitrivial with respect to the full competition model. The breakeven quotas ðsRi ; sSi Þ are deﬁned in relation to the dilution rate (D) that governs population loss.

The trivial equilibrium is unstable, and the persistence equilibrium is stable, if the quotas obtained at the supply concentrations Rð0Þ and Sð0Þ exceed the breakeven quotas. Equivalently, persistence is stable if the supply concentrations exceed the corresponding breakeven concentrations (λRi , λSi , Theorem 3.2), which are often referred to as “Rn” quantities in the ecological literature (Grover, 1997; Tilman, 1982). Persistence of a single species in the model analyzed here is essentially equivalent to persistence in the corresponding model without allelopathy (Li and Smith, 2007), and differs only in that a portion (ki) of the biomass that would otherwise contribute to population growth is apportioned to toxin production. Generically, semitrivial equilibria are unique and the growth rate of the population is limited either by resource R or S with the quota for the limiting resource approaching the breakeven quota and the concentration of the limiting resource approaching the breakeven quota. The quota and the concentration for the non-limiting resource are then determined from an implicit relationship with population size (Eq. (6)). For the full model representing competition between two species, the trivial equilibrium is unstable if at least one of the species can increase when rare. This requires that the species obtain quotas exceeding its breakeven quotas under the nutrient

J.P. Grover, F.-B. Wang / Journal of Theoretical Biology 351 (2014) 9–24

supply conditions, which in turn requires that nutrient supply concentrations exceed the corresponding breakeven concentrations. Of course, the most interesting results arise when both species have sufﬁcient resources to meet their breakeven requirements. In that case, either species can invade an empty habitat and establish a persistent population that approaches a single species, or semitrivial equilibrium. The ultimate outcome of competition can then be understood, partly, from the stability of these semitrivial equilibria to invasion by the missing species. Invasibility depends on whether the missing species obtains sufﬁcient quotas to exceed its breakeven quotas under the nutrient conditions established by the resident species, enabling population growth at a rate exceeding the loss to dilution and the mortality induced by the toxin concentration produced by the resident species. When there is mutual invasibility of both semitrivial equilibria, there is robust coexistence in the form of uniform persistence, or permanence. Each species' population will eventually increase when rare, and has a ﬁnite lower bound protecting it from extinction. In our numerical studies, mutual invasibility was always associated with a unique and globally stable coexistence equilibrium. However, we cannot rule out more complex possibilities involving multiple coexistence equilibria. For the cases of robust coexistence we observed, the classical conditions for coexistence in competition for two resources are fulﬁlled. Each species has a lower breakeven concentration for a different resource, each species' growth rate is limited by a different resource, and each species consumes relatively more of the resource that limits its own growth rate (León and Tumpson, 1975; Tilman, 1982). Essentially, when such robust coexistence arises, the enhancement of interspeciﬁc competition by allelopathic toxins is too weak to destabilize the coexistence equilibrium. Interestingly, mutual invasibility is sufﬁcient but not necessary for the existence and local stability of coexistence equilibria. Our numerical studies uncovered instances of multiple stable equilibria in which a locally stable coexistence equilibrium coexists with one or two locally stable semitrivial equilibria. We emphasize that coexistence in such cases is not likely to be robust, because real world disturbances could be large enough to displace populations from the basin of attraction for coexistence, to that for a semitrivial equilibrium. The cases of multiple stable equilibria found here also demonstrate that invasibility of semitrivial equilibria must be interpreted carefully, when the interspeciﬁc competition associated with allelopathy is sufﬁciently strong. While mutual invasibility implies robust coexistence for the model studied here, other patterns of invasibility could have complex outcomes. For example, in many other competition models, invasibility of one semitrivial equilibrium but not the other implies that one species is competitively excluded independent of initial conditions. In the model analyzed here, this pattern of invasibility is only necessary, but not sufﬁcient for unconditional competitive exclusion. In some cases, this pattern of invasibility nevertheless permits a coexistence equilibrium to be approached from some initial conditions (e.g. Fig. 2e). Similarly, mutual non-invasibility of semitrivial equilibria is necessary but not sufﬁcient for bistability in which one or the other competitor is excluded depending on initial conditions. Again, some cases of mutual invasibility permit a coexistence equilibrium to be approached for some initial conditions (e.g. Fig. 2f). When coexistence is possible in the absence of allelopathy, cases of multiple stable equilibria (with or without a locally stable coexistence equilibrium) arise only when supplies of both nutrients are sufﬁcient (e.g. Fig. 1b). Likewise, when bistability arises in the absence of allelopathy, high supplies of both nutrients expand the range of conditions under which bistability occurs (e.g. Fig. 1d). Thus high nutrient supply is generally associated with outcomes that depend on initial conditions. High nutrient supply supports higher

19

biomass, and in turn, greater production of toxins. Toxins contribute only to interspeciﬁc competition, and high toxin production thus destabilizes (at least some) coexistence equilibria, generating unstable manifolds that separate multiple stable equilibria. Ecologically, this ﬁnding suggests that when both allelopathy and resource competition occur in phytoplankton communities, population dynamics will be less predictable in eutrophic systems with high nutrient supply. In such systems, unpredictable natural disturbances can be expected to displace conditions from the basin of attraction for one equilibrium to another. Many who study harmful algal blooms have noted the potential for allelopathy (Granéli and Hansen, 2006; Granéli et al., 2008), and difﬁculty of predicting such blooms (Hallegraeff, 2010). We echo our earlier suggestion that the intrinsic dynamics of allelopathic competition could contribute to this lack of predictability (Grover and Wang, 2013). The single population model (2) is basically the same as the system (4.1) in Li and Smith (2007) except for the factor (1 k). Our mathematical analysis in Section 3 follows Li and Smith (2007, Section 4) very closely. In fact, we can also establish Theorem 3.2 by appealing to the theory of monotone dynamical systems and Zhao (2003, Theorem 2.3.2). For the two species model (1), one can derive a mass conservation and hence, we ﬁrst focus on the study of the dynamics of the limiting system (14). We are able to utilize theory of uniform persistence (e.g. Smith and Thieme, 2011; Zhao, 2003) to prove that coexistence occurs when both semitrivial equilibria are unstable (that is, (20) and (23) are met). Finally, we can lift the dynamics of the limiting system (14) to the full system (1) by using the concepts of chain transitive sets. Mathematically, it is very difﬁcult to prove that competitive exclusion occurs when one semitrivial equilibrium is (locally) asymptotically stable and the other one is unstable. We point out that the competition model without inhibitor effects in Li and Smith (2007) can be transformed into a monotone dynamical system, so the complete analysis can be obtained. There are some limitations of the models studied here, especially of their numerical realization. For simplicity, resource consumption rates depended only on nutrient concentrations (Cunningham and Nisbet, 1980), neglecting the negative dependence on nutrient quota assumed in some models (e.g Morel, 1987). The consequence of this assumption is that competitors in the model consume more of the non-limiting nutrient than is observed to be the case for phytoplankton (Klausmeier et al., 2004). Consumption of a non-limiting nutrient tends to destabilize coexistence equilibria, thus our numerical studies may be somewhat biased towards outcomes of bistability. Our numerical studies are also limited to toxins that have a low cost of production (low ki), and contain little or none of the nutrients for which competition occurs. These assumptions are reasonable for the polyketide toxins produced by many species of eukaryotic algae (Hallegraeff, 1993; Murata and Yasumoto, 2000), which contain no phosphorus and little or no nitrogen, but they would not apply to many of the toxins produced by cyanobacteria, which often contain much nitrogen (Berry et al., 2006). Moreover, the parameter ki as deﬁned here accounts only for the cost of the toxin released outside the cell, and not for costs within the cell, such as that of maintaining the synthetic machinery for toxin production. The latter costs could be represented by adjusting the parameters of the growth functions so that toxin production reduces growth. The numerical results presented are not unduly compromised by such assumptions, since the analytical results concerning stability, invasibility and the necessary and sufﬁcient conditions for the competitive outcomes illustrated would apply to a much broader set of models, choosing various functions and parameters consistent with the basic hypotheses of our theorems. Our ﬁndings of bistability and multiple stable equilibria are likely to be qualitatively robust to another biologically reasonable

20

J.P. Grover, F.-B. Wang / Journal of Theoretical Biology 351 (2014) 9–24

relaxation of assumptions. Namely, we have not allowed for the chemical degradation or inactivation of toxins, although such processes likely occur and might sometimes be rapid. If so, then a toxic population would have to produce more toxin in total to have the same effect on its competitor. In turn, this would likely increase the total nutrient supply at which multiple stable equilibria would emerge. The models we studied also assume that allelopathic toxins act interspeciﬁcally, only against other competing species, not intraspeciﬁcally. Reports of intraspeciﬁc autotoxicity are infrequent, but do exist (e.g. Olli and Trunov, 2007; Wang et al., 2013). Possibly few researchers have sought evidence of autotoxicity, but theoretically it would enhance intraspeciﬁc competition and thus tend to stabilize coexistence (McPeek, 2012). More broadly, the resource competition components of the models we study are well validated. The dependence of competitive outcomes on nutrient supplies and supply ratios has been experimentally validated with numerous laboratory culture experiments with algae (Sommer, 1983, 2002; Tilman, 1981) and is supported by broad observations of natural phytoplankton communities (Sommer, 1989, 1990). But the allelopathy components of the models we study lack such foundations. We assumed that toxin production is proportional to biomass production, but alternative formulations have been proposed (Chakraborty et al., 2008; Grover et al., 2011; Martines et al., 2009). Likewise, we assumed that toxins are fatal, inducing a mortality rate that saturates with concentration, while again alternative modes of toxic effects have been proposed, such as growth inhibition (Hsu and Waltman, 2004). There is no consensus on such matters among theoreticians who model allelopathy, nor do available experimental results strongly constrain the likely alternatives. The development of a mechanistic theory of resource competition did much to enhance our understanding of how nutrients inﬂuence phytoplankton communities (Sommer, 2002; Tilman et al., 1982), and has helped inform our understanding of how to manage eutrophication problems (Smith et al., 1999; Sterner, 2008). Our theoretical and empirical understanding of allelopathic interactions among phytoplankton is less advanced, but we suggest that its development will advance our understanding of harmful algal blooms and their management.

Appendix In this appendix, we provide the technical proofs of some Lemmas and Theorems. Proof of Lemma 3.1. We ﬁrst prove Parts (i) and (ii). It is not hard to see that E 0 is unstable if and only if minfμS ðQ 0S Þ; μR ðQ 0R Þg 4ðD=ð1 kÞÞ, or equivalently if and only if Q 0S 4 sS and Q 0R 4 sR . In view of (5), it follows that f S ðSð0Þ ; sS Þ DsS 4 0 and f R ðRð0Þ ; sR Þ DsR 40, that is, λS o Sð0Þ and λR o Rð0Þ . Parts (i) and (ii) are proved. For Part (iii), we discuss the following two cases: Case 1: If the population is limited by resource S at the equilibrium, then μS ðQ S Þ o μR ðQ R Þ at the equilibrium. From the ﬁrst equation of (6), it follows that Q^ S ¼ sS . From the second equation of (6), it follows that N^ ¼ ðSð0Þ λS Þ=sS . Note that N^ 40 if and only if Sð0Þ 4 λS . In this case, Q^ R is uniquely determined by (9). On the other hand, μR ðQ^ R Þ 4 μS ðQ^ S Þ ¼ D=ð1 kÞ ¼ μR ðsR Þ, that is, ^ Q R 4 sR . Thus, Q^ R 4 sR if and only if f R Rð0Þ

Sð0Þ λS

sS

!

sR ; sR DsR 4 f R Rð0Þ

! Sð0Þ λS ^ ^ Q R ; Q R DQ^ R ¼ 0

sS

In view of (7), it follows that Q^ R 4 sR if and only if Rð0Þ ððSð0Þ λS Þ=sS ÞsR 4 λR , or equivalently if and only if (8) holds. Therefore, E ¼ E S ¼ ððSð0Þ λS Þ=sS Þ; sS ; Q^ R ; ðk=ð1 kÞÞðSð0Þ λS Þ=sS Þ exists if λS o Sð0Þ , λR o Rð0Þ and (8) holds.

Case 2: If the population is limited by resource R at the equilibrium, then μR ðQ R Þ o μS ðQ S Þ at the equilibrium. By the same arguments as those in Case 1, it follows that E ¼ E R ¼ ððRð0Þ λR Þ=sR ÞÞ; Q^ S ; sR ; ðk=ð1 kÞÞððRð0Þ λR Þ=sR Þ exists if λS o Sð0Þ , λR o Rð0Þ and ðSð0Þ λS Þ=ðRð0Þ λR Þ 4 sS =sR , where Q^ S is uniquely determined by (10). □ Proof of Theorem 3.2. The ﬁrst three equations in (4) are independent of P, so we ﬁrst study the following system: 8 dN > > ¼ ½ð1 kÞ minfμS ðQ S Þ; μR ðQ R Þg DN; > > > dt > > > > < dQ S ¼ f ðSð0Þ Q N; Q Þ ð1 kÞ minfμ ðQ Þ; μ ðQ ÞgQ ; S S S R S S S R dt ð25Þ > > dQ R > ð0Þ > > ¼ f R ðR Q R N; Q R Þ ð1 kÞ minfμS ðQ S Þ; μR ðQ R ÞgQ R ; > > dt > > : Nð0Þ Z 0; Q ð0Þ Z Q ; Q ð0Þ Z Q ; S

min;S

R

min;R

with initial values in the domain fðN; Q S ; Q R Þ A R3þ : Q H Z Q min;H ; Q H N r H ð0Þ ; for H ¼ S; Rg: One can use the same arguments as those in Li and Smith (2007, Section 4) to determine the long-term behavior of the system (25): (a) If ð1 kÞ minfμS ðQ 0S Þ; μR ðQ 0R Þg D o0, then every solution of (25) satisﬁes lim ðNðtÞ; Q S ðtÞ; Q R ðtÞÞ ¼ ð0; Q 0S ; Q 0R Þ;

t-1

where ðQ 0S ; Q 0R Þ satisﬁes (5). (b) If ðSð0Þ λS Þ=ðRð0Þ λR Þ a ðsS =sR Þ and ð1 kÞ minfμS ðQ 0S Þ; μR ðQ 0R Þg D 4 0, then every solution of (25) with Nð0Þ 4 0 satisﬁes ^ Q^ ; Q^ Þ; lim ðNðtÞ; Q S ðtÞ; Q R ðtÞÞ ¼ ðN; S R

t-1

^ Q^ S ; Q^ R Þ satisﬁes the ﬁrst three equations in (6). where ðN; For the case (a), the equation of P in (4) is asymptotic to dP=dt ¼ DP. This implies that lim PðtÞ ¼ 0

t-1

by the theory for asymptotically autonomous semiﬂows (see, e.g., Thieme, 1992, Corollary 4.3). For the case (b), the equation of P in (4) is asymptotic to n o dP ¼ k min μS ðQ^ S Þ; μR ðQ^ R Þ N^ DP: dt This implies that n o k k ^ N: lim PðtÞ ¼ min μS ðQ^ S Þ; μR ðQ^ R Þ N^ ¼ t-1 D 1k The proof is complete.

□

Proof of Lemma 4.1. The local stability of E0 is determined by the Jacobian matrix of (14) at E0, denoted by 0 1 a11 0 0 0 0 0 0 0 B C 0 0 0 0 C B a21 a22 a23 a24 B C B a31 a32 a33 a34 0 0 0 0 C B C B C 0 0 0 0 C 0 0 a44 B 0 C; J0 ¼ B B a51 0 0 a54 a55 a56 0 0 C B C B C 0 0 a64 a65 a66 0 0 C B a61 B C B a71 0 0 0 0 0 D 0 C @ A 0 0 0 D 0 0 0 a84 where n o a11 ¼ ð1 k1 Þ min μS1 ðQ 0S1 Þ; μR1 ðQ 0R1 Þ D;

J.P. Grover, F.-B. Wang / Journal of Theoretical Biology 351 (2014) 9–24

∂f S1 ðSð0Þ ; Q 0S1 Þ ; ∂S ∂f ðRð0Þ ; Q 0R1 Þ ∂f ðSð0Þ ; Q 0S2 Þ ; a51 ¼ Q 0S1 S2 ; a31 ¼ Q 0R1 R1 ∂R ∂S 0 ð0Þ n o ∂f ðR ; Q R2 Þ ; a71 ¼ k1 min μS1 ðQ 0S1 Þ; μR1 ðQ 0R1 Þ ; a61 ¼ Q 0R1 R2 ∂R ∂f S1 ðSð0Þ ; Q 0S1 Þ ∂½minfμS1 ðQ S1 Þ; μR1 ðQ R1 ÞgQ S1 a22 ¼ ð1 k1 Þ 0 0 o 0; ∂Q S1 ∂Q S1 ðQ S1 ;Q R1 Þ ∂ minf μ ðQ Þ; μ ðQ Þg S1 R1 S1 R1 a32 ¼ ð1 k1 ÞQ 0R1 0 0 r 0; ∂Q S1 ðQ ;Q Þ a21 ¼ Q 0S1

S1

∂ minfμS1 ðQ S1 Þ; μR1 ðQ R1 Þg a23 ¼ ð1 k1 ÞQ 0S1 ∂Q R1

a33 ¼

R1

ðQ 0S1 ;Q 0R1 Þ

r 0;

0 S1 ;Q R1 Þ

o 0;

∂f S1 ðSð0Þ ; Q 0S1 Þ ; ∂S n

∂f ðSð0Þ ; Q 0R1 Þ a34 ¼ Q 0R2 R1 ; ∂R o a44 ¼ ð1 k2 Þmin μS2 ðQ 0S2 Þ; μR2 ðQ 0R2 Þ D; ∂f S2 ðSð0Þ ; Q 0S2 Þ ; ∂S 0 ð0Þ ∂f ðR ; Q R2 Þ ; a64 ¼ Q 0R2 R2 ∂R

n o a84 ¼ k2 min μS2 ðQ 0S2 Þ; μR2 ðQ 0R2 Þ ;

∂f ðSð0Þ ; Q 0S2 Þ ∂½minfμS2 ðQ S2 Þ; μR2 ðQ R2 ÞgQ S2 a55 ¼ S2 ð1 k2 Þ 0 0 o 0; ∂Q S2 ∂Q S2 ðQ S2 ;Q R2 Þ ∂ minf μ ðQ Þ; μ ðQ Þg S2 R2 S2 R2 a65 ¼ ð1 k2 ÞQ 0R2 0 0 r 0; ∂Q S2 ðQ ;Q Þ S2 R2 0 ∂ minfμS2 ðQ S2 Þ; μR2 ðQ R2 Þg a56 ¼ ð1 k2 ÞQ S2 0 0 r 0; ∂Q R2 ðQ S2 ;Q R2 Þ ∂f R2 ðRð0Þ ; Q 0R2 Þ ∂½minfμS2 ðQ S2 Þ; μR2 ðQ R2 ÞgQ R2 a66 ¼ ð1 k2 Þ 0 0 o 0: ∂Q R2 ∂Q R2 ðQ ;Q Þ S2

R2

From (15), it follows that either a23 or a32 must be equal to 0. Thus, it is easy to see that a11 and a44 determine the stability of E0. □ Proof of Lemma 4.3. From our previous discussions before the statements of Lemma 4.3, Part (i) is true. For Part (ii), we only consider the case where E1 ¼ ES1 , where ES1 is deﬁned in (17). The local stability of ES1 is determined by the Jacobian matrix of (14) at ES1, denoted by 0 1 0 c12 0 0 0 0 0 m012 ð0ÞN n1 B C 0 0 0 B c21 c22 0 c24 0 C B C B c31 c32 c33 c34 0 C 0 0 0 B C B C B C 0 0 0 0 0 0 0 c 44 S C; J1 ¼ B Bc C c 0 c c c 0 0 55 52 54 56 B 51 C B C B c61 0 c63 c64 c65 c66 C 0 0 B C Bc C c 0 0 0 0 D 0 72 @ 71 A 0 0 D 0 0 0 c84 0 c11 ¼ ð1 k1 Þ min μS1 ðQ nS1 Þ; μR1 ðQ nR1 Þ D ¼ ð1 k1 ÞμS1 ðQ nS1 Þ D ¼ 0; ∂f S1 ðSð0Þ Q nS1 N n1 ; Q nS1 Þ ; ∂S n n ð0Þ n ∂f ðR Q R1 N 1 ; Q R1 Þ ; c31 ¼ Q nR1 R1 ∂R ∂f ðSð0Þ Q nS1 N n1 ; Q nS2 Þ ; c51 ¼ Q nS1 S2 ∂S n n ð0Þ n ∂f ðR Q R1 N 1 ; Q R2 Þ ; c61 ¼ Q nR1 R2 ∂R n 0 c71 ¼ k1 μS1 ðQ S1 Þ; c12 ¼ ð1 k1 ÞμS1 ðQ nS1 ÞNn1 ; c21 ¼ Q nS1

ð1 k1 Þ½μS1 ðQ nS1 Þ þ μ0S1 ðQ nS1 ÞQ nS1 o0; c32 ¼ ð1 k1 ÞQ nR1 μ0S1 ðQ nS1 Þ;

c52 ¼ N n1

c72 ¼ k1 Nn1 μ0S1 ðQ nS1 Þ;

∂f S2 ðSð0Þ Q nS1 N n1 ; Q nS2 Þ ; ∂S

∂f R1 ðRð0Þ Q nR1 Nn1 ; Q nR1 Þ ∂f R1 ðRð0Þ Q nR1 N n1 ; Q nR1 Þ þ ∂R ∂Q R1 n ð1 k1 ÞμS1 ðQ S1 Þ o 0;

c33 ¼ N n1

∂f R2 ðRð0Þ Q nR1 Nn1 ; Q nR2 Þ ; ∂R ð0Þ n n n ∂f ðS Q S1 N1 ; Q S1 Þ c24 ¼ Q nS2 S1 ; ∂S ∂f ðRð0Þ Q nR1 N n1 ; Q nR1 Þ ; c34 ¼ Q nR2 R1 ∂R c44 ¼ ð1 k2 Þ min μS2 ðQ nS2 Þ; μR2 ðQ nR2 Þ m21 ðP n1 Þ D; ∂f S2 ðSð0Þ Q nS1 Nn1 ; Q nS2 Þ ; ∂S n n ð0Þ n ∂f ðR Q R1 N 1 ; Q R2 Þ ; c64 ¼ Q nR2 R2 ∂R n n c84 ¼ k2 min μS2 ðQ S2 Þ; μR2 ðQ R2 Þ ;

c54 ¼ Q nS2

a54 ¼ Q 0S2

where

∂f S1 ðSð0Þ Q nS1 Nn1 ; Q nS1 Þ ∂f S1 ðSð0Þ Q nS1 N n1 ; Q nS1 Þ þ ∂S ∂Q S1

c63 ¼ N n1

∂f R1 ðRð0Þ ; Q 0R1 Þ ∂½minfμS1 ðQ S1 Þ; μR1 ðQ R1 ÞgQ R1 ð1 k1 Þ 0 ∂Q R1 ∂Q R1 ðQ

a24 ¼ Q 0S2

c22 ¼ N n1

21

c55 ¼

∂f S2 ðSð0Þ Q nS1 N n1 ; Q nS2 Þ ∂Q S2 ∂½minfμS2 ðQ S2 Þ; μR2 ðQ R2 ÞgQ S2 ð1 k2 Þ n ∂Q S2 ðQ

n

S2 ;Q R2 Þ

o 0;

∂½minfμS2 ðQ S2 Þ; μR2 ðQ R2 Þg n n ; ∂Q S2 ðQ ;Q Þ S2 R2 n ∂½minfμS2 ðQ S2 Þ; μR2 ðQ R2 Þg c56 ¼ ð1 k2 ÞQ S2 n n ; ∂Q c65 ¼ ð1 k2 ÞQ nR2

R2

ðQ S2 ;Q R2 Þ

∂f ðRð0Þ Q nR1 N n1 ; Q nR2 Þ c66 ¼ R2 ∂Q R2 ∂½minfμS2 ðQ S2 Þ; μR2 ðQ R2 ÞgQ R2 ð1 k2 Þ n ∂Q R2 ðQ

n

S2 ;Q R2 Þ

o0:

Since only one nutrient is limiting, it follows that either c56 or c65 must be equal to 0. Thus, the stability of ES1 is determined by c44 ¼ ð1 k2 ÞminfμS2 ðQ nS2 Þ; μR2 ðQ nR2 Þg m21 ðP n1 Þ D:

□

Proof of Theorem 4.6. Suppose Ψ t : Σ -Σ are the solution semiﬂows associated with system (14), that is,

Ψ t ðN 1 ð0Þ; Q S1 ð0Þ; Q R1 ð0Þ; N 2 ð0Þ; Q S2 ð0Þ; Q R2 ð0Þ; P 1 ð0Þ; P 2 ð0ÞÞ ¼ ðN 1 ðtÞ; Q S1 ðtÞ; Q R1 ðtÞ; N 2 ðtÞ; Q S2 ðtÞ; Q R2 ðtÞ; P 1 ðtÞ; P 2 ðtÞÞ; where x≔ðN 1 ð0Þ; Q S1 ð0Þ; Q R1 ð0Þ; N 2 ð0Þ; Q S2 ð0Þ; Q R2 ð0Þ; P 1 ð0Þ; P 2 ð0ÞÞ A Σ . Let ωðxÞ be the omega-limit set of the orbit of Ψt with initial values x A Σ . It is easy to see that Ψ t ðΣ 0 Þ Σ 0 . Since solutions of the system (14) are ultimately bounded, it follows that Ψt is point dissipative and compact. Recall that E0, E1 and E2 are ﬁxed points of Ψt. Further, fE0 g, fE1 g and fE2 g are pairwise disjoint, compact and isolated invariant sets for Ψt in ∂Σ0. We are going to show the following property: ⋃ ωðxÞ fE0 ; E1 ; E2 g:

ð26Þ

x A ∂Σ 0

In the case where N 1 ð0Þ 4 0 and N 2 ð0Þ ¼ 0, we have N1 ðtÞ 4 0 and N2 ðtÞ ¼ 0, 8 t Z 0. Then the equation for P 2 ðtÞ in (14) satisﬁes the following system: dP 2 ¼ DP 2 ; dt

and hence; lim P 2 ðtÞ ¼ 0: t-1

22

J.P. Grover, F.-B. Wang / Journal of Theoretical Biology 351 (2014) 9–24

Thus, ðN 1 ðtÞ; Q S1 ðtÞ; Q R1 ðtÞ; P 1 ðtÞÞ is asymptotic to the following system: 8 dN 1 > > ¼ ½ð1 k1 Þ minfμS1 ðQ S1 Þ; μR1 ðQ R1 Þg DN 1 ; > > > dt > > > > dQ S1 > > ¼ f S1 ðSð0Þ Q S1 N 1 ; Q S1 Þ ð1 k1 Þ minfμS1 ðQ S1 Þ; μR1 ðQ R1 ÞgQ S1 ; < dt > > dQ R1 ¼ f ðRð0Þ Q N 1 ; Q Þ ð1 k1 Þ minfμ ðQ Þ; μ ðQ ÞgQ ; > R1 R1 S1 R1 R1 R1 > S1 R1 > dt > > > > dP > 1 > ¼ k1 minfμS1 ðQ S1 Þ; μR1 ðQ R1 ÞgN 1 DP 1 ; : dt

Proof of Theorem 4.7. Since systems (1) and (12) are equivalent, it sufﬁces to study the system (12). Let

From (18) and Remark 4.2, it follows that ð1 k1 Þ minfμS1 ðQ 0S1 Þ; μR1 ðQ 0R1 Þg D 4 0. By Theorem 3.2 and the theory for asymptotically autonomous semiﬂows (see, e.g., Thieme, 1992, Corollary 4.3), it follows that n

Σ~ 0 ¼ fðN1 ; Q S1 ; Q R1 ; N2 ; Q S2 ; Q R2 ; P 1 ; P 2 ; W S ; W R Þ A Σ~ : N1 40; N2 4 0g; ∂Σ~ 0 ≔Σ~ \Σ~ 0 ; where Σ~ is deﬁned in (13). Assume that Ψ~ t : Σ~ -Σ~ are the solution semiﬂows associated with system (12), that is,

Ψ~ t ðN1 ð0Þ; Q S1 ð0Þ; Q R1 ð0Þ; N2 ð0Þ; Q S2 ð0Þ; Q R2 ð0Þ; P 1 ð0Þ; P 2 ð0Þ; W S ð0Þ; W R ð0ÞÞ

n

¼ ðN 1 ðtÞ; Q S1 ðtÞ; Q R1 ðtÞ; N 2 ðtÞ; Q S2 ðtÞ; Q R2 ðtÞ; P 1 ðtÞ; P 2 ðtÞ; W S ðtÞ; W R ðtÞÞ;

lim ðN1 ðtÞ; Q S1 ðtÞ; Q R1 ðtÞ; P 1 ðtÞÞ ¼ ðN 1 ; Q S1 ; Q R1 ; P 1 Þ:

t-1

Then, the equations for Q S2 ðtÞ and Q R2 ðtÞ in (14) are asymptotic to 8 dQ S2 > > ¼ f S2 ðSð0Þ Q nS1 N n1 ; Q S2 Þ ð1 k2 Þ minfμS2 ðQ S2 Þ; μR2 ðQ R2 ÞgQ S2 ; < dt dQ R2 > > : ¼ f R2 ðRð0Þ Q nR1 N n1 ; Q R2 Þ ð1 k2 Þ minfμS2 ðQ S2 Þ; μR2 ðQ R2 ÞgQ R2 : dt

ð27Þ By the same arguments as those in Li and Smith (2007, Lemma 9), we can show that ðQ nS2 ; Q nR2 Þ is the unique equilibrium of the system (27), where ðQ nS2 ; Q nR2 Þ is deﬁned in (16). It is easy to see that the solution semiﬂows associated with system (27) is monotone with respect to the partial order r K (see, e.g., Smith, 1995), which is induced by the positive cone K≔R þ R in R2 . Consequently, we can conclude that ðQ nS2 ; Q nR2 Þ is globally attractive for (27) by appealing to the theory of monotone systems. From the theory for asymptotically autonomous semiﬂows (see, e.g., Thieme, 1992, Corollary 4.3), it follows that n

ð28Þ

k minfμSi ðQ~ Si Þ; μRi ðQ~ Ri ÞgN~ i 4 0; i ¼ 1; 2: P~ i ¼ i D It then follows that ðN~ 1 ; Q~ S1 ; Q~ R1 ; N~ 2 ; Q~ S2 ; Q~ R2 ; P~ 1 ; P~ 2 Þ is a positive steady-state solution for (14). This completes the proof. □

fðN 1 ; Q S1 ; Q R1 ; P 1 Þ A R4þ : Q H Z Q min;H1 ; Q H1 N1 r H ð0Þ ; for H ¼ S; Rg:

n

Ψt in Σ0 and Ψt has at least one ﬁxed point:

ðN~ 1 ; Q~ S1 ; Q~ R1 ; N~ 2 ; Q~ S2 ; Q~ R2 ; P~ 1 ; P~ 2 Þ A Σ 0 ; In fact,

with initial values in the domain

n

attractor A0 for

n

lim ðQ S2 ðtÞ; Q S2 ðtÞÞ ¼ ðQ S2 ; Q R2 Þ:

where ~ x≔ðN 1 ð0Þ; Q S1 ð0Þ; Q R1 ð0Þ; N 2 ð0Þ; Q S2 ð0Þ; Q R2 ð0Þ; P 1 ð0Þ; P 2 ð0Þ; W S ð0Þ; W R ð0ÞÞ A Σ~ . We also recall that Ψ t : Σ -Σ are the solution semi~ be the omega-limit ~ ≔ω ~ ðxÞ ﬂows associated with system (14). Let ω set of the orbit of Ψ~ t with initial values x~ A Σ~ . From the ninth and tenth equations of the system (12), it follows that lim W S ðtÞ ¼ lim W R ðtÞ ¼ 0:

t-1

t-1

~ ¼ I fð0; 0Þg: Since Σ~ Thus, there exists a set I R8þ such that ω ~ Σ~ . For any given ðN 1 ; Q S1 ; is closed, it follows that ω Q R1 ; N 2 ; Q S2 ; Q R2 ; P 1 ; P 2 Þ A I , we have ðN 1 ; Q S1 ; Q R1 ; N 2 ; Q S2 ; Q R2 ; ~ Σ~ . By the deﬁnition of Σ~ , it follows that P 1 ; P 2 ; 0; 0Þ A ω ðN 1 ; Q S1 ; Q R1 ; N 2 ; Q S2 ; Q R2 ; P 1 ; P 2 Þ A Σ . Thus, I Σ . ~ is a compact, invariant and By Zhao (2003, Lemma 1.2.10 ), ω 0 internal chain transitive set for Ψ~ t . Moreover, if ϕ ≔ðN 01 ; Q 0S1 ; 0 0 0 0 0 0 0 8 ~ , there holds Q R1 ; N 2 ; Q S2 ; Q R2 ; P 1 ; P 2 Þ A R þ with ðϕ ; 0; 0Þ A ω

Ψ~ t ∣ω~ ðϕ0 ; 0; 0Þ ¼ ðΨ t ðϕ0 Þ; 0; 0Þ;

In the case where N 1 ð0Þ ¼ 0 and N2 ð0Þ 4 0, we can use the similar arguments to show that

where Ψ t ðϕ Þ are the solution semiﬂows associated with (14) on Σ . It then follows from the deﬁnition of internally chain transitive sets that I is a compact, invariant and internal chain transitive set for Ψ t : Σ -Σ . In order to use Zhao (2003, Theorem 1.3.1) with L ¼ I , we ﬁrst verify that I 2 = ffE0 g; fE1 g; fE2 gg. Suppose, by contradiction, that I ¼ fE1 g, then

lim Ψ t ðN 1 ð0Þ; Q S1 ð0Þ; Q R1 ð0Þ; N 2 ð0Þ; Q S2 ð0Þ; Q R2 ð0Þ; P 1 ð0Þ; P 2 ð0ÞÞ ¼ E2 :

~ ω~ ¼ ðE1 ; 0; 0Þ≔E:

t-1

It then follows that lim Ψ t ðN 1 ð0Þ; Q S1 ð0Þ; Q R1 ð0Þ; N 2 ð0Þ; Q S2 ð0Þ; Q R2 ð0Þ; P 1 ð0Þ; P 2 ð0ÞÞ ¼ E1 :

t-1

t-1

In the case where N 1 ð0Þ ¼ 0 and N 2 ð0Þ ¼ 0, we can also show that lim Ψ t ðN 1 ð0Þ; Q S1 ð0Þ; Q R1 ð0Þ; N 2 ð0Þ; Q S2 ð0Þ; Q R2 ð0Þ; P 1 ð0Þ; P 2 ð0ÞÞ ¼ E0 :

t-1

Consequently, Ψ t : Σ -Σ satisﬁes the property (26). It is obvious that no subset of fE0 ; E1 ; E2 g forms a cycle in ∂Σ0. By (18), (21), (20) and (23) and the same arguments to those in Zhao (2003, Lemma 5.1.1), there exists a δ 4 0 such that each Ej is a uniform weak repeller for Σ0 in the sense that lim sup J Ψ t ðxÞ Ej J Z δ; t-1

for any x A Σ 0 : Therefore, each Ej is isolated in Σ and W s ðEj Þ \ Σ 0 ¼ ∅, where W s ðEj Þ is the stable set of Ej (see Zhao, 2003). Since Ψ t : Σ -Σ is point dissipative and compact, we conclude from Zhao (2003, Theorem 1.1.3) that there exists a global attractor A for Ψt in Σ . By Zhao (2003, Theorem 1.3.1) on strong repellers, Ψ t : Σ -Σ is uniformly persistent with respect to ðΣ 0 ; ∂Σ 0 Þ. It follows from Zhao (2003, Theorem 1.3.6) that there exists a global

0

Thus, we have lim Ψ~ t ðN 01 ; Q 0S1 ; Q 0R1 ; N02 ; Q 0S2 ; Q 0R2 ; P 01 ; P 02 ; W 0S ; W 0R Þ ¼ ðE1 ; 0; 0Þ:

t-1

From this, we have that lim J ðQ S2 ðtÞ; Q R2 ðtÞÞ ðQ nS2 ; Q nR2 Þ J ¼ 0;

t-1

lim ∣P 1 ðtÞ P n1 ∣ ¼ 0:

t-1

Since minfμS2 ð♠Þ; μR2 ð♣Þg and m21 ðÞ are continuous functions, for any ϵ 4 0, there is a T 4 0 such that for all t Z T ϵ n n ; ∣min μS2 ðQ S2 ðtÞÞ; μR2 ðQ R2 ðtÞÞ min μS2 ðQ S2 Þ; μR2 ðQ R2 Þ ∣ o 2ð1 k2 Þ 8 t Z T; and

ϵ

∣m21 ðP 1 ðtÞÞ m21 ðP n1 Þ∣ o ; 2

8 t Z T:

It then follows that ∣½ð1 k2 Þmin μS2 ðQ S2 ðtÞÞ; μR2 ðQ R2 ðtÞÞ m21 ðP 1 ðtÞÞ D ½ð1 k2 Þmin μS2 ðQ nS2 Þ; μR2 ðQ nR2 Þ m21 ðP n1 Þ D∣

J.P. Grover, F.-B. Wang / Journal of Theoretical Biology 351 (2014) 9–24

rð1 k2 Þ∣min μS2 ðQ S2 ðtÞÞ; μR2 ðQ R2 ðtÞÞ min μS2 ðQ nS2 Þ; μR2 ðQ nR2 Þ ∣ þ ∣m21 ðP 1 ðtÞÞ m21 ðP n1 Þ∣

ϵ ϵ

o þ ¼ ϵ; 2 2

8 t Z T:

Therefore, ð1 k2 Þ minfμS2 ðQ S2 ðtÞÞ; μR2 ðQ R2 ðtÞÞg m21 ðP 1 ðtÞÞ D

4ð1 k2 Þ minfμS2 ðQ nS2 Þ; μR2 ðQ nR2 Þg m21 ðP n1 Þ D ϵ;

8 t Z T:

ð29Þ min μS2 ðQ S2 Þ; μR2 ðQ R2 Þ By (20), it follows that ϵ n m21 ðP 1 Þ D 4 0. Substituting ϵ ¼ ϵ0 into (29), we get 1 0 ≔2½ð1 k2 Þ

n

ð1 k2 Þ minfμS2 ðQ S2 ðtÞÞ; μR2 ðQ R2 ðtÞÞg m21 ðP 1 ðtÞÞ D 4 ϵ0 ;

n

8 t Z T:

From the fourth equation of the system (12), it follows that dN 2 Z ϵ0 N 2 ; dt

8 t Z T:

Thus, limt-1 N2 ðtÞ ¼ 1 is a contradiction. If we suppose, by contradiction, that I ¼ fE2 g, it follows from (23) and the same arguments that limt-1 N1 ðtÞ ¼ 1 is a contradiction. If we suppose, by contradiction, that I ¼ fE0 g, it follows from (18), (21) and Remark 4.2 that limt-1 N1 ðtÞ ¼ 1, limt-1 N2 ðtÞ ¼ 1, which is a contradiction. Thus, I A=ffE0 g; fE1 g; fE2 gg. By using Zhao (2003, Theorem 1.3.1) with L ¼ I , it follows that there exists a δ 4 0 such that inf dðx; ∂Σ 0 Þ Z δ:

xAI

Since ðN1 ðtÞ; Q S1 ðtÞ; Q R1 ðtÞ; N 2 ðtÞ; Q S2 ðtÞ; Q R2 ðtÞ; P 1 ðtÞ; P 2 ðtÞ; W S ðtÞ; W R ðtÞÞ ~ ¼ I fð0; 0Þg as t-1; -ω it follows that there exists an η, such that lim inf N i ðtÞ Z η; t-1

i ¼ 1; 2:

This implies that the solution ﬂows Ψ~ t : Σ~ -Σ~ are uniformly persistent with respect to ðΣ~ 0 ; ∂Σ~ 0 Þ. By Zhao (2003, Theorem 1.3.6), it follows that system (12) admits at least one steady-state ^ S; W ^ R Þ A Σ~ 0 . It is easy solution ðN^ 1 ; Q^ S1 ; Q^ R1 ; N^ 2 ; Q^ S2 ; Q^ R2 ; P^ 1 ; P^ 2 ; W to see that ^ S; W ^ RÞ ðN^ 1 ; Q^ S1 ; Q^ R1 ; N^ 2 ; Q^ S2 ; Q^ R2 ; P^ 1 ; P^ 2 ; W ¼ ðN~ 1 ; Q~ S1 ; Q~ R1 ; N~ 2 ; Q~ S2 ; Q~ R2 ; P~ 1 ; P~ 2 ; 0; 0Þ; where ðN~ 1 ; Q~ S1 ; Q~ R1 ; N~ 2 ; Q~ S2 ; Q~ R2 ; P~ 1 ; P~ 2 Þ is deﬁned in (28). From (11), we deﬁne two positive number S~ and R~ which satisfy Sð0Þ S~ Q~ S1 N~ 1 Q~ S2 N~ 2 ¼ 0 and Rð0Þ R~ Q~ R1 N~ 1 Q~ R2 N~ 2 ¼ 0, ~ ~ ~ ~ ~ respectively. Indeed, ðS; R; N 1 ; Q S1 ; Q R1 ; N~ 2 ; Q~ S2 ; Q~ R2 ; P~ 1 ; P~ 2 Þ is a positive (coexistence) steady-state solution for the system (1). Since systems (1) and (12) are equivalent, the statement in (24) is true for the system (1). We complete our proof. □ References Berry, J.P., Gantar, M., Perez, M.H., Berry, G., Noriega, F.G., 2006. Cyanobacterial toxins as allelochemicals with potential applications as algaecides, herbicides and insecticides. Mar. Drugs 6, 117–146. Bhattacharyya, J., Smith, H.L., Pal, S., 2012. Allelopathy of plasmid-bearing and plasmid-free organisms competing for two complementary resources in a chemostat. J. Biol. Dyn. 6, 628–644. Braselton, R.J.P., Waltman, Paul, 2001. A competition model with dynamically allocated inhibitor production. Math. Biosci. 173, 55–84. Chakraborty, S., Roy, S., Chattopadhyay, J., 2008. Nutrient-limited toxin production and the dynamics of two phytoplankton in culture media: a mathematical model. Ecol. Model. 213, 191–201. Chesson, P., 2000. Mechanisms of maintenance of species diversity. Annu. Rev. Ecol. Syst. 31, 343–366. Cunningham, A., Nisbet, R.M., 1980. Time lag and co-operativity in the transient growth dynamics of microalgae. J. Theor. Biol. 84, 189–203.

23

Cunningham, A., Nisbet, R.M., 1983. Transient and oscillation in continuous culture. In: Bazin, M.J. (Ed.), Mathematics in Microbiology. Academic Press, New York Defreitas, M.J., Frederickson, A.G., 1978. Inhibition as a factor in the maintenance of the diversity of microbial ecosystems. J. Gen. Microbiol. 106, 307–320. Droop, M., 1968. Vitamin B12 and marine ecology. IV. The kinetics of uptake, growth and inhibition in Monochrysis Lutheri. J. Mar. Biol. Assoc. UK 48, 689–733. Droop, M., 1973. Some thoughts on nutrient limitation in algae. J. Phycol. 9, 264–272. Droop, M., 1974. The nutrient status of algal cells in continuous culture. J. Mar. Biol. Assoc. UK 54, 825–855. Fistarol, G.O., Legrand, C., Granéli, E., 2003. Allelopathic effect of Prymnesium parvum on a natural phytoplankton community. Mar. Ecol. Prog. Ser. 255, 115–125. Fu, F.X., Tatters, A.O., Hutchins, D.A., 2012. Global change and the future of harmful algal blooms in the ocean. Mar. Ecol. Prog. Ser. 470, 207–233. Granéli, E., Hansen, P.J., 2006. Allelopathy in harmful algae: a mechanism to compete for resources?. In: Granéli, E., Turner, J.T. (Eds.), Ecology of Harmful Algae. Springer-Verlag, Berlin, pp. 189–201 Granéli, E., Salomon, P.S., Fistarol, G.O., 2008. The role of allelopathy for harmful algae bloom formation. In: Evangelista, V., Barsanti, L., Frassanto, A.M., Passarelli, V., Gualtieri, P. (Eds.), Algal Toxins: Nature, Occurrence, Effect and Detection. Springer, Berlin, pp. 159–178 Grover, J.P., 1992. Constant- and variable-yield models of population growth: responses to environmental variability and implications for competition. J. Theor. Biol. 158, 409–428. Grover, J.P., 1997. Resource Competition. Chapman and Hall, London Grover, J.P., Crane, K.W., Baker, J.W., Brooks, B.W., Roelke, D.L., 2011. Spatial variation of harmful algae and their toxins in ﬂowing-water habitats: a theoretical exploration. J. Plankton Res. 33, 211–227. Grover, J.P., Wang, F.-B., 2013. Competition for one nutrient with internal storage and toxin mortality. Math. Biosci. 244, 82–90. Guildford, S.J., Hecky, R.E., 2000. Total nitrogen, total phosphorus, and nutrient limitation in lakes and oceans: is there a common relationship? Limnol. Oceanogr. 45, 1213–1223. Hallegraeff, G.M., 1993. A review of harmful algal blooms and their apparent global increase. Phycologia 32, 79–99. Hallegraeff, G.M., 2010. Ocean climate change, phytoplankton community responses, and harmful algal blooms: a formidable predictive challenge. J. Phycol. 46, 220–235. Hattenrath-Lehman, T.K., Gobler, C.J., 2011. Allelopathic inhibition of competing phytoplankton by North American strains of the toxic dinoﬂagellate Alexandrium fundyense: evidence from ﬁeld experiments, laboratory experiments, and bloom events. Harmful Algae 11, 106–116. Heil, C.A., Glibert, P.M., Al-Sarawi, M.A., Faraj, M., Behbehani, M., Husain, M., 2001. First record of a ﬁsh-killing Gymnodinium sp. bloom in Kuwait Bay, Arabian Sea: chronology and potential causes. Mar. Ecol. Prog. Ser. 214, 15–23. Heisler, J., Glibert, P.M., Burkholder, J.M., Anderson, D.M., Cochlan, W., Dennison, W. C., Dortch, Q., Gobler, C.J., Heil, C.A., Humphries, E., Lewitus, A., Magnien, R., Marshall, H.G., Sellner, K., Stockwell, D.A., Stoecker, D.K., Suddleson, J., 2008. Eutrophication and harmful algal blooms: a scientiﬁc consensus. Harmful Algae 8, 3–13. Hsu, S.B., Waltman, P., 2004. A survey of mathematical models of competition with an inhibitor. Math. Biosci. 187, 53–97. Jonsson, P.R., Pavia, H., Toth, G., 2009. Formation of harmful algal blooms cannot be explained by allelopathic interaction, Proc. Natl. Acad. Sci. USA 106, 11177–11182. Klausmeier, C.A., Litchman, E., Levin, S.A., 2004. Phytoplankton growth and stoichiometry under multiple nutrient limitation. Limnol. Oceanogr. 49, 1463–1470. Legrand, C., Rengefors, K., Fistarol, G.O., Granéli, E., 2003. Allelopathy in phytoplankton-biochemical, ecological and evolutionary aspects. Phycologia 42, 409–419. León, J.A., Tumpson, D.B., 1975. Competition between two species for two complementary or substitutable resources. J. Theor. Biol. 50, 185–201. Li, B., Smith, H.L., 2007. Global dynamics of microbial competition for two resources with internal storage. J. Math. Biol. 55 (4), 481–515. Litchman, E., Klausmeier, C.A., Miller, J.R., Schoﬁeld, O.M., Falkowski, P.G., 2006. Multi-nutrient, multi-group model of present and future oceanic phytoplankton communities. Biogeosciences 3, 585–606. Martines, I.P., Kojouharov, H.V., Grover, J.P., 2009. A chemostat model of resource competition and allelopathy. Appl. Math. Comput. 215, 573–583. McPeek, M.A., 2012. Intraspeciﬁc density dependence and a guild of consumers coexisting on one resource. Ecology 93, 2728–2735. Morel, F.M.M., 1987. Kinetics of nutrient uptake and growth in phytoplankton. J. Phycol. 23, 137–150. Murata, M., Yasumoto, T., 2000. The structure elucidation and biological activities of high molecular weight algal toxins: maitotoxin, prymnesins and zooxanthellatoxins. Nat. Prod. Rep. 17, 293–314. Olli, K., Trunov, K., 2007. Self-toxicity of Prymnesium parvum (Prymnesiophyceae). Phycologia 46, 109–112. Remmel, E.J., Kohmescher, N., Larson, J.H., Hambright, K.D., 2011. An experimental analysis of harmful algae-zooplankton interactions and the ultimate defense. Limnol. Oceanogr. 56, 461–470. Rhee, G.-Y., 1974. Phosphate uptake under nitrate limitation by Scenedesmus sp. and its ecological implications. J. Phycol. 10, 470–475.

24

J.P. Grover, F.-B. Wang / Journal of Theoretical Biology 351 (2014) 9–24

Rhee, G.-Y., 1978. Effects of N: P atomic ratios and nitrate limitation on algal growth, cell composition, and nitrate uptake. Limnol. Oceanogr. 27, 1101–1112. Smayda, T.J., 1990. Novel and nuisance phytoplankton blooms in the sea: evidence for a global epidemic. In: Granéli, E., Sandstrom, B., Edler, L., Anderson, D.M. (Eds.), Toxic Marine Phytoplankton. Elsevier Science Publishing Co., New York, pp. 29–40 Smith, H.L., 1995. Monotone dynamical systems: an introduction to the theory of competitive and cooperative systems. In: Mathematical Surveys and Monographs, vol. 41. American Mathematical Society Providence, RI. Smith, H.L., Thieme, H.R., 2011. Dynamical Systems and Population Persistence, Graduate Studies in Mathematics, vol. 118. American Mathematical Society. Smith, V.H., Bennett, S.J., 1999. Nitrogen: phosphorus supply ratios and phytoplankton community structure in lakes. Arch. Hydrobiol. 146, 37–53. Smith, V.H., Tilman, G.D., Nekola, J.C., 1999. Eutrophication: impacts of excess nutrient inputs on freshwater, marine, and terrestrial ecosystems. Environ. Pollut. 100, 179–196. Sommer, U., 1983. Nutrient competition between phytoplankton species in multispecies chemostat experiments. Arch. Hydrobiol. 96, 399–416. Sommer, U., 1989. The role of competition for resources in phytoplankton succession. In: Sommer, U. (Ed.), Plankton Ecology: Succession in Plankton Communities. Springer-Verlag, Berlin, pp. 57–106 Sommer, U., 1990. Phytoplankton nutrient competition-from laboratory to lake. In: Grace, J., Tilman, D. (Eds.), Perspectives on Plant Competition. Academic Press, New York, pp. 193–213 Sommer, U., 2002. Competition and coexistence in plankton communities. In: Sommer, B., Worm, B. (Eds.), Competition and Coexistence. Springer-Verlag, Berlin, pp. 79–108

Sopanen, S., Koski, M., Kuuppo, P., Uronen, P., Legrand, C., Tamminen, T., 2006. Toxic haptophyte Prymnesium parvum affects grazing, survival, egestion and egg production of the calanoid copepods Eurytemora afﬁnis and Acartia biﬁlosa. Mar. Ecol. Prog. Ser. 327, 223–232. Southard, G.M., Fries, L.T., Barkoh, A., 2010. Prymnesium parvum: the Texas experience. J. Am. Water Resour. Assoc. 46, 14–23. Sterner, R.W., 2008. On the phosphorus limitation paradigm for lakes. Int. Rev. Hydrobiol. 93, 433–445. Thieme, H.R., 1992. Convergence results and a Poincare–Bendixson trichotomy for asymptotically autonomous differential equations. J. Math. Biol., 755–763 Tillman, U., 2003. Kill and eat your predator: a winning strategy of the planktonic ﬂagellate Prymnesium parvum. Aquat. Microb. Ecol. 32, 73–84. Tilman, D., 1981. Tests of resource competition theory using four species of Lake Michigan algae. Ecology 62, 802–815. Tilman, D., 1982. Resource Competition and Community Structure. Princeton University Press, Princeton Tilman, D., Kilham, S.S., Kilham, P., 1982. Phytoplankton community ecology: the role of limiting nutrients. Annu. Rev. Ecol. Syst. 13, 349–372. Turpin, D.H., 1988. Physiological mechanisms in phytoplankton resource competition. In: Sandgren, C.D. (Ed.), Growth and Reproductive Strategies of Freshwater Phytoplankton. Cambridge University Press, Cambridge, pp. 316–368 Wang, J., Zhang, Y., Li, H., Cao, J., 2013. Competitive interaction between diatom Skeletonema costatum and dinoﬂagellate Prorocentrum donghaiense in laboratory culture. J. Plankton Res. 35, 367–378. Zhai, C., Song, S., Zou, S., Liu, C., Xue, Y., 2013. The mechanism of competition between two bloom-forming Microcystis species. Freshw. Biol. 58, 1831–1839. Zhao, X.-Q., 2003. Dynamical Systems in Population Biology. Springer, New York