Research article

Coexistence and extinction in a data-based ratio-dependent model of an insect community

  • Received: 24 February 2020 Accepted: 23 April 2020 Published: 27 April 2020
  • In theory, pure competition often leads to competitive exclusion of species. However, what we often see in nature is a large number of distinct predator or consumer species coexist in a community consisting a smaller number of prey or plant species. In an effort of dissecting how indirect competition and selective predation may have contributed to the coexistence of species in an insect community, according to the replicated cage experiments (two aphid species and a specialist parasitoid that attacks only one of the aphids) and proposed mathematical models, van Veen et. al. [5] conclude that the coexistence of the three species is due to a combination of density-mediated and trait-mediated indirect interactions. In this paper, we formulate an alternative model that observes the conventional law of mass conservation and provides a better fitting to their experimental data sets. Moreover, we present an initial attempt in studying the stabilities of the nonnegative steady states of this model.

    Citation: Yang Kuang, Kaifa Wang. Coexistence and extinction in a data-based ratio-dependent model of an insect community[J]. Mathematical Biosciences and Engineering, 2020, 17(4): 3274-3293. doi: 10.3934/mbe.2020187

    Related Papers:

    [1] Juan Ye, Yi Wang, Zhan Jin, Chuanjun Dai, Min Zhao . Dynamics of a predator-prey model with strong Allee effect and nonconstant mortality rate. Mathematical Biosciences and Engineering, 2022, 19(4): 3402-3426. doi: 10.3934/mbe.2022157
    [2] Eric M. Takyi, Charles Ohanian, Margaret Cathcart, Nihal Kumar . Dynamical analysis of a predator-prey system with prey vigilance and hunting cooperation in predators. Mathematical Biosciences and Engineering, 2024, 21(2): 2768-2786. doi: 10.3934/mbe.2024123
    [3] Eric Ruggieri, Sebastian J. Schreiber . The Dynamics of the Schoener-Polis-Holt model of Intra-Guild Predation. Mathematical Biosciences and Engineering, 2005, 2(2): 279-288. doi: 10.3934/mbe.2005.2.279
    [4] Yuanfu Shao . Bifurcations of a delayed predator-prey system with fear, refuge for prey and additional food for predator. Mathematical Biosciences and Engineering, 2023, 20(4): 7429-7452. doi: 10.3934/mbe.2023322
    [5] Kawkab Al Amri, Qamar J. A Khan, David Greenhalgh . Combined impact of fear and Allee effect in predator-prey interaction models on their growth. Mathematical Biosciences and Engineering, 2024, 21(10): 7211-7252. doi: 10.3934/mbe.2024319
    [6] Yun Kang, Sourav Kumar Sasmal, Amiya Ranjan Bhowmick, Joydev Chattopadhyay . Dynamics of a predator-prey system with prey subject to Allee effects and disease. Mathematical Biosciences and Engineering, 2014, 11(4): 877-918. doi: 10.3934/mbe.2014.11.877
    [7] Saheb Pal, Nikhil Pal, Sudip Samanta, Joydev Chattopadhyay . Fear effect in prey and hunting cooperation among predators in a Leslie-Gower model. Mathematical Biosciences and Engineering, 2019, 16(5): 5146-5179. doi: 10.3934/mbe.2019258
    [8] Sourav Kumar Sasmal, Jeet Banerjee, Yasuhiro Takeuchi . Dynamics and spatio-temporal patterns in a prey–predator system with aposematic prey. Mathematical Biosciences and Engineering, 2019, 16(5): 3864-3884. doi: 10.3934/mbe.2019191
    [9] Lazarus Kalvein Beay, Agus Suryanto, Isnani Darti, Trisilowati . Hopf bifurcation and stability analysis of the Rosenzweig-MacArthur predator-prey model with stage-structure in prey. Mathematical Biosciences and Engineering, 2020, 17(4): 4080-4097. doi: 10.3934/mbe.2020226
    [10] Mengya Huang, Anji Yang, Sanling Yuan, Tonghua Zhang . Stochastic sensitivity analysis and feedback control of noise-induced transitions in a predator-prey model with anti-predator behavior. Mathematical Biosciences and Engineering, 2023, 20(2): 4219-4242. doi: 10.3934/mbe.2023197
  • In theory, pure competition often leads to competitive exclusion of species. However, what we often see in nature is a large number of distinct predator or consumer species coexist in a community consisting a smaller number of prey or plant species. In an effort of dissecting how indirect competition and selective predation may have contributed to the coexistence of species in an insect community, according to the replicated cage experiments (two aphid species and a specialist parasitoid that attacks only one of the aphids) and proposed mathematical models, van Veen et. al. [5] conclude that the coexistence of the three species is due to a combination of density-mediated and trait-mediated indirect interactions. In this paper, we formulate an alternative model that observes the conventional law of mass conservation and provides a better fitting to their experimental data sets. Moreover, we present an initial attempt in studying the stabilities of the nonnegative steady states of this model.



    Understanding the rich diversity or coexistence of populations observed in the nature is a main goal in many ecological research activities. The major interactions among community of populations are competition and predation which were the focuses of the classical book on mathematical models of population ecology authored by Freedman in 1980 [1]. In a theoretical setting where populations competing for a single resource, only one population can persist. This theoretical finding is often termed as the competitive exclusion principle [2]. Indeed, Gause observed that two species competing for the same limiting resource cannot coexist at constant population values [3]. Pure predator-prey interactions give rise to predator-prey communities called simple food chains. A simple food chain often consists of only a few species in ecological settings which prevents it to be the framework to host many species. However, predator-prey interaction often results in oscillatory population densities which may provide variable resources levels to enable the coexistence of many competing species [4]. Therefore coexistence of many species is mostly likely the consequence of concurrent competition and predation interactions among species.

    In a remarkable effort of dissecting how density-mediated (indirect competition) and trait-mediated (selective predation) interactions among species may have contributed to the coexistence of species, van Veen et. al. studied an insect community consisting of two aphid species (Acyrthosiphon pisum and Megoura viciae) and a specialist parasitoid (Aphidius ervi) that attacks only one of the aphids (A. pisum). In extensive experiments, they found that the two aphid species alone were unable to coexist, with A. pisum competitively excluding M. viciae. Moreover, they observed that the interaction between A. pisum and the parasitoid was unstable. However, the three-species community persisted for at least 50 weeks [5]. It is observed that parasitoid attack on the susceptible host reduces the interspecific competition experienced by the non-host (a density-mediated effect), and the presence of the non-host reduces the searching efficiency of the parasitoid (a trait-mediated effect). To study this experiment analytically, van Veen et. al. [5] proposed the following three-species interaction dynamic model:

    {dN1dt=r1N1(1α11N1α12N2)N1α1pP1+b1N1+b2N2,dN2dt=r2N2(1α21N1α22N2),dPdt=N1sα1pP1+b1N1+b2N2+cPμP, (1.1)

    where N1 and N2 are numbers of the two aphid species, A. pisum and M. viciae, respectively. P denotes the number of the specialist parasitoid, A. ervi, which attacks only N1. The parameter ri is the per capita growth rate (or the intrinsic rate of increase) of aphid species i, i=1,2. The parameter αii is the intraspecific competition coefficient (also be interpreted as the reciprocal of carrying capacity) and αij,ij, is the interspecific competition coefficient, i,j=1,2. 1/μ represents the average lifespan of the parasitoid. Parameter α1p is the parasitoid per capita attack rate of N1, while b1 is the parameter controlling the reduction in parasitoid attack rate with increasing host number (N1), which may be interpreted as correlated with the searching and handling time, and b2 can be interpreted as the time wasted when a parasitoid encounters an unsuitable aphid (N2), which is equivalent to the "recognition time" in classical diet models from foraging theory. s is the parasitoid sex ratio, and c denotes the effect of parasitoid number on parasitoid recruitment. All these parameters are positive constants.

    Note that the functional response and incidence functions of parasitoid are not proportional in model (1.1), i.e., they don't obey the usual conservation law of mass. While proportionality is not a must in reality, the lack of it implies the conservation of some weighted sum of masses is lost which may artificially add complexity to an already complicated population growth process. In addition, Figure 1 shows that there is a big discrepancy between the fitted population trajectories of (1.1) and the experimental data sets in [5]. In particular, parasitoid numbers are sometimes much larger. In fact, from the full three-species cage experiments in [5], people can infer that parasitoid (A. ervi) have to search for the susceptible host (A. pisum), and the other aphids (M. viciae) may interfere the search of parasitoid at the same time. In an effort to reduce this discrepancy, we employ a standard incidence function to represent both the functional response and incidence functions of parasitoid. This yields the following modified model:

    {dN1dt=r1N1(1α11N1α12N2)N1α1pPb1N1+b2N2+cP,dN2dt=r2N2(1α21N1α22N2),dPdt=N1sα1pPb1N1+b2N2+cPμP. (1.2)
    Figure 1.  Fit of model (1.1) and (1.2) to the experimental data sets for A. pisum, M. viciae and A. ervi in the full-community experiment in [5]. Since r1,r2,α11,α22 and μ are species-specific growth parameters, we adopt the estimated values in [5]. Specifically, r1=3.22/week,r2=2.82/week,α11=3.82×104,α22=3.82×104. For other parameters, according to [5], we take α12=3.70×104,α21=3.97×104,α1p=2.81×101/week,μ=6.34×101/week,b1=2.33×102,b2=4.34×102,c=1.26 and s=0.50 in (1.1). Minimizing the squared deviations of the model and experimental data, we take α12=2.39×104,α21=5.7116×104,α1p=3.5190×102/week,μ=0.50/week,b1=5.2925×105,b2=6.4178×105,c=1.7527×102 and s=7.9556×103 in (1.2).

    From Table 1, we observe that the model (1.2) provides a better fitting than that of model (1.1) in Figure 1. Most noteworthy is the sharp reduction of the aforementioned big discrepancy between the data and the model (1.2) fitting of the parasitoid numbers compare to that of the data and the model (1.1) fitting of the parasitoid. Thus, (1.2) may be more suitable to describe the full three-species cage experiments in [5]. In the present paper, we will study the dynamics of model (1.2) with a focus on the local and global stabilities of its nonnegative steady states. Observe that the function

    R(N1,N1,P)α1pPb1N1+b2N2+cP
    Table 1.  Square sum of errors (SSE) of model fitting to the experimental data sets in [5] for model (1.1) and (1.2).
    Model (1.1) Model (1.2)
    SSEN1=1.0467×106 SSEN1=1.6146×106
    SSEN2=4.8173×106 SSEN2=3.4856×106
    SSEP=1.3738×105 SSEP=7.7021×102
    SSETotal=6.0014×106 SSETotal=5.1010×106

     | Show Table
    DownLoad: CSV

    is not defined at (0,0,0). However, lim(N1,N1,P)(0+,0+,0+)R(N1,N1,P)=α1p/c. Hence, in the following, we define R(0,0,0)=α1p/c. With this definition, we see that (0,0,0) is a steady state of model (1.2). However, the model is not differentiable at (0,0,0). Indeed, model (1.2) is a ratio-dependent model with two preys and one specialist predator. A ratio-dependent model with one prey and two predators was studied in [6] while a ratio-dependent food chain was studied in [7]. In these studies, through a nonlinear transformation (blow-up transformation), the authors were able to reveal the rich dynamics often observed in ratio-dependent models due to the non-smoothness of the trivial steady state (0,0,0). Indeed, rich dynamics, such as global stability, limit cycles and extinction, was also found in ratio-dependent models with only two populations [8,9,10,11].

    Since we are interested in the long term dynamics of the population interactions, we assume that the initial condition of (1.2) has the form

    N1(0)>0,N2(0)>0andP(0)>0. (2.1)

    We also assume that parameters of (1.2) are all positive. It is easy to see that the right hand side functions of the equations of (1.2) are uniformly Lipschitzian and hence the solution of model (1.2) with initial conditions (2.1) is unique and exist for all t>0 [12]. We first show that that model (1.2) produces solutions that are biologically plausible. Mathematically, this is equivalent to establish the following proposition.

    Proposition 2.1. For model (1.2), solutions with initial conditions (2.1) are positive and eventually uniformly bounded.

    Before proving the above proposition, we would like to establish the following general result.

    Lemma 2.1. Assume the function F(t,x) is uniformly Lipschitzian with respect to x and there is a continuous function C(t) such that F(t,x)C(t)x for x0. Then the solution of

    dxdt=F(t,x),x(0)>0

    exists, is unique and positive for all t>0.

    Proof. The existence and uniqueness follow from the uniform Lipschitzian property of the function F(t,x) and Theorem 5.2 in [12]. Since F(t,x)C(t)x for x0 we see that dxdtC(t)x which implies that for all t>0, we have x(t)x(0)et0C(s)ds>0.

    We are now in a position to prove the Proposition 2.1.

    Proof. We first establish the positiviness of P(t). Notice that dPdt=FP(N1,N2,P)=PfP(N1,N2,P) where fP(N1,N2,P)μ for positive values of N1(t),N2(t) and P(t). By Lemma 2.1 with C(t)=μ, we see that P(t) stays positive.

    We now consider the positivity of N1(t). Observe that

    dN1dt=FN1(N1,N2,P)=N1fN1(N1,N2,P)C1(t)N1(t),

    where

    fN1=r1(1α11N1α12N2)α1pPb1N1+b2N2+cP.

    Note that α1pPb1N1+b2N2+cPα1pc is always valid because all variables are nonnegative. As a result, we have

    C1(t)=min0st{r1[α11N1(s)+α12N2(s)]α1pc.}

    By Lemma 2.1 with C(t)=C1(t), we see that N1(t) stays positive. The positivity of N2(t) can be established similarly.

    Now we present the arguments for ultimate boundedness of solutions of (1.2) in R3+. Since all components of a solution of (1.2) are positive, using the first and second equations of (1.2), we have

    dN1dt<r1N1(1α11N1)anddN2dt<r2N2(1α22N2).

    As a result, we have lim supt+N1(t)1/α11N1 and lim supt+N2(t)1/α22N2, respectively.

    Using the monotonicity of sα1pN1Pb1N1+cP on N1 and above analysis, we have

    sα1pN1Pb1N1+b2N2+cP<sα1pN1Pb1N1+cPsα1pN1Pb1N1+cP.

    Thus, using the third equation of (1.2), we have

    dPdt<(sα1pN1b1N1+cPμ)P=[(sα1pb1μ)N1cμP]Pb1N1+cP.

    As a result, we have lim supt+P(t)sα1pb1μcμP if sα1pb1μ>0, and limt+P(t)=0 if sα1pb1μ0.

    A hallmark of the rich dynamics of ratio-dependent population models is the possibility of the origin as an attractor, implying the collapse of the population community [7,10,11,13]. Contrast to many ratio-dependent population models, it is easy to see from the N2 equation that the origin as an equilibrium of model (1.2) can not be an attractor. This maybe intuitive biologically since the parasitoid species only attack one of the species which may indirectly help the persistence of the other aphid species. In addition, we have the following much stronger result.

    Proposition 2.2. Assume that α1p<cr1 or α1pb1μ/s. If (N1(t),N2(t),P(t)) is a solution of model (1.2) with initial conditions (2.1), then

    lim inft+(N1+N2)>0. (2.2)

    Proof. We have two cases to consider: 1) α1p<cr1 and 2) α1pb1μ/s. In the following, we let Z(t)N1(t)+N2(t).

    Consider first case 1. We prove it by contradiction. If the proposition is false, then there is a strictly increasing sequence of positive values ti such that Z(ti)<1/i is a strictly decreasing sequence and Z(ti)0, i=1,2,.... Since all components of the solution stay positive, we have

    α1pPb1N1+b2N2+cP<α1pc.

    Hence

    dZdt>r1N1(1α1pcr1α11N1α12N2)+r2N2(1α21N1α22N2).

    Since α1pcr1, we see for large enough values of i, we will have both 1α1pcr1α11N1α12N2>0 and 1α21N1α22N2>0. Hence dZ(ti)dt>0 which contradicts the fact that ti is selected to have Z(ti)0.

    Consider now the case 2 with α1pb1μ/s. From the proof of the case 1, the conclusion of the proposition is true if α1p<cr1. Hence we assume below that cr1α1pb1μ/s. The proof of Proposition 2.1 shows that α1p<b1μ/s implies limt+P(t)=0. Hence for any δ>0 there is a time tδ>0 such that t>tδ implies that P(t)<δ. Let

    0<δ0<1(α11+α12+α21+α22),

    and δ be so small such that

    0<α1pδr1(min{b1,b2}δ0+cδ)<1(α11+α12)δ0. (2.3)

    We claim that if for some t0>tδ such that Z(t0)δ0, then Z(t)δ0 for all t>t0. If this claim is false, then there is a t1to such that Z(t1)=δ0 and Z(t1)0. Observe that

    dZdt=r1N1(1α1pPr1(b1N1+b2N2+cP)α11N1α12N2)+r2N2(1α21N1α22N2). (2.4)

    However,

    dZdt(t1)r1N1(t1)(1α1pδr1(min{b1,b2}δ0+cδ(α11+α12)δ0)+r2N2(t1)(1(α21+α22)δ0)>0. (2.5)

    This is a contradiction which proves our claim. Clearly this claim implies the conclusion of the proposition.

    If the statement that Z(t0)δ0 for some t0>tδ is not true, then Z(t)<δ0 for all ttδ. In this situation, we have

    dN2dt=r2N2(1α21N1α22N2)>r2(1(α21+α22)δ0)N2.

    This implies that for all ttδ, we have Z(t)>N2(t)>N2(tδ)egt where g=1(α21+α22)δ0>0. Therefore Z(t) will be unbounded which contradicts the statement that Z(t)<δ0 for all ttδ. The proof of the proposition is now complete.

    The proof of the second case for the Proposition 2.2 implies that if limt+P(t)=0, then the Proposition 2.2 hold. Hence we have the following result.

    Theorem 2.1. If (N1(t),N2(t),P(t)) is a solution of model (1.2) with initial conditions (2.1), then

    lim inft+(N1+N2)+lim supt+P>0. (2.6)

    We conjecture that Proposition 2.2 remains true even if the condition α1p<cr1 or α1pb1μ/s is removed.

    Since sα1pb1μ0 implies limt+P(t)=0, this together with Proposition 2.2 imply that if α1p<b1μ/s in system (1.2), then its limit system is

    {dN1dt=r1N1(1α11N1α12N2),dN2dt=r2N2(1α21N1α22N2). (2.7)

    System (2.7) is the well studied Lotka-Volterra competition model (see related sections in [1], [12] or [14]). Let D1=α22α12,D2=α11α21 and D=α11α22α12α21. We have the following global results for system (1.2).

    Proposition 2.3. Suppose that α1p<b1μ/s in system (1.2), then limt+P(t)=0 and the following statements are true.

    (i) When α11<α21 and α22>α12, equilibrium E1=(1α11,0,0) is globally asymptotically stable and E2=(0,1α22,0) is unstable.

    (ii) When α11>α21 and α22<α12, system (1.2) has no coexist equilibrium, equilibrium E1 is always unstable and E2 is globally asymptotically stable.

    (iii) When α11>α21 and α22>α12, system (1.2) has a unique coexist equilibrium E12=(D1D,D2D,0), which is globally asymptotically stable, both equilibria E1 and E2 are unstable.

    (iv) When α11<α21 and α22<α12, the unique coexist equilibrium E12 is a saddle, both equilibria E1 and E2 are stable, i.e., the bistability phenomenon of initial value dependence will occur.

    In the rest of this paper, we assume that sα1pb1μ>0 is valid in system (1.2).

    In this section, we study the existence, location and numbers of equilibrium of system (1.2). Let

    {r1N1(1α11N1α12N2)N1α1pPb1N1+b2N2+cP=0,r2N2(1α21N1α22N2)=0,N1sα1pPb1N1+b2N2+cPμP=0. (3.1)

    It is easy to see that system (1.2) has boundary equilibria E0=(0,0,0),E1=(1/α11,0,0) and E2=(0,1/α22,0). Recall that D1=α22α12,D2=α11α21 and D=α11α22α12α21. After straightforward algebraic manipulation on (3.1), we obtain the following results.

    Proposition 3.1. The following are true for system (1.2).

    (i) When D1,D2 and D have the same sign, i.e., α11>α21 and α22>α12, or α11<α21 and α22<α12, equilibrium E12=(D1D,D2D,0) exists.

    (ii) When r1s>sα1pb1μc, equilibrium E13=(˜N1,0,˜P) exists, and

    (˜N1,0,˜P)=(r1sc(sα1pb1μ)r1scα11,0,sα1pb1μcμ˜N1).

    Now, we consider the existence of positive equilibrium of system (1.2). According to (3.1), we know that the positive equilibrium must be the solution of the following algebraic equations

    {r1(1α11N1α12N2)α1pPb1N1+b2N2+cP=0,N2=1α22α21α22N1ϕ(N1),α1pb1N1+b2N2+cP=μsN1. (3.2)

    For (3.2), substituting the second and third equations into the first equation, and substituting the second equation into the third equation, we have

    {μP=r1sN1(θ1θN1)φ1(N1),μP=ω1N1ω2φ2(N1), (3.3)

    where

    θ1=D1α22,θ=Dα22,ω1=sα1pb1μc+b2μα21cα22,andω2=b2μcα22.

    According φ1(N1)=φ2(N1), we have

    φ(N1)r1sθN21+(ω1r1sθ1)N1ω2=0. (3.4)

    Note that =(r1sθ1ω1)2+4r1sθω2 is the discriminant of (3.4). According to the formula of roots of a quadric equation, we know the roots of (3.4) are

    N11=r1sθ1ω1+2r1sθ,N21=r1sθ1ω12r1sθ. (3.5)

    Furthermore, we have

    φ(ω2ω1)=r1sω2(θω2θ1ω1)ω21,φ(θ1θ)=θ1ω1θω2θ.

    Combining (3.3), (3.4) and (3.5), similar to Proposition 2.3, we split the analysis into four subsections.

    In this case, we know D1>0, i.e., θ1>0, but the sign of θ is undetermined.

    When θ>0, we know >|r1sθ1ω1| is always valid. So N11>0 and N21<0 in (3.5), i.e., N11 is the unique positive root of (3.4). Note that θ1/θ>1/α21 due to α11<α21. We know ω2ω1<N11<1α21 can ensure ϕ(N11)>0,φ1(N11)>0 and φ2(N11)>0.

    When θ<0, we will analyze the existence of positive equilibrium of system (1.2) according to the sign of and r1sθ1ω1. If <0 or r1sθ1ω1>0, it is clear that there is no positive equilibrium in system (1.2). If =0 and r1sθ1ω1<0, we have N11=N21=r1sθ1ω12r1sθN121>0. To ensure ϕ(N121)>0,φ1(N121)>0 and φ2(N121)>0, we also need ω2ω1<N121<1α21. If >0 and r1sθ1ω1<0, we know N11>0,N21>0 in (3.5). Similarly, we know ϕ(Nk1)>0,φ1(Nk1)>0 and φ2(Nk1)>0 if ω2ω1<Nk1<1α21,k=1,2.

    Summarize the above analyses, we have

    Proposition 3.2. Suppose that sα1pb1μ>0,α11<α21 and α22>α12.

    (i) If θ>0 and ω2ω1<N11<1α21, then system (1.2) has a unique positive equilibrium E+=(N11,N2,P), in which N2=ϕ(N11),P=φ1(N11)μ. Otherwise, there is no positive equilibrium in system (1.2).

    (ii) If θ<0,=0 and r1sθ1ω1<0, then system (1.2) has a unique positive equilibrium E+=(N121,N2,P) when ω2ω1<N121<1α21, in which N2=ϕ(N121),P=φ1(N121)μ. If θ<0,>0 and r1sθ1ω1<0, system (1.2) may have two positive equilibria Ek+=(Nk1,Nk2,Pk), which depends on whether ω2ω1<Nk1<1α21 is satisfied, here Nk2=ϕ(Nk1),Pk=φ1(Nk1)μ,k=1,2. Otherwise, there is no positive equilibrium in system (1.2).

    In this case, we know D1<0, i.e., θ1<0, but the sign of θ is also undetermined.

    When θ>0, similar to the above analyses, we know N11 is the unique positive root of (3.4). However, because φ1(x)<0 is always valid for x>0 if θ1<0 and θ>0, there is no positive equilibrium in system (1.2).

    When θ<0, since r1sθ1ω1<0 is always valid, we will analyze the existence of positive equilibrium of system (1.2) according to the sign of . If <0, it is clear that there is no positive equilibrium in system (1.2). If =0, we have N11=N21=r1sθ1ω12r1sθN121>0. Note that θ1/θ>1/α21 due to α11>α21 and the sign of θ,θ1. We know ω2ω1<N121<1α21 can ensure ϕ(N121)>0,φ1(N121)>0 and φ2(N121)>0. If >0, we know N11>0,N21>0 in (3.5). Similarly, we know ϕ(Nk1)>0,φ1(Nk1)>0 and φ2(Nk1)>0 if ω2ω1<Nk1<1α21,k=1,2.

    Summarize the above analyses, we have

    Proposition 3.3. Suppose that sα1pb1μ>0,α11>α21 and α22<α12.

    (i) If θ>0, then there is no positive equilibrium in system (1.2).

    (ii) If θ<0 and =0, then system (1.2) has a unique positive equilibrium E+=(N121,N2,P) when ω2ω1<N121<1α21, in which N2=ϕ(N121),P=φ1(N121)μ. If θ<0 and >0, system (1.2) may have two positive equilibria Ek+=(Nk1,Nk2,Pk), which depends on whether ω2ω1<Nk1<1α21 is satisfied, here Nk2=ϕ(Nk1),Pk=φ1(Nk1)μ,k=1,2. Otherwise, there is no positive equilibrium in system (1.2).

    In this case, we know D1>0,D>0, i.e., θ1>0,θ>0. Thus, >|r1sθ1ω1| is always valid, and N11 is the unique positive root of (3.4). Note that θ1/θ<1/α21 in this case. We know ω2ω1<N11<θ1θ can ensure ϕ(N11)>0,φ1(N11)>0 and φ2(N11)>0. As a result, we have

    Proposition 3.4. Suppose that sα1pb1μ>0,α11>α21 and α22>α12. If ω2ω1<N11<θ1θ, then system (1.2) has a unique positive equilibrium E+=(N11,N2,P), in which N2=ϕ(N11),P=φ1(N11)μ. Otherwise, there is no positive equilibrium in system (1.2).

    In this case, we know D1<0,D<0, i.e., θ1<0,θ<0. Since r1sθ1ω1<0 is always true, we will analyze the existence of positive equilibrium of system (1.2) according to the sign of . If <0, it is clear that there is no positive equilibrium in system (1.2). If =0, we have N11=N21=r1sθ1ω12r1sθN121>0. Note that θ1/θ<1/α21 in this case. We know ω2ω1<N121<θ1θ can ensure ϕ(N121)>0,φ1(N121)>0 and φ2(N121)>0. If >0, we have N11>0,N21>0 in (3.5). Similarly, we know ϕ(Nk1)>0,φ1(Nk1)>0 and φ2(Nk1)>0 if ω2ω1<Nk1<θ1θ,k=1,2.

    Summarize the above analyses, we have

    Proposition 3.5. Suppose that sα1pb1μ>0,α11<α21 and α22<α12. If =0, then system (1.2) has a unique positive equilibrium E+=(N121,N2,P) when ω2ω1<N121<1α21, in which N2=ϕ(N121),P=φ1(N121)μ. If >0, then system (1.2) may have two positive equilibria Ek+=(Nk1,Nk2,Pk), which depends on whether ω2ω1<Nk1<1α21 is satisfied, here Nk2=ϕ(Nk1),Pk=φ1(Nk1)μ,k=1,2. Otherwise, there is no positive equilibrium in system (1.2).

    In order to gain a global understanding of the rich and complex dynamics of (1.2), we would like to study the stability of all its equilibria. Note that (1.2) is not differentiable at origin equilibrium E0, and the standard linearization approach cannot be used at E0. To avoid this difficulty for now, this section will focus on the non-origin equilibria, including boundary equilibria E1,E2,E12,E13 in Proposition 3.1, and the possible positive equilibria E+,Ek+,k=1,2.

    The Jacobian matrix of (1.2) at E1 and E2 are

    J(E1)=[r1r1α12α11α1pb10r2(α11α21)α11000sα1pb1μb1] (4.1)

    and

    J(E2)=[r1(α22α12)α2200r2α21α22r2000μ], (4.2)

    respectively. Note that the eigenvalues of (4.1) and (4.2) lie on the diagonal. Hence, we have the following results.

    Proposition 4.1. Suppose that sα1pb1μ>0 in (1.2).

    (i) Equilibrium E1 is a saddle.

    (ii) Equilibrium E2 is locally asymptotically stable if α22<α12.

    The Jacobian matrix of (1.2) at E12 is

    J(E12)=[r1α11D1Dr1α12D1Dα1pD1b1D1+b2D2r2α21D2Dr2α22D2D000(sα1pb1μ)D1b2μD2b1D1+b2D2]. (4.3)

    Clearly, the element in the lower right-hand corner of (4.3) is one of eigenvalue of (4.3). When sα1pb1μ>0, (sα1pb1μ)D1b2μD2b1D1+b2D2<0 is equivalent to sα1pb1μb2μ<D2D1. Furthermore, recall D=α11α22α12α21, we can conclude that the trace of the upper left-hand 2×2 matrix is always negative and its determinant is positive if all D1,D2,D are positive. Thus, we have

    Proposition 4.2. Suppose that sα1pb1μ>0 in (1.2). Equilibrium E12 is locally asymptotically stable if D1>0,D2>0,D>0 and sα1pb1μb2μ<D2D1.

    The Jacobian matrix of (1.2) at E13 is

    J(E13)=[B1B2B30r2(1α21˜N1)0B4B5B6]. (4.4)

    Here

    B1=r1(12α11˜N1)cα1p˜P2(b1˜N1+c˜P)2,B2=r1α12˜N1+b2α1p˜N1˜P(b1˜N1+c˜P)2,B3=b1α1p˜N21(b1˜N1+c˜P)2,B4=scα1p˜P2(b1˜N1+c˜P)2,B5=sα1pb2˜N1˜P(b1˜N1+c˜P)2,B6=sα1pb1˜N21(b1˜N1+c˜P)2μ.

    Clearly, λ1=r2(1α21˜N1) is one of eigenvalue of (4.4), and the other two eigenvalues of (4.4) satisfying

    H(λ)λ2+A1λ+A2=0. (4.5)

    Here

    A1=μsα1pb1˜N21(b1˜N1+c˜P)2+cα1p˜P2(b1˜N1+c˜P)2r1(12α11˜N1),A2=(sα1pb1˜N21(b1˜N1+c˜P)2μ)(cα1p˜P2(b1˜N1+c˜P)2+r1(12α11˜N1))+scb1α21p˜N21˜P2(b1˜N1+c˜P)4.

    According to the existence of E13 in Proposition 3.1, we know

    r1(1α11˜N1)=α1p˜Pb1˜N1+c˜P,μ=sα1p˜N1b1˜N1+c˜P,1α11˜N1=sα1pb1μr1sc,˜P˜N1=sα1pb1μcμ. (4.6)

    As a result,

    A1=μc˜Pb1˜N1+c˜P+r1α11˜N1r1(1α11˜N1)b1˜N1b1˜N1+c˜P=r1+(sα1pb1μ)(scμb1μsα1p)s2cα1p,A2=r1(1α11˜N1)b1μ2sα1pr1μ(1α11˜N1)r1α11˜N1b1μ2sα1p+r1μα11˜N1+r21(1α11˜N1)2cμα1p=μ(sα1pb1μ)(r1sc(sα1pb1μ))s2cα1p>0.

    Thus, when r1=(sα1pb1μ)(b1μ+sα1pscμ)s2cα1pr1, we have A1=0. In addition, A1r1=10 is always valid.

    Note that λ1<0 is equivalent to ˜N1>1α21. Using the expression of ˜N1 in Proposition 3.1, we have

    r1sc(sα1pb1μ)r1scα11>1α21,

    i.e., r1sc(sα1pb1μ)>r1scα11α21, which is equivalent to r1>α21(sα1pb1μ)scD2r1 if D2<0. As a result, according to the Routh-Hurwitz criteria (rephrased to be user-friendly for the cases of polynomial of lower orders in Kuang et al. [15]) we have

    Proposition 4.3. Assume that r1>sα1pb1μsc and α11<α21 in (1.2). E13 is locally asymptotically stable if r1>max{r1,r1}.

    According to the the criterion in [16] (which was rephrased to be user-friendly and employed in the work of Beretta and Kuang [13]), without using eigenvalues, we obtain

    Proposition 4.4. Assume that r1>sα1pb1μsc and α11<α21 in (1.2). If r1>r1, there is a simple Hopf bifurcation at r1=r1.

    Remark 4.1. For three-species interaction model (1.2) in full-community experiment, Proposition 12 indicates that fluctuation scenarios can be found in host-parasitoid plane (N1-P plane) under some suitable parameters. In fact, taking r2=2.82,α11=1.82×105,α12=3.70×102,α21=3.97×102,α22=3.82×104,α1p=5.81,μ=0.5,b1=5.2925,b2=6.4178×105,c=3.7527 and s=7.9556×101 artificially, we have r1=0.8270,r1=0.6619. Let r1=0.8270. Figure 2 shows that the patterns of fluctuation is indeed on the N1-P plane.

    Figure 2.  Simulations of Hopf bifurcation at equilibrium E13 when r1=0.8270. Here r2=2.82,α11=1.82×105,α12=3.70×102,α21=3.97×102,α22=3.82×104,α1p=5.81,μ=0.5,b1=5.2925,b2=6.4178×105,c=3.7527 and s=7.9556×101. Black pentagram denotes the equilibrium E13=(1.0999×104,0,1.1583×104), and red circle is the attractor.

    Now, we study the stability of positive equilibrium of (1.2). For the convenience of notation, we record the positive equilibrium as E+=(N1,N2,P). Based on (3.2), after a series of calculations, the Jacobian matrix at E+ is given by

    J(E+)=[B1B2B3r2α21N2r2α22N20B4B5B6]. (4.7)

    Here

    B1=r1α11N1+α1pb1N1P(b1N1+b2N2+cP)2,B2=r1α12N1+b2α1pN1P(b1N1+b2N2+cP)2,B3=α1pN1(b1N1+b2N2)(b1N1+b2N2+cP)2,B4=sα1p(b2N2+cP)P(b1N1+b2N2+cP)2,B5=sα1pb2N1P(b1N1+b2N2+cP)2,B6=cμPb1N1+b2N2+cP.

    Furthermore, the characteristic equation of (4.7) is given by

    H(E+)(λ)λ3+C1λ2+C2λ+C3=0, (4.8)

    in which

    C1=r2α22N2B1B6,C2=B2r2α21N1(B1+B6)r2α22N2B3B4+B1B6,C3=(B1B6B3B4)r2α22N2+(B3B5B2B6)r2α21N1. (4.9)

    According to the Routh-Hurwitz criteria [15], we have

    Proposition 4.5. Suppose that sα1pb1μ>0 and the positive equilibrium E+ exists in (1.2). E+ is asymptotically stable if and if C2>0,C3>0 and C1C2C3>0.

    Remark 4.2. For three-species interaction model (1.2) in full-community experiment, Proposition 4.5 indicates that coexistence scenarios can be found under some suitable parameters, which is consistent with the experimental results in [5]. In fact, taking r1=2.90,r2=2.82×101,α11=4.92×103,α12=3.70×104,α21=3.97×104,α22=3.92×104,α1p=2.81,μ=0.50×103,b1=2.2925×105,b2=6.4178×105,c=1.7527 and s=7.9556×103 artificially, we have positive equilibrium E+=(7.6586×101,2.5502×103,1.8603), and Figure 3 shows that E+ is stable.

    Figure 3.  Simulations of the stability of positive equilibrium E+. Here r1=2.90,r2=2.82×101,α11=4.92×103,α12=3.70×104,α21=3.97×104,α22=3.92×104,α1p=2.81,μ=0.50×103,b1=2.2925×105,b2=6.4178×105,c=1.7527 and s=7.9556×103. In these cases, the positive equilibrium E+=(7.6586×101,2.5502×103,1.8603), which is the attractor and is expressed as a red pentagram.

    According to the criterion in [16], we have

    Proposition 4.6. Suppose that sα1pb1μ>0 and the positive equilibrium E+ exists. If existing a parameter value, note as θ0, such that C2(θ0)>0,C3(θ0)>0,(θ0)=0 and d(theta)dθ|θ=θ00, then a simple Hopf bifurcation occurs at θ=θ0.

    As mentioned previously, the origin equilibrium E0 in (1.2) can not be asymptotically stable since no positive solution can tend to it due to the N2 equation. Since the system is not differentiable at E0, the standard linearized techniques cannot be used to study the solution properties near it. To overcome this difficulty, we apply the technique of ratio-dependent transformations (also referred as blow-up transformations) which were effectively used in [9,10,17]. Even in the case of much simpler ratio-dependent predator-prey models, the origin is known as the source and difficulty in revealing and understanding its rich and complex dynamics [9,10,13].

    We first perform the transformation (N1,N2,P)(y1,y2,y3) where y1=N1N2,y2=N2 and y3=PN2. This transforms (1.2) to the following system:

    {dy1dt=y1(r1(1α12y2α11y1y2)α1py3b2+b1y1+cy3),dy2dt=r2y2(1α22y2α21y1y2),dy3dt=y3(sα1py1b2+b1y1+cy3μr2(1α22y2α21y1y2)). (5.1)

    The equilibria of the transformed system (5.1) include the following boundary equilibria

    V0=(0,0,0),V1=(ˆy1,0,ˆy3),V2=(0,1α22,0),V3=(D1D2,D2D,0),

    and possible positive equilibrium V4=(N1N2,N2,PN2) under corresponding conditions. Here N1,N2,P are given in Proposition 3.2-3.5. Interestingly, there is a singular line r1sy1(r2+μ)y3=0 in system (5.1), and all the points on it are the steady states (V1) of the system. Clearly, V2,V3 and V4 correspond to the E2,E12 and E+ of system (1.2) respectively, while E0 has been blown up into infinite equilibria: V0 and points (V1) on the singular line.

    The variational matrix of (5.1) evaluated at V0 is

    J(V0)=[r1000r2000r2μ].

    Note that there are two eigenvalues λ1=r1>0,λ2=r2>0. V0 is always unstable. Clearly, the y3-axis is the stable manifold of V0 while the y1y2 plane is the unstable manifold of V0.

    We can see that solution starting near the y3-axis and close to the origin will decline along y3-axis and leave the origin in a fashion tangent to the y1y2-plane.

    The variational matrix of (5.1) evaluated at V1 is

    J(V1)=[B1B2B30r20B4B5B6]. (5.2)

    Here

    B1=r1α1pˆy3(b2+cˆy3)(b2+b1ˆy1+cˆy3)2,B2=r2α12ˆy1r1α11ˆy21,B3=α1pˆy1(b2+b1ˆy1)(b2+b1ˆy1+cˆy3)2,B4=sα1pˆy3(b2+cˆy3)(b2+b1ˆy1+cˆy3)2,B5=r2α22ˆy3+r2α21ˆy1ˆy3,B6=r2μ+sα1pˆy1(b2+b1ˆy1)(b2+b1ˆy1+cˆy3)2.

    Clearly, λ1=r2>0 is one of eigenvalue of (5.2). Hence V1 is also always unstable.

    The second transformation is (N1,N2,P)(z1,z2,z3) where z1=N1P,y2=N2P and y3=P. This transforms (1.2) to the following system:

    {dz1dt=z1(r1+μr1z3(α11z1+α12z2)α1p(1+sz1)c+b1z1+b2z2),dz2dt=z2(r2+μr2z3(α21z1+α22z2)sα1pz1c+b1z1+b2z2),dz3dt=z3(sα1pz1c+b1z1+b2z2μ). (5.3)

    The equilibria of the transformed system (5.3) include the following boundary equilibria

    W0=(0,0,0),W1=(ˆz1,0,0),W2=(˜z1,˜z2,0),W3=(˜N1˜P,0,˜P),

    and possible positive equilibrium W4=(N1P,N2P,P) under corresponding conditions. Here ˆz1=c(r1+μ)α1psα1pb1μr1b1, ˜N1,˜P are given in Proposition 3.1, N1,N2,P are given in Proposition 3.2-3.5. Interestingly, there is also a singular line c+b1z1+b2z2=α1pr1r2 in system (5.3), and all the points on it are the steady states (W2) of the system. Clearly, W3 corresponds to the E13 of system (1.2), while E0 has been blown up into infinite equilibria: W0,W1 and points (W2) on the singular line.

    The variational matrix of (5.3) evaluated at W0 is

    J(W0)=[r1+μα1pc000r2+μ000μ].

    Note that one of eigenvalues λ1=r2+μ>0. W0 is always unstable.

    The variational matrix of (5.3) evaluated at W1 is

    J(W1)=[A1α1pb2ˆz1(1+sˆz1)(c+b1ˆz1)2r1α11ˆz210A2000A3], (5.4)

    and three corresponding eigenvalues of (5.4) are

    A1=r1+μα1p(c+sˆz1(2c+b1ˆz1))(c+b1ˆz1)2=(c(r1+μ)α1p)(r1b1(sα1pb1μ))α1p(scb1),A2=r2+μsα1pˆz1c+b1ˆz1=sα1psc(r1r2)b1(r2+μ)scb1,A3=μ+sα1pˆz1c+b1ˆz1=r1sc(sα1pb1μ)scb1.

    In order to make all eigenvalues of (5.4) are negative, combining with the nonnegativity of ˆz1, we split the analysis into two cases based on the precondition of sα1pb1μ>0.

    (i) Suppose that scb1>0. According to the expressions, we know ˆz1>0,A1<0,A2<0 and A3<0 means

    {c(r1+μ)α1p<0,sα1pb1μr1b1<0,sα1psc(r1r2)b1(r2+μ)<0,r1sc(sα1pb1μ)<0, (5.5)

    or

    {c(r1+μ)α1p>0,sα1pb1μr1b1>0,sα1psc(r1r2)b1(r2+μ)<0,r1sc(sα1pb1μ)<0. (5.6)

    (ii) Suppose that scb1<0. According to the expressions, we know ˆz1>0,A1<0,A2<0 and A3<0 means

    {c(r1+μ)α1p<0,sα1pb1μr1b1>0,sα1psc(r1r2)b1(r2+μ)>0,r1sc(sα1pb1μ)>0, (5.7)

    or

    {c(r1+μ)α1p>0,sα1pb1μr1b1<0,sα1psc(r1r2)b1(r2+μ)>0,r1sc(sα1pb1μ)>0. (5.8)

    However, after some analyses, we find that conditions (5.5)-(5.6) and (5.7)-(5.8) are contradict to the assumptions scb1>0 and scb1<0 respectively. Thus, W1 is always unstable.

    The variational matrix of (5.3) evaluated at W2 is

    J(W2)=[A1A2A3A4A5A600A7]. (5.9)

    Here

    A1=r1+μα1p(c+b2˜z2)(1+2s˜z1)+sα1pb1˜z21(c+b1˜z1+b2˜z2)2,A2=b2α1p˜z1(1+s˜z1)(c+b1˜z1+b2˜z2)2,A3=r1˜z1(α11˜z1+α12˜z2),A4=sα1p˜z2(c+b2˜z2)(c+b1˜z1+b2˜z2)2,A5=r2+μsα1p˜z1(c+b1˜z1)(c+b1˜z1+b2˜z2)2,A6=r2˜z2(α21˜z1+α22˜z2),

    and A7=μ+sα1p˜z1c+b1˜z1+b2˜z2=r2 based on the expression of W2. Clearly, λ1=r2>0 is one of eigenvalue of (5.9). W2 is also always unstable.

    The main purposes of this paper includes 1) the formulation and validation of a biologically more realistic and mathematically more coherent model (1.2) and 2) an initial attempt in studying the rich dynamics of model (1.2).

    As mentioned in the introduction, a main focus of population ecology studies is to identify some plausible mechanisms that leading to the extinction and coexistence of community populations. This includes total extinction or the collapse of all populations in a community. Unfortunately, existing food web models that based on prey-dependent functional response functions totally excludes the possibility of the extinction of all community populations as the bottom prey or produce population will never go extinct in those models. The key assumption of ratio-dependent models is that the growth of consumer species is a function of per-capita resource level instead of the total resource level assumed in most prey-dependent models [1,13,18]. Admittedly, ratio-dependent models are often less tractable than prey-dependent models due to the fact the origin is a non-smooth equilibrium. For this reason, there are many remaining mathematical open questions related to the global and nonlinear dynamics of model (1.2).

    An even more realistic but mathematically less tractable model of a given community population growth and interactions may involve time delays in population growth [19,20,21,22]. A first attempt on systematic studying of a special delayed ratio-dependent population model based on classic delay-independent parameter method was presented in Beretta and Kuang [23]. Realistic population models with time delays often include survival rate parameters which are usually functions of time delays [21,22]. For these delay models, most popular traditional methods of studying characteristic equations for delay models no longer effective and their study require the applications of the geometric stability switch method or its extensions presented in [20]. Not surprisingly, for the experimental results of van Veen et al. [5], we were able to show numerically that a delayed version of model (1.2) can improve the model fitting. We hope to present a mathematical study of such delayed model in the future.

    Per capita growth rate often correlates negatively with population density. The well known logistic equation for the growth of a single species incorporates this intraspecific competition. Multi-trophic models often ignore self-limitation of the consumers leading to the survival of only the fittest species. Kuang et al. found that intraspecific competition can account for the stable coexistence of many consumer species on a single resource in a homogeneous environment and resource growth rates may also play an important role in promoting coexistence of consumer species [24].

    In nature, a single resource species may contain many limiting nutrients for consumer species. In addition, producer and consumer may competing for the same essential chemical elements such as carbon, nitrogen and phosphorous [25]. These limiting nutrient contents can vary over time which affect the resource quality and hence impact the consumer growth in a complicated manner [26,27,28]. Population growth models that incorporating both resource quantity and quality are likely better match experimental observations, generating richer dynamics and promote population coexistence [29,30,31]. Such models are often termed stoichiometric population models, especially when limiting nutrients are chemical elements [32,33]. One promising and timely direction of extending the current work is to incorporating the quality dynamics of the aphid species as resource to the parasitoid.

    Yang Kuang's research is partially supported by NSF grants DMS-1615879, DEB-1930728 and an NIH grant 5R01GM131405-02. Kaifa Wang's research is partially supported by an NSFC grant 11771448.

    The authors declare that they have no competing interests.



    [1] H. I. Freedman, Deterministic Mathematical Models in Population Ecology, Marcel Dekker, Inc., New York, (1980).
    [2] G. Hardin, The competitive exclusion principle, Science, 131 (1960), 1292-1297. doi: 10.1126/science.131.3409.1292
    [3] G. F. Gause, The Struggle For Existence, Baltimore, Williams & Wilkins, (1934).
    [4] Y. Kuang, H. I. Freedman, Uniqueness of limit cycles in Gause-type models of predator-prey systems, Math. Biosci., 88 (1988), 67-84. doi: 10.1016/0025-5564(88)90049-1
    [5] F. J. F. van Veen, P. D. van Holland, H. C. J. Godfray, Stable coexistence in insect communities due to density- and trait-mediated indirect effects, Ecology, 86 (2005), 3182-3189. doi: 10.1890/04-1590
    [6] S. B. Hsu, T. W. Hwang, Y. Kuang, Rich dynamics of a ratio-dependent one prey two predator model, J. Math. Biol., 43 (2001), 377-396. doi: 10.1007/s002850100100
    [7] S. B. Hsu, T. W. Hwang, Y. Kuang, A ratio-dependent food chain model and its applications to biological control, Math. Biosci., 181 (2003), 55-83. doi: 10.1016/S0025-5564(02)00127-X
    [8] Y. Kuang, E. Beretta, Global qualitative analysis of a ratio-dependent predator-prey system, J. Math. Biol., 36 (1998), 389-406. doi: 10.1007/s002850050105
    [9] S. B. Hsu, T. W. Hwang, Y. Kuang, Global analysis of the Michaelis-Menten-type ratio-dependent predator-prey system, J. Math. Biol., 42 (2001), 489-506. doi: 10.1007/s002850100079
    [10] T. W. Hwang, Y. Kuang, Deterministic extinction effect of parasites on host populations, J. Math. Biol., 46 (2003), 17-30. doi: 10.1007/s00285-002-0165-7
    [11] T. W. Hwang, Y. Kuang, Host extinction dynamics in a simple parasite-host interaction model, Math. Biosci. Eng., 2 (2005), 743-751. doi: 10.3934/mbe.2005.2.743
    [12] P. Waltman, A Second Course in Elementary Differential Equations, Reprinted, Dover Publications, (2004).
    [13] E. Beretta, Y. Kuang, Modeling and analysis of a marine bacteriophage infection, Math. Biosci., 149 (1998), 57-76. doi: 10.1016/S0025-5564(97)10015-3
    [14] Z. E. Ma, Mathematical modeling and research of population ecology, Hefei: Anhui Education Press, (1996).
    [15] Y. Kuang, J. D. Nagy, S. E. Eikenberry, Introduction to mathematical oncology, Chapman and HallCRC, (2016).
    [16] W. M. Liu, Criterion of Hopf bifurcations without using eigenvalues, J. Math. Anal. Appl., 182 (1994), 250-256. doi: 10.1006/jmaa.1994.1079
    [17] S. Hews, S. Eikenberry, J. D. Nagy, Y. Kuang, Rich dynamics of a hepatitis B viral infection model with logistic hepatocyte growth, J. Math. Biol., 60 (2010), 573-590. doi: 10.1007/s00285-009-0278-3
    [18] H. R. Thieme, Mathematics in population biology, Princeton, (2003).
    [19] E. Beretta, Y. Kuang, Modeling and analysis of a marine bacteriophage infection with latency period, Nonlinear Anal. RWA, 2 (2001), 35-74. doi: 10.1016/S0362-546X(99)00285-0
    [20] E. Beretta, Y. Kuang, Geometric stability switch criteria in delay differential systems with delay dependent parameters, SIAM J. Math. Anal., 33 (2002), 1144-1165. doi: 10.1137/S0036141000376086
    [21] S. A. Gourley, Y. Kuang, A stage structured predator-prey model and its dependence on maturation delay and death rate, J. Math. Biol., 49 (2004), 188-200.
    [22] Y. Kuang, Delay differential equations with applications in population dynamics, Academic Press, (1993).
    [23] E. Beretta, Y. Kuang, Global analyses in some delayed ratio-dependent predator-prey systems, Nonlinear Anal. TMA, 32 (1998), 381-408. doi: 10.1016/S0362-546X(97)00491-4
    [24] Y. Kuang, W. Fagan, I. Loladze, Biodiversity, habitat area, resource growth rate and interference competition, Bull. Math. Biol., 65 (2003), 497-518. doi: 10.1016/S0092-8240(03)00008-9
    [25] R. Sterner, J. J. Elser, Ecological Stoichiometry, Princeton University Press, Princeton, NJ, (2002).
    [26] I. Loladze, Y. Kuang, J. J. Elser, Stoichiometry in producer-grazer systems: linking energy flow and element cycling, Bull. Math. Biol., 62 (2000), 1137-1162. doi: 10.1006/bulm.2000.0201
    [27] I. Loladze, Y. Kuang, J. J. Elser, W. F. Fagan, Coexistence of two predators on one prey mediated by stoichiometry, Theor. Popul. Biol., 65 (2004), 1-15. doi: 10.1016/S0040-5809(03)00105-9
    [28] X. Yang, X. Li, H. Wang, Y. Kuang, Stability and bifurcation in a stoichiometric producer-grazer model with knife edge, SIAM J. Appl. Dyn. Syst., 15 (2016), 2051-2077. doi: 10.1137/15M1023610
    [29] M. Chen, M. Fan, Y. Kuang, Global dynamics in a stoichiometric food chain model with two limiting nutrients, Math. Biosci., 289 (2017), 9-19. doi: 10.1016/j.mbs.2017.04.004
    [30] J. J. Elser, I. Loladze, A. L. Peace, Y. Kuang, Lotka re-loaded: Modeling trophic interactions under stoichiometric constraints, Ecol. Model., 245 (2012), 3-11. doi: 10.1016/j.ecolmodel.2012.02.006
    [31] A. Peace, H. Wang, Y. Kuang, Dynamics of a producer-grazer model incorporating the effects of excess food-nutrient content on grazer's growth, Bull. Math. Biol., 76 (2014), 2175-2197. doi: 10.1007/s11538-014-0006-z
    [32] J. J. Elser, Y. Kuang, Ecological stoichiometry. In: Hastings, A., Gross, L. (eds.), Encyclopedia of Theoretical Ecology, University of California Press, (2012), 718-722.
    [33] D. O. Hessen, J. J. Elser, R. W. Sterner, J. Urabe, Ecological stoichiometry: An elementary approach using basic principles, Limnol. Oceanogr., 58 (2013), 2219-2236. doi: 10.4319/lo.2013.58.6.2219
  • This article has been cited by:

    1. Zhongcai Zhu, Yuanxian Hui, Linchao Hu, The impact of predators of mosquito larvae on Wolbachia spreading dynamics , 2023, 17, 1751-3758, 10.1080/17513758.2023.2249024
    2. Yingwei Song, Tie Zhang, Jinpeng Li, Dynamical behavior in a reaction-diffusion system with prey-taxis, 2022, 2022, 1072-6691, 37, 10.58997/ejde.2022.37
  • Reader Comments
  • © 2020 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(4058) PDF downloads(306) Cited by(2)

Figures and Tables

Figures(3)  /  Tables(1)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog