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

1.   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 $\alpha_{ii}$ is the intraspecific competition coefficient (also be interpreted as the reciprocal of carrying capacity) and $\alpha_{ij}, \, i \neq j$, is the interspecific competition coefficient, $i, j = 1,2$. $1/\mu$ represents the average lifespan of the parasitoid. Parameter $\alpha_{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:

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) \rightarrow (0^+, 0^+, 0^+)} \, R(N_{1}, N_{1}, P) = \alpha_{1p}/c$. Hence, in the following, we define $R(0, 0, 0) = \alpha_{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].

2.   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

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) \geq C(t) x$ for $x \geq 0$. Then the solution of

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) \geq C(t)x$ for $x\geq 0$ we see that $\frac{d x}{d t} \geq C(t) x$ which implies that for all $t>0$, we have $x(t) \geq x(0)e^{\int^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 $\frac{d P}{d t} = F_P(N_{1}, N_{2}, P) = Pf_P(N_{1}, N_{2}, P)$ where $f_P(N_{1}, N_{2}, P) \geq -\mu$ for positive values of $N_{1}(t), N_{2}(t)$ and $P(t).$ By Lemma 2.1 with $C(t) = -\mu$, we see that $P(t)$ stays positive.

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

where

Note that $\frac{\alpha_{1p}P}{b_{1}N_{1} + b_{2}N_{2} + cP}\leq \frac{\alpha_{1p}}{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 $\mathbf{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 $\limsup\limits_{t\rightarrow +\infty}N_{1}(t)\leq 1/\alpha_{11}\equiv N_{1}^{\infty}$ and $\limsup\limits_{t\rightarrow +\infty}N_{2}(t)\leq 1/\alpha_{22}\equiv N_{2}^{\infty}$, respectively.

Using the monotonicity of $\frac{s\alpha_{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 $\limsup\limits_{t\rightarrow +\infty}P(t)\leq \frac{s\alpha_{1p}-b_{1}\mu}{c\mu} \equiv P^{\infty}$ if $s\alpha_{1p}-b_{1}\mu > 0$, and $\lim\limits_{t\rightarrow +\infty}P(t) = 0$ if $s\alpha_{1p}-b_{1}\mu \leq 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 $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.

Proposition 2.2. Assume that $\alpha_{1p}< cr_1$ or $\alpha_{1p} \leq b_{1}\mu/s$. If $(N_{1}(t), N_{2}(t), P(t))$ is a solution of model (1.2) with initial conditions (2.1), then

Proof. We have two cases to consider: 1) $\alpha_{1p}< cr_1$ and 2) $\alpha_{1p} \leq b_{1}\mu/s$. In the following, we let $Z(t)\equiv N_1(t)+N_2(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 $t_i$ such that $Z(t_i)<1/i$ is a strictly decreasing sequence and $Z'(t_i) \leq 0$, $i = 1, 2, ...$. Since all components of the solution stay positive, we have

Hence

Since $\alpha_{1p} \leq cr_1$, we see for large enough values of $i$, we will have both $1-\frac{\alpha_{1p}}{cr_{1}}-\alpha_{11}N_{1}-\alpha_{12}N_{2}>0$ and $1-\alpha_{21}N_{1}-\alpha_{22}N_{2}>0.$ Hence $\frac{d Z (t_i)}{d t}>0$ which contradicts the fact that $t_i$ is selected to have $Z'(t_i) \leq 0$.

Consider now the case 2 with $\alpha_{1p} \leq b_{1}\mu/s$. From the proof of the case 1, the conclusion of the proposition is true if $\alpha_{1p} < cr_1.$ Hence we assume below that $cr_1 \leq \alpha_{1p} \leq b_{1}\mu/s.$ The proof of Proposition 2.1 shows that $\alpha_{1p} < b_{1}\mu/s$ implies $\lim\limits_{t\rightarrow +\infty}P(t) = 0.$ Hence for any $\delta >0$ there is a time $t_{\delta}>0$ such that $t> t_{\delta}$ implies that $P(t) < \delta.$ Let

and $\delta$ be so small such that

We claim that if for some $t_0>t_{\delta}$ such that $Z(t_0) \geq \delta_0$, then $Z(t) \geq \delta_0$ for all $t>t_0$. If this claim is false, then there is a $t_1 \geq t_o$ such that $Z(t_1) = \delta_0$ and $Z'(t_1)\leq 0$. Observe that

However,

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

If the statement that $Z(t_0) \geq \delta_0$ for some $t_0>t_{\delta}$ is not true, then $Z(t) < \delta_0$ for all $t \geq t_{\delta}$. In this situation, we have

This implies that for all $t\geq t_{\delta}$, we have $Z(t)>N_2(t)>N_2(t_{\delta})e^{gt}$ where $g = 1-(\alpha_{21}+\alpha_{22})\delta_0 >0$. Therefore $Z(t)$ will be unbounded which contradicts the statement that $Z(t) < \delta_0$ for all $t \geq t_{\delta}$. The proof of the proposition is now complete.

The proof of the second case for the Proposition 2.2 implies that if $\lim\limits_{t\rightarrow +\infty}P(t) = 0$, then the Proposition 2.2 hold. Hence we have the following result.

Theorem 2.1. If $(N_{1}(t), N_{2}(t), P(t))$ is a solution of model (1.2) with initial conditions (2.1), then

We conjecture that Proposition 2.2 remains true even if the condition $\alpha_{1p}<cr_1$ or $\alpha_{1p} \leq b_{1}\mu/s$ is removed.

Since $s\alpha_{1p}-b_{1}\mu \leq 0$ implies $\lim\limits_{t\rightarrow +\infty}P(t) = 0$, this together with Proposition 2.2 imply that if $\alpha_{1p} < b_{1}\mu/s$ in system (1.2), then its limit system is

System (2.7) is the well studied Lotka-Volterra competition model (see related sections in ^[1], ^[12] or ^[14]). Let $D_{1} = \alpha_{22}-\alpha_{12}, D_{2} = \alpha_{11}-\alpha_{21}$ and $D = \alpha_{11}\alpha_{22}-\alpha_{12}\alpha_{21}$. We have the following global results for system (1.2).

Proposition 2.3. Suppose that $\alpha_{1p} < b_{1}\mu/s$ in system (1.2), then $\lim\limits_{t\rightarrow +\infty}P(t) = 0$ and the following statements are true.

(i) When $\alpha_{11}<\alpha_{21}$ and $\alpha_{22}>\alpha_{12}$, equilibrium $E_{1} = (\frac{1}{\alpha_{11}},0, 0)$ is globally asymptotically stable and $E_{2} = (0, \frac{1}{\alpha_{22}}, 0)$ is unstable.

(ii) When $\alpha_{11}>\alpha_{21}$ and $\alpha_{22}<\alpha_{12}$, system (1.2) has no coexist equilibrium, equilibrium $E_{1}$ is always unstable and $E_{2}$ is globally asymptotically stable.

(iii) When $\alpha_{11}>\alpha_{21}$ and $\alpha_{22}>\alpha_{12}$, system (1.2) has a unique coexist equilibrium $E_{12} = (\frac{D_{1}}{D}, \frac{D_{2}}{D}, 0)$, which is globally asymptotically stable, both equilibria $E_{1}$ and $E_{2}$ are unstable.

(iv) When $\alpha_{11}<\alpha_{21}$ and $\alpha_{22}<\alpha_{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\alpha_{1p}-b_{1}\mu > 0$ is valid in system (1.2).

3.   Equilibria of system (1.2)

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

It is easy to see that system (1.2) has boundary equilibria $E_{0} = (0, 0, 0), E_{1} = (1/\alpha_{11}, 0, 0)$ and $E_{2} = (0, 1/\alpha_{22}, 0)$. Recall that $D_{1} = \alpha_{22}-\alpha_{12}, D_{2} = \alpha_{11}-\alpha_{21}$ and $D = \alpha_{11}\alpha_{22}-\alpha_{12}\alpha_{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 $D_{1}, D_{2}$ and $D$ have the same sign, i.e., $\alpha_{11}>\alpha_{21}$ and $\alpha_{22}>\alpha_{12}$, or $\alpha_{11}<\alpha_{21}$ and $\alpha_{22}<\alpha_{12}$, equilibrium $E_{12} = (\frac{D_{1}}{D}, \frac{D_{2}}{D},0)$ exists.

(ii) When $r_{1}s > \frac{s\alpha_{1p}-b_{1}\mu}{c}$, equilibrium $E_{13} = (\tilde{N}_{1}^{*},0,\tilde{P}^{*})$ exists, and

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

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

where

According $\varphi_{1}(N_{1}) = \varphi_{2}(N_{1})$, we have

Note that $\triangle = (r_{1}s\theta_{1}-\omega_{1})^{2}+4r_{1}s\theta\omega_{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

Furthermore, we have

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

3.1.   $\alpha_{11}<\alpha_{21}$ and $\alpha_{22}>\alpha_{12}$

In this case, we know $D_{1}>0$, i.e., $\theta_{1}>0$, but the sign of $\theta$ is undetermined.

When $\theta > 0$, we know $\triangle>|r_{1}s\theta_{1}-\omega_{1}|$ is always valid. So $N_{1}^{1*}>0$ and $N_{1}^{2*}<0$ in (3.5), i.e., $N_{1}^{1*}$ is the unique positive root of (3.4). Note that $\theta_{1}/\theta > 1/\alpha_{21}$ due to $\alpha_{11}<\alpha_{21}$. We know $\frac{\omega_{2}}{\omega_{1}}<N_{1}^{1*}<\frac{1}{\alpha_{21}}$ can ensure $\phi(N_{1}^{1*})>0, \varphi_{1}(N_{1}^{1*})>0$ and $\varphi_{2}(N_{1}^{1*})>0$.

When $\theta < 0$, we will analyze the existence of positive equilibrium of system (1.2) according to the sign of $\triangle$ and $r_{1}s\theta_{1}-\omega_{1}$. If $\triangle < 0$ or $r_{1}s\theta_{1}-\omega_{1}>0$, it is clear that there is no positive equilibrium in system (1.2). If $\triangle = 0$ and $r_{1}s\theta_{1}-\omega_{1}<0$, we have $N_{1}^{1*} = N_{1}^{2*} = \frac{r_{1}s\theta_{1}-\omega_{1}}{2r_{1}s\theta}\equiv N_{1}^{12*}>0$. To ensure $\phi(N_{1}^{12*})>0, \varphi_{1}(N_{1}^{12*})>0$ and $\varphi_{2}(N_{1}^{12*})>0$, we also need $\frac{\omega_{2}}{\omega_{1}}<N_{1}^{12*}<\frac{1}{\alpha_{21}}$. If $\triangle > 0$ and $r_{1}s\theta_{1}-\omega_{1}<0$, we know $N_{1}^{1*}>0, N_{1}^{2*}>0$ in (3.5). Similarly, we know $\phi(N_{1}^{k*})>0, \varphi_{1}(N_{1}^{k*})>0$ and $\varphi_{2}(N_{1}^{k*})>0$ if $\frac{\omega_{2}}{\omega_{1}}<N_{1}^{k*}< \frac{1}{\alpha_{21}}, k = 1,2$.

Summarize the above analyses, we have

Proposition 3.2. Suppose that $s\alpha_{1p}-b_{1}\mu > 0, \alpha_{11}<\alpha_{21}$ and $\alpha_{22}>\alpha_{12}$.

(i) If $\theta > 0$ and $\frac{\omega_{2}}{\omega_{1}}<N_{1}^{1*}<\frac{1}{\alpha_{21}}$, then system (1.2) has a unique positive equilibrium $E_{+} = (N_{1}^{1*}, N_{2}^{*}, P^{*})$, in which $N_{2}^{*} = \phi(N_{1}^{1*}), P^{*} = \frac{\varphi_{1}(N_{1}^{1*})}{\mu}$. Otherwise, there is no positive equilibrium in system (1.2).

(ii) If $\theta < 0, \triangle = 0$ and $r_{1}s\theta_{1}-\omega_{1}<0$, then system (1.2) has a unique positive equilibrium $E_{+} = (N_{1}^{12*}, N_{2}^{*}, P^{*})$ when $\frac{\omega_{2}}{\omega_{1}}<N_{1}^{12*}<\frac{1}{\alpha_{21}}$, in which $N_{2}^{*} = \phi(N_{1}^{12*}), P^{*} = \frac{\varphi_{1}(N_{1}^{12*})}{\mu}$. If $\theta < 0, \triangle > 0$ and $r_{1}s\theta_{1}-\omega_{1}<0$, system (1.2) may have two positive equilibria $E_{+}^{k} = (N_{1}^{k*}, N_{2}^{k*}, P^{k*})$, which depends on whether $\frac{\omega_{2}}{\omega_{1}}<N_{1}^{k*}<\frac{1}{\alpha_{21}}$ is satisfied, here $N_{2}^{k*} = \phi(N_{1}^{k*}), P^{k*} = \frac{\varphi_{1}(N_{1}^{k*})}{\mu}, k = 1,2$. Otherwise, there is no positive equilibrium in system (1.2).

3.2.   $\alpha_{11}>\alpha_{21}$ and $\alpha_{22}<\alpha_{12}$

In this case, we know $D_{1}<0$, i.e., $\theta_{1}<0$, but the sign of $\theta$ is also undetermined.

When $\theta > 0$, similar to the above analyses, we know $N_{1}^{1*}$ is the unique positive root of (3.4). However, because $\varphi_{1}(x)<0$ is always valid for $x>0$ if $\theta_{1}<0$ and $\theta>0$, there is no positive equilibrium in system (1.2).

When $\theta < 0$, since $r_{1}s\theta_{1}-\omega_{1}<0$ is always valid, we will analyze the existence of positive equilibrium of system (1.2) according to the sign of $\triangle$. If $\triangle < 0$, it is clear that there is no positive equilibrium in system (1.2). If $\triangle = 0$, we have $N_{1}^{1*} = N_{1}^{2*} = \frac{r_{1}s\theta_{1}-\omega_{1}}{2r_{1}s\theta}\equiv N_{1}^{12*}>0$. Note that $\theta_{1}/\theta > 1/\alpha_{21}$ due to $\alpha_{11}>\alpha_{21}$ and the sign of $\theta, \theta_{1}$. We know $\frac{\omega_{2}}{\omega_{1}}<N_{1}^{12*}<\frac{1}{\alpha_{21}}$ can ensure $\phi(N_{1}^{12*})>0, \varphi_{1}(N_{1}^{12*})>0$ and $\varphi_{2}(N_{1}^{12*})>0$. If $\triangle > 0$, we know $N_{1}^{1*}>0, N_{1}^{2*}>0$ in (3.5). Similarly, we know $\phi(N_{1}^{k*})>0, \varphi_{1}(N_{1}^{k*})>0$ and $\varphi_{2}(N_{1}^{k*})>0$ if $\frac{\omega_{2}}{\omega_{1}}<N_{1}^{k*}< \frac{1}{\alpha_{21}}, k = 1,2$.

Summarize the above analyses, we have

Proposition 3.3. Suppose that $s\alpha_{1p}-b_{1}\mu > 0, \alpha_{11}>\alpha_{21}$ and $\alpha_{22}<\alpha_{12}$.

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

(ii) If $\theta < 0$ and $\triangle = 0$, then system (1.2) has a unique positive equilibrium $E_{+} = (N_{1}^{12*}, N_{2}^{*}, P^{*})$ when $\frac{\omega_{2}}{\omega_{1}}<N_{1}^{12*}<\frac{1}{\alpha_{21}}$, in which $N_{2}^{*} = \phi(N_{1}^{12*}), P^{*} = \frac{\varphi_{1}(N_{1}^{12*})}{\mu}$. If $\theta < 0$ and $\triangle > 0$, system (1.2) may have two positive equilibria $E_{+}^{k} = (N_{1}^{k*}, N_{2}^{k*}, P^{k*})$, which depends on whether $\frac{\omega_{2}}{\omega_{1}}<N_{1}^{k*}<\frac{1}{\alpha_{21}}$ is satisfied, here $N_{2}^{k*} = \phi(N_{1}^{k*}), P^{k*} = \frac{\varphi_{1}(N_{1}^{k*})}{\mu}, k = 1,2$. Otherwise, there is no positive equilibrium in system (1.2).

3.3.   $\alpha_{11}>\alpha_{21}$ and $\alpha_{22}>\alpha_{12}$

In this case, we know $D_{1}>0, D>0$, i.e., $\theta_{1}>0, \theta>0$. Thus, $\triangle>|r_{1}s\theta_{1}-\omega_{1}|$ is always valid, and $N_{1}^{1*}$ is the unique positive root of (3.4). Note that $\theta_{1}/\theta < 1/\alpha_{21}$ in this case. We know $\frac{\omega_{2}}{\omega_{1}}<N_{1}^{1*}<\frac{\theta_{1}}{\theta}$ can ensure $\phi(N_{1}^{1*})>0, \varphi_{1}(N_{1}^{1*})>0$ and $\varphi_{2}(N_{1}^{1*})>0$. As a result, we have

Proposition 3.4. Suppose that $s\alpha_{1p}-b_{1}\mu > 0, \alpha_{11}>\alpha_{21}$ and $\alpha_{22}>\alpha_{12}$. If $\frac{\omega_{2}}{\omega_{1}}<N_{1}^{1*}<\frac{\theta_{1}}{\theta}$, then system (1.2) has a unique positive equilibrium $E_{+} = (N_{1}^{1*}, N_{2}^{*}, P^{*})$, in which $N_{2}^{*} = \phi(N_{1}^{1*}), P^{*} = \frac{\varphi_{1}(N_{1}^{1*})}{\mu}$. Otherwise, there is no positive equilibrium in system (1.2).

3.4.   $\alpha_{11}<\alpha_{21}$ and $\alpha_{22}<\alpha_{12}$

In this case, we know $D_{1}<0, D<0$, i.e., $\theta_{1}<0, \theta<0$. Since $r_{1}s\theta_{1}-\omega_{1}<0$ is always true, we will analyze the existence of positive equilibrium of system (1.2) according to the sign of $\triangle$. If $\triangle < 0$, it is clear that there is no positive equilibrium in system (1.2). If $\triangle = 0$, we have $N_{1}^{1*} = N_{1}^{2*} = \frac{r_{1}s\theta_{1}-\omega_{1}}{2r_{1}s\theta}\equiv N_{1}^{12*}>0$. Note that $\theta_{1}/\theta < 1/\alpha_{21}$ in this case. We know $\frac{\omega_{2}}{\omega_{1}}<N_{1}^{12*}<\frac{\theta_{1}}{\theta}$ can ensure $\phi(N_{1}^{12*})>0, \varphi_{1}(N_{1}^{12*})>0$ and $\varphi_{2}(N_{1}^{12*})>0$. If $\triangle > 0$, we have $N_{1}^{1*}>0, N_{1}^{2*}>0$ in (3.5). Similarly, we know $\phi(N_{1}^{k*})>0, \varphi_{1}(N_{1}^{k*})>0$ and $\varphi_{2}(N_{1}^{k*})>0$ if $\frac{\omega_{2}}{\omega_{1}}<N_{1}^{k*}<\frac{\theta_{1}}{\theta}, k = 1,2$.

Summarize the above analyses, we have

Proposition 3.5. Suppose that $s\alpha_{1p}-b_{1}\mu > 0, \alpha_{11}<\alpha_{21}$ and $\alpha_{22}<\alpha_{12}$. If $\triangle = 0$, then system (1.2) has a unique positive equilibrium $E_{+} = (N_{1}^{12*}, N_{2}^{*}, P^{*})$ when $\frac{\omega_{2}}{\omega_{1}}<N_{1}^{12*}<\frac{1}{\alpha_{21}}$, in which $N_{2}^{*} = \phi(N_{1}^{12*}), P^{*} = \frac{\varphi_{1}(N_{1}^{12*})}{\mu}$. If $\triangle > 0$, then system (1.2) may have two positive equilibria $E_{+}^{k} = (N_{1}^{k*}, N_{2}^{k*}, P^{k*})$, which depends on whether $\frac{\omega_{2}}{\omega_{1}}<N_{1}^{k*}<\frac{1}{\alpha_{21}}$ is satisfied, here $N_{2}^{k*} = \phi(N_{1}^{k*}), P^{k*} = \frac{\varphi_{1}(N_{1}^{k*})}{\mu}, k = 1,2$. Otherwise, there is no positive equilibrium in system (1.2).

4.   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.

Proposition 4.1. Suppose that $s\alpha_{1p}-b_{1}\mu > 0$ in (1.2).

(i) Equilibrium $E_{1}$ is a saddle.

(ii) Equilibrium $E_{2}$ is locally asymptotically stable if $\alpha_{22}<\alpha_{12}$.

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\alpha_{1p}-b_{1}\mu > 0$, $\frac{(s\alpha_{1p}-b_{1}\mu)D_{1}-b_{2}\mu D_{2}}{b_{1}D_{1}+b_{2}D_{2}}<0$ is equivalent to $\frac{s\alpha_{1p}-b_{1}\mu}{b_{2}\mu}<\frac{D_{2}}{D_{1}}$. Furthermore, recall $D = \alpha_{11}\alpha_{22}-\alpha_{12}\alpha_{21}$, we can conclude that the trace of the upper left-hand 2$\times$2 matrix is always negative and its determinant is positive if all $D_{1},D_{2}, D$ are positive. Thus, we have

Proposition 4.2. Suppose that $s\alpha_{1p}-b_{1}\mu > 0$ in (1.2). Equilibrium $E_{12}$ is locally asymptotically stable if $D_{1}>0, D_{2}>0, D>0$ and $\frac{s\alpha_{1p}-b_{1}\mu}{b_{2}\mu}<\frac{D_{2}}{D_{1}}$.

The Jacobian matrix of (1.2) at $E_{13}$ is

Here

Clearly, $\lambda_{1} = r_{2}(1-\alpha_{21}\tilde{N}_{1}^{*})$ is one of eigenvalue of (4.4), and the other two eigenvalues of (4.4) satisfying

Here

According to the existence of $E_{13}$ in Proposition 3.1, we know

As a result,

Thus, when $r_{1} = \frac{(s\alpha_{1p}-b_{1}\mu)(b_{1}\mu+s\alpha_{1p}-sc\mu)}{s^{2}c\alpha_{1p}}\equiv r_{1}^{*}$, we have $A_{1} = 0$. In addition, $\frac{\partial A_{1}}{\partial r_{1}} = 1\neq 0$ is always valid.

Note that $\lambda_{1}<0$ is equivalent to $\tilde{N}_{1}^{*}>\frac{1}{\alpha_{21}}$. Using the expression of $\tilde{N}_{1}^{*}$ in Proposition 3.1, we have

i.e., $r_{1}sc-(s\alpha_{1p}-b_{1}\mu)>\frac{r_{1}sc\alpha_{11}}{\alpha_{21}}$, which is equivalent to $r_{1}>- \frac{\alpha_{21}(s\alpha_{1p}-b_{1}\mu)}{scD_{2}}\equiv r_{1}^{**}$ if $D_{2}<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 $r_{1}>\frac{s\alpha_{1p}-b_{1}\mu}{sc}$ and $\alpha_{11}<\alpha_{21}$ in (1.2). $E_{13}$ is locally asymptotically stable if $r_{1}>\max\{r_{1}^{*},r_{1}^{**}\}$.

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 $r_{1}>\frac{s\alpha_{1p}-b_{1}\mu}{sc}$ and $\alpha_{11}<\alpha_{21}$ in (1.2). If $r_{1}^{*}>r_{1}^{**}$, there is a simple Hopf bifurcation at $r_{1} = r_{1}^{*}$.

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 ($N_{1}$-$P$ plane) under some suitable parameters. In fact, taking $r_{2} = 2.82, \alpha_{11} = 1.82\times 10^{-5}, \alpha_{12} = 3.70\times 10^{2}, \alpha_{21} = 3.97\times 10^{2}, \alpha_{22} = 3.82\times 10^{-4}, \alpha_{1p} = 5.81, \mu = 0.5, b_{1} = 5.2925, b_{2} = 6.4178\times 10^{-5},c = 3.7527$ and $s = 7.9556\times 10^{-1}$ artificially, we have $r_{1}^{*} = 0.8270, r_{1}^{**} = 0.6619$. Let $r_{1} = 0.8270$. Figure 2 shows that the patterns of fluctuation is indeed on the $N_{1}$-$P$ plane.

Now, we study the stability of positive equilibrium of (1.2). For the convenience of notation, we record the positive equilibrium as $E_{+} = (N_{1}^{*},N_{2}^{*},P^{*})$. Based on (3.2), after a series of calculations, the Jacobian matrix at $E_{+}$ is given by

Here

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

in which

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

Proposition 4.5. Suppose that $s\alpha_{1p}-b_{1}\mu > 0$ and the positive equilibrium $E_{+}$ exists in (1.2). $E_{+}$ is asymptotically stable if and if $C_{2}>0, C_{3}>0$ and $\triangle\equiv C_{1}C_{2}-C_{3}>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 $r_{1} = 2.90,r_{2} = 2.82\times 10^{-1}, \alpha_{11} = 4.92\times 10^{-3}, \alpha_{12} = 3.70\times 10^{-4}, \alpha_{21} = 3.97\times 10^{-4}, \alpha_{22} = 3.92\times 10^{-4}, \alpha_{1p} = 2.81, \mu = 0.50\times 10^{-3}, b_{1} = 2.2925\times 10^{-5}, b_{2} = 6.4178\times 10^{-5},c = 1.7527$ and $s = 7.9556\times 10^{-3}$ artificially, we have positive equilibrium $E_{+} = (7.6586\times 10^{-1}, 2.5502\times 10^{3}, 1.8603)$, and Figure 3 shows that $E_{+}$ is stable.

According to the criterion in ^[16], we have

Proposition 4.6. Suppose that $s\alpha_{1p}-b_{1}\mu > 0$ and the positive equilibrium $E_{+}$ exists. If existing a parameter value, note as $\theta_{0}$, such that $C_{2}(\theta_{0})>0,\, C_{3}(\theta_{0})>0,\, \triangle(\theta_{0}) = 0$ and $\frac{d \triangle(theta)}{d \theta}|_{\theta = \theta_{0}}\neq 0$, then a simple Hopf bifurcation occurs at $\theta = \theta_{0}$.

5.   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].

We first perform the transformation $(N_{1},N_{2},P)\rightarrow (y_{1},y_{2},y_{3})$ where $y_{1} = \frac{N_{1}}{N_{2}}, y_{2} = N_{2}$ and $y_{3} = \frac{P}{N_{2}}$. This transforms (1.2) to the following system:

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

and possible positive equilibrium $V_{4} = (\frac{N_{1}^{*}}{N_{2}^{*}},N_{2}^{*},\frac{P^{*}}{N_{2}^{*}})$ 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}+\mu)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 $\lambda_{1} = r_{1}>0, \lambda_{2} = r_{2}>0$. $V_{0}$ is always unstable. Clearly, the $y_3$-axis is the stable manifold of $V_{0}$ while the $y_1y_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_1y_2$-plane.

The variational matrix of (5.1) evaluated at $V_{1}$ is

Here

Clearly, $\lambda_{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)\rightarrow (z_{1},z_{2},z_{3})$ where $z_{1} = \frac{N_{1}}{P}, y_{2} = \frac{N_{2}}{P}$ and $y_{3} = P$. This transforms (1.2) to the following system:

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

and possible positive equilibrium $W_{4} = (\frac{N_{1}^{*}}{P^{*}},\frac{N_{2}^{*}}{P^{*}},P^{*})$ under corresponding conditions. Here $\hat{z}_{1} = \frac{c(r_{1}+\mu)-\alpha_{1p}}{s\alpha_{1p}-b_{1}\mu-r_{1}b_{1}}$, $\tilde{N}_{1}^{*},\tilde{P}^{*}$ are given in Proposition 3.1, $N_{1}^{*},N_{2}^{*},P^{*}$ are given in Proposition 3.2-3.5. Interestingly, there is also a singular line $c+b_{1}z_{1}+b_{2}z_{2} = \frac{\alpha_{1p}}{r_{1}-r_{2}}$ in 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 $\lambda_{1} = r_{2}+\mu>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 $\hat{z}_{1}$, we split the analysis into two cases based on the precondition of $s\alpha_{1p}-b_{1}\mu>0$.

(i) Suppose that $sc-b_{1}>0$. According to the expressions, we know $\hat{z}_{1}>0,A_{1}<0,A_{2}<0$ and $A_{3}<0$ means

or

(ii) Suppose that $sc-b_{1}<0$. According to the expressions, we know $\hat{z}_{1}>0,A_{1}<0,A_{2}<0$ and $A_{3}<0$ means

or

However, after some analyses, we find that conditions (5.5)-(5.6) and (5.7)-(5.8) are contradict to the assumptions $sc-b_{1}>0$ and $sc-b_{1}<0$ respectively. Thus, $W_{1}$ is always unstable.

The variational matrix of (5.3) evaluated at $W_{2}$ is

Here

and $A_{7} = -\mu+\frac{s\alpha_{1p}\tilde{z}_{1}}{c+b_{1}\tilde{z}_{1}+b_{2}\tilde{z}_{2}} = r_{2}$ based on the expression of $W_{2}$. Clearly, $\lambda_{1} = r_{2}>0$ is one of eigenvalue of (5.9). $W_{2}$ is also always unstable.

6.   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 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.

Acknowledgments

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.

Conflict of interest

The authors declare that they have no competing interests.