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

Abstract: 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.


Introduction
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: where N 1 and N 2 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 N 1 . The parameter r i 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 α i j , i j, 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 N 1 , while b 1 is the parameter controlling the reduction in parasitoid attack rate with increasing host number (N 1 ), which may be interpreted as correlated with the searching and handling time, and b 2 can be interpreted as the time wasted when a parasitoid encounters an unsuitable aphid (N 2 ), 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: (  viciae and A. ervi in the full-community experiment in [5]. Since r 1 , r 2 , α 11 , α 22 and µ are species-specific growth parameters, we adopt the estimated values in [5]. Specifically, r 1 = 3.22/week, r 2 = 2.82/week, α 11 = 3.82 × 10 −4 , α 22 = 3.82 × 10 −4 . For other parameters, according to [5], we take α 12 = 3.70 × 10 −4 , α 21 = 3.97 × 10 −4 , α 1p = 2.81 × 10 −1 /week, µ = 6.34 × 10 −1 /week, b 1 = 2.33 × 10 −2 , b 2 = 4.34 × 10 −2 , 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 × 10 −4 , α 21 = 5.7116 × 10 −4 , α 1p = 3.5190 × 10 −2 /week, µ = 0.50/week, b 1 = 5.2925 × 10 −5 , b 2 = 6.4178 × 10 −5 , c = 1.7527 × 10 −2 and s = 7.9556 × 10 −3 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 is not defined at (0, 0, 0). However, lim (N 1 ,N 1 ,P)→(0 + ,0 + ,0 + ) R(N 1 , N 1 , 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 ratiodependent models with only two populations [8,9,10,11]. S S E N1 = 1.6146 × 10 6 S S E N2 = 4.8173 × 10 6 S S E N2 = 3.4856 × 10 6 S S E P = 1.3738 × 10 5 S S E P = 7.7021 × 10 2 S S E T otal = 6.0014 × 10 6 S S E T otal = 5.1010 × 10 6

Dissipativity
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 N 1 (0) > 0, N 2 (0) > 0 and P(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. Before proving the above proposition, we would like to establish the following general result. 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 x ≥ 0 we see that dx dt ≥ C(t)x which implies that for all t > 0, we have x(t) ≥ x(0)e t 0 C(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 dP dt = F P (N 1 , N 2 , P) = P f P (N 1 , N 2 , P) where f P (N 1 , N 2 , P) ≥ −µ for positive values of N 1 (t), N 2 (t) and P(t). By Lemma 2.1 with C(t) = −µ, we see that P(t) stays positive.
We now consider the positivity of N 1 (t). Observe that where Note that c is always valid because all variables are nonnegative. As a result, we have By Lemma 2.1 with C(t) = C 1 (t), we see that N 1 (t) stays positive. The positivity of N 2 (t) can be established similarly. Now we present the arguments for ultimate boundedness of solutions of (1.2) in R 3 + . Since all components of a solution of (1.2) are positive, using the first and second equations of (1.2), we have As a result, we have lim sup t→+∞ N 1 (t) ≤ 1/α 11 ≡ N ∞ 1 and lim sup t→+∞ N 2 (t) ≤ 1/α 22 ≡ N ∞ 2 , respectively. Using the monotonicity of sα 1p N 1 P b 1 N 1 +cP on N 1 and above analysis, we have Thus, using the third equation of (1.2), we have As a result, we have lim sup 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 N 2 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.
Proof. We have two cases to consider: 1) α 1p < cr 1 and 2) α 1p ≤ b 1 µ/s. In the following, we let Consider first case 1. We prove it by contradiction. If the proposition is false, then there is a strictly increasing sequence of positive values t i such that Z(t i ) < 1/i is a strictly decreasing sequence and Z (t i ) ≤ 0, i = 1, 2, .... Since all components of the solution stay positive, we have Since α 1p ≤ cr 1 , we see for large enough values of i, we will have both 1 − α 1p and δ be so small such that We claim that if for some t 0 > t δ such that Z(t 0 ) ≥ δ 0 , then Z(t) ≥ δ 0 for all t > t 0 . If this claim is false, then there is a t 1 ≥ t o such that Z(t 1 ) = δ 0 and Z (t 1 ) ≤ 0. Observe that This is a contradiction which proves our claim. Clearly this claim implies the conclusion of the proposition.
If the statement that Z(t 0 ) ≥ δ 0 for some t 0 > t δ is not true, then Z(t) < δ 0 for all t ≥ t δ . In this situation, we have This implies that for all t ≥ t δ , we have Therefore Z(t) will be unbounded which contradicts the statement that Z(t) < δ 0 for all t ≥ t δ . The proof of the proposition is now complete.
The proof of the second case for the Proposition 2.2 implies that if lim t→+∞ P(t) = 0, then the Proposition 2.2 hold. Hence we have the following result.
We conjecture that Proposition 2.2 remains true even if the condition α 1p < cr 1 or α 1p ≤ b 1 µ/s is removed. Since System (2.7) is the well studied Lotka-Volterra competition model (see related sections in [1], [12] or [14]). Let We have the following global results for system (1.2).
(iv) When α 11 < α 21 and α 22 < α 12 , the unique coexist equilibrium E 12 is a saddle, both equilibria E 1 and E 2 are stable, i.e., the bistability phenomenon of initial value dependence will occur.
In the rest of this paper, we assume that sα 1p − b 1 µ > 0 is valid in system (1.2).
Proposition 3.1. The following are true for system (1.2).
(i) If θ > 0, then there is no positive equilibrium in system (1.2).

Stability of non-origin equilibria
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 E 0 , and the standard linearization approach cannot be used at E 0 . To avoid this difficulty for now, this section will focus on the non-origin equilibria, including boundary equilibria E 1 , E 2 , E 12 , E 13 in Proposition 3.1, and the possible positive equilibria E + , E k + , k = 1, 2. The Jacobian matrix of (1.2) at E 1 and E 2 are and respectively. Note that the eigenvalues of (4.1) and (4.2) lie on the diagonal. Hence, we have the following results.
The Jacobian matrix of (1.2) at E 12 is Clearly, the element in the lower right-hand corner of (4.3) is one of eigenvalue of (4.3). When sα 1p − b 1 µ > 0, 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 D 1 , D 2 , D are positive. Thus, we have

Dynamics near the origin equilibrium E 0
As mentioned previously, the origin equilibrium E 0 in (1.2) can not be asymptotically stable since no positive solution can tend to it due to the N 2 equation. Since the system is not differentiable at E 0 , 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].
This transforms (1.2) to the following system: The equilibria of the transformed system (5.1) include the following boundary equilibria ) under corresponding conditions. Here N * 1 , N * 2 , P * are given in Proposition 3.2-3.5. Interestingly, there is a singular line r 1 sy 1 − (r 2 + µ)y 3 = 0 in system (5.1), and all the points on it are the steady states (V 1 ) of the system. Clearly, V 2 , V 3 and V 4 correspond to the E 2 , E 12 and E + of system (1.2) respectively, while E 0 has been blown up into infinite equilibria: V 0 and points (V 1 ) on the singular line.
The variational matrix of (5.1) evaluated at V 0 is Note that there are two eigenvalues λ 1 = r 1 > 0, λ 2 = r 2 > 0. V 0 is always unstable. Clearly, the y 3 -axis is the stable manifold of V 0 while the y 1 y 2 plane is the unstable manifold of V 0 .
We can see that solution starting near the y 3 -axis and close to the origin will decline along y 3 -axis and leave the origin in a fashion tangent to the y 1 y 2 -plane.
The variational matrix of (5.1) evaluated at V 1 is Here Clearly, λ 1 = r 2 > 0 is one of eigenvalue of (5.2). Hence V 1 is also always unstable.
The second transformation is (N 1 , N 2 , P) → (z 1 , z 2 , z 3 ) where z 1 = N 1 P , y 2 = N 2 P and y 3 = P. This transforms (1.2) to the following system:  (5.3), and all the points on it are the steady states (W 2 ) of the system. Clearly, W 3 corresponds to the E 13 of system (1.2), while E 0 has been blown up into infinite equilibria: W 0 , W 1 and points (W 2 ) on the singular line.
The variational matrix of (5.3) evaluated at W 0 is Note that one of eigenvalues λ 1 = r 2 + µ > 0. W 0 is always unstable. The variational matrix of (5.3) evaluated at W 1 is and three corresponding eigenvalues of (5.4) are In order to make all eigenvalues of (5.4) are negative, combining with the nonnegativity ofẑ 1 , we split the analysis into two cases based on the precondition of sα 1p − b 1 µ > 0.
The variational matrix of (5.3) evaluated at W 2 is Here c+b 1z1 +b 2z2 = r 2 based on the expression of W 2 . Clearly, λ 1 = r 2 > 0 is one of eigenvalue of (5.9). W 2 is also always unstable.

Discussion
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 delayindependent 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.