Persistence and Oscillations of Plant–Pollinator–Herbivore Systems

Abstract

This paper considers plant–pollinator–herbivore systems where the plant produces food for the pollinator, the pollinator provides pollination service for the plant in return, while the herbivore consumes both the food and the plant itself without providing pollination service. Based on these resource–consumer interactions, we form a plant–pollinator–herbivore model which includes the intermediary food. Using qualitative method and Kuznetsov theorem, we show global dynamics of the subsystems, uniform persistence of the whole system and periodic oscillation by Hopf bifurcation. Rigorous analysis on the system demonstrates mechanisms by which varying parameters could make the system transition between extinction of herbivore, coexistence of the three species at steady states, coexistence in periodic oscillations and extinction of pollinator. It is shown that (i) in plant–pollinator interactions, the plant would produce food; (ii) in plant–herbivore interactions, the plant would produce toxin; (iii) in the presence of both pollinator and herbivore, the plant would produce both food and toxin, and intermediate productions are analytically given by which the plant can reach its maximal density; and (iv) an appropriate toxin production could drive the herbivore into extinction, an unappropriate one would drive the pollinator into extinction, while too much toxin production will drive the plant itself into extinction. The analysis leads to explanations for experimental observations and provides new insights.

Appendices

Appendix A

Proof of Proposition 2.1

From the first equation of (6), we have ${\dot{u}}_1< u_1(r_1+\sigma _1-d_1u_1)$. Then the comparison theorem (Cosner 1996) implies

$$\begin{aligned} \lim \sup _{t\rightarrow \infty }u_1(t)\le \frac{r_1+\sigma _1}{d_1}. \end{aligned}$$

Thus, for fixed $\varepsilon > 0$, there exists $T>0$, such that $u_1(t)<\dfrac{r_1+\sigma _1}{d_1}+\varepsilon =M_1$ when $t>T$.

From the second equation of (6), we have

$$\begin{aligned} {\dot{u}}_2< u_2\left( -r_2+\frac{\sigma _2M_1}{1+u_2}\right) < -r_2u_2+\sigma _2M_1 \end{aligned}$$

when $t>T$. Then the comparison theorem implies $\lim \sup _{t\rightarrow \infty }u_2(t)\le \dfrac{\sigma _2M_1}{r_2}$.

From the first and third equations of (6), we have

$$\begin{aligned} \frac{\sigma _{31}}{\beta _{31}} {\dot{u}}_1+{\dot{u}}_3&=\frac{\sigma _{31}}{\beta _{31}}u_1\left( r_1-d_1u_1+\sigma _1\frac{\beta _0+u_2}{1+u_2+u_3}-\beta _{31}u_3\right) \\&\quad +u_3 \left( -r_3+\sigma _{31}u_1+\frac{\sigma _{32}u_1}{1+u_2+u_3}\right) \\&<\frac{\sigma _{31}}{\beta _{31}}u_1(r_1-d_1u_1+\sigma _1)-r_3u_3+\sigma _{32}u_1 \\&=-pu_1^2+qu_1-r_3u_3, \end{aligned}$$

where $p=\dfrac{\sigma _{31}d_1}{\beta _{31}}$, $q=\dfrac{\sigma _{31}r_1+\sigma _{31}\sigma _1+\sigma _{32}\beta _{31}}{\beta _{31}}$. It is obvious $-pu_1^2+qu_1\le \dfrac{q^2}{4p}$. Recall $u_1(t)<M_1$ when $t>T$. If

$$\begin{aligned} \dfrac{\sigma _{31}}{\beta _{31}} u_1(t)+u_3(t)\ge \dfrac{\sigma _{31}M_1}{\beta _{31}}+\dfrac{q^2}{4pr_3}+1 =M_3, \end{aligned}$$

then $u_3(t)\ge \dfrac{q^2}{4pr_3}+1$, which implies $\dfrac{\sigma _{31}}{\beta _{31}} {\dot{u}}_1+{\dot{u}}_3\le -r_3<0$. Thus, we have $\lim \sup _{t\rightarrow \infty }(\dfrac{\sigma _{31}}{\beta _{31}} u_1(t)+u_3(t))\le M_3$. Since solutions of (6) are nonnegative, we obtain $\lim \sup _{t\rightarrow \infty }u_3(t) \le M_3$, which implies that system (6) is dissipative. $\square$

Fig. 4

Phase plane diagrams of system (10). Stable and unstable equilibria are identified by solid and open circles, respectively. Vector fields are shown by gray arrows. The isoclines of $u_1$ and $u_3$ are represented by red and blue lines, respectively. In model (10), we fix $r_1=0.04, r_3=0.3, d_1= 0.01, \beta _0 =1/3$, $\sigma _1= 0.5$, $\sigma _{31}= 0.5, \sigma _{32}=3/16$. Numerical simulations display that the unique positive equilibrium $E_{13}(0.5435,0.1897)$ is globally asymptotically stable (Color figure online)

Appendix B

Proof of Theorem 2.4

Denote isoclines of (10) as follows:

$$\begin{aligned} \begin{aligned} l_1:~u_1&=\dfrac{1}{d_1}\left( r_1+\dfrac{\sigma _1\beta _0}{1+u_3}-\beta _{31}u_3\right) =f_1(u_3)\\ l_3:~u_1&=\dfrac{r_3(1 +u_3)}{\sigma _{31}(1+u_3)+\sigma _{32}}=f_3(u_3).\\ \end{aligned} \end{aligned}$$

Then we have

$$\begin{aligned} \dfrac{\mathrm{{d}}f_1}{\mathrm{{d}}u_3}= & {} -\dfrac{1}{d_1}\left( \dfrac{\sigma _1\beta _0}{(1+u_3)^2}+\beta _{31}\right)<0,~~ \dfrac{\mathrm{{d}}^2f_1}{\mathrm{{d}}u_3^2}= \dfrac{2\sigma _1\beta _0}{d_1(1+u_3)^3}>0 \\ \dfrac{\mathrm{{d}}f_3}{\mathrm{{d}}u_3}= & {} \dfrac{r_3\sigma _{32}}{(\sigma _{31}1+\sigma _{32}+\sigma _{31}u_3)^2}>0,~~ \dfrac{\mathrm{{d}}^2f_3}{\mathrm{{d}}u_3^2}= -\dfrac{2r_3\sigma _{31}\sigma _{32}}{(\sigma _{31}1+\sigma _{32}+\sigma _{31}u_3)^3}<0 . \end{aligned}$$

Thus, $f_1(u_3)$ is monotone decreasing and convex leftward, and $f_3(u_3)$ is monotone increasing and convex rightward as shown in Fig. 4. Let $u_3=0$. We have $f_1(0)={{\bar{u}}}_1$, and $f_3(0)=\dfrac{r_3}{\sigma _{31}+\sigma _{32}}$, which implies that $f_1(0)> f_3(0)$ if and only if $\lambda _1^{(3)} > 0$.

(i) Since $\lambda _1^{(3)} \le 0$, we have $f_1(0)\le f_3(0)$. The monotonicity and convexity of $l_1$ and $l_3$ mean that there is no positive equilibrium in system (10). Since equilibrium ${{\bar{O}}}$ has no stable manifold in int$R_+^2$, equilibrium ${{\bar{E}}}_1$ is globally asymptotically stable in int$R_+^2$.

(ii) Since $\lambda _1^{(3)} > 0$, equilibrium ${{\bar{E}}}_1$ is a saddle point and has no stable manifold in int$R_+^2$. The monotonicity and convexity of $l_1$ and $l_3$ mean that there is a unique positive equilibrium $E_{13}$ in system(10), which is asymptotically stable. By Proposition (2.3), $E_{13}$ is globally asymptotically stable in int$R_+^2$. $\square$

Appendix C

Proof of Proposition 4.1

From the second equation of (18), we obtain that $u_2^*,u_3^*$ satisfy

$$\begin{aligned} L_1:~ u_3=g_1(u_2)=-u_2+\dfrac{\sigma _2u_1^*}{r_2}-1. \end{aligned}$$

Moreover, from the first and third equations of (18), we obtain that $u_2^*,u_3^*$ satisfy

$$\begin{aligned} L_2:~ u_3=g_2(u_2)=\dfrac{\sigma _1r_2}{\beta _{31}\sigma _2u_1^*}u_2 +\dfrac{1}{\beta _{31}}\left( \dfrac{r_2\sigma _1\beta _0}{\sigma _2u_1^*}+r_1-d_1u_1^*\right) . \end{aligned}$$

On the $u_2u_3$ plane, $g_1(u_2)$ is monotonically decreasing, while $g_2(u_2)$ is monotonically increasing. The straight line $L_1$ passes through points $(\dfrac{\sigma _2u_1^*}{r_2}-1,0)$ and $(0,\dfrac{\sigma _2u_1^*}{r_2}-1)$, while $L_2$ passes through points $(-\dfrac{\sigma _2u_1^*}{\sigma _1r_2}(\dfrac{r_2\sigma _1\beta _0}{\sigma _2u_1^*}+r_1-d_1u_1^*),0)$ and $(0,\dfrac{1}{\beta _{31}}(\dfrac{r_2\sigma _1\beta _0}{\sigma _2u_1^*}+r_1-d_1u_1^*))$. Thus, there is a positive intersection point between $L_1$ and $L_2$ if and only if

$$\begin{aligned} \left\{ \begin{aligned}&\dfrac{\sigma _2u_1^*}{r_2}-1>0, \\&-\dfrac{\sigma _2u_1^*}{\sigma _1r_2}\left( \dfrac{r_2\sigma _1\beta _0}{\sigma _2u_1^*}+r_1-d_1u_1^*\right)<\dfrac{\sigma _2u_1^*}{r_2}-1, \\&\dfrac{1}{\beta _{31}}\left( \dfrac{r_2\sigma _1\beta _0}{\sigma _2u_1^*}+r_1-d_1u_1^*\right) <\dfrac{\sigma _2u_1^*}{r_2}-1. \end{aligned} \right. \end{aligned}$$

Let $u_1^* > r_2/\sigma _2$. Then we obtain the proof by (19). $\square$

Appendix D

Explicit calculation of the first Lyapunov coefficient.

Given that the explicit expression for computing the first Lyapunov coefficient involves linear A, quadratic B and cubic C terms, we are going to calculate them. The explicit form of the matrix A is in (26). The expression for B at the points $x(x_1,x_2,x_3)$ and $y(y_1,y_2,y_3)$ is

$$\begin{aligned} B(x,y)= (B_1(x,y),B_2(x,y),B_3(x,y)), \end{aligned}$$

where

$$\begin{aligned} B_1(x,y)= & {} -d_1 (1 +u_2+u_3) x_1 y_1 +(\sigma _1+ r_1 -2d_1 u_1 - \beta _{31} u_3 ) ( x_1 y_2 + x_2 y_1),\\ B_2(x,y)= & {} \sigma _2 ( x_1 y_2 + x_2 y_1) - r_2 ( x_2 y_3 + x_3 y_2+ x_2 y_2),\\ B_3(x,y)= & {} \sigma _{31} u_3 ( x_1 y_2 + x_2 y_1) +[\sigma _{31}(1 +u_2+2u_3) +\sigma _{32} ]( x_1 y_3 + x_3 y_1) \\&+ (-r_3 + \sigma _{31} u_1) ( x_2 y_3 + x_3 y_2+ x_3 y_3). \end{aligned}$$

Meanwhile, the corresponding expression for C at the points $x(x_1,x_2,x_3)$, $y(y_1,y_2,y_3)$ and $z(z_1,z_2,z_3)$ is

$$\begin{aligned} C(x,y,z)= (C_1(x,y,z),C_2(x,y,z),C_3(x,y,z)), \end{aligned}$$

where

$$\begin{aligned} C_1(x,y,z)= & {} -d_1 ( x_1^2 y_2 + x_2 y_1^2 + x_1^2 z_3 + x_3 z_1^2) \\&-\beta _{31} ( x_1 z_3^2 + x_3^2 z_1 + \sum _{i\ne j,j \ne k,k \ne i,i,j,k=1}^3 x_i y_j z_k), \\ C_2(x,y,z)= & {} 0,~~~~ C_3(x,y,z)= \sigma _{31} ( x_1 z_3^2 + x_3^2 z_1 + \sum _{i\ne j,j \ne k,k \ne i,i,j,k=1}^3 x_i y_j z_k). \end{aligned}$$

Now we proceed to compute the normalized vectors $q, {{\bar{q}}}, p$ and ${{\bar{p}}}$. The eigenvector q of A corresponding to the eigenvalue iw is

$$\begin{aligned} q=\dfrac{1}{s_0} \left( \begin{aligned}&-w^2 + (r_2u_2 +\sigma _{32} \dfrac{r_2}{\sigma _2} u_3)wi \\&\sigma _2u_2 (-\sigma _{32} u_1u_3 + wi)\\&u_3[\sigma _2 \sigma _{31} u_1u_2 + (\sigma _{31} \dfrac{\sigma _2}{r_2} u_1 + \sigma _{32})wi ]\\ \end{aligned} \right) , \end{aligned}$$

where

$$\begin{aligned} s_0= \sqrt{w^4 +(r_2u_2 +\sigma _{32} \dfrac{r_2}{\sigma _2} u_3)^2w^2 +\sigma _2^2u_2^2 (\sigma _{32}^2 u_1^2u_3^2 + w^2) +u_3^2 [\sigma _2^2 \sigma _{31}^2 u_1^2u_2^2 + (\sigma _{31} \dfrac{\sigma _2}{r_2} u_1 + \sigma _{32})^2w^2 ]}. \end{aligned}$$

The adjoint eigenvector p associated with $A^T$ (the transpose matrix of A) is

$$\begin{aligned} p= \left( \begin{aligned}&u_3\left[ \sigma _2 \sigma _{31} u_1u_2 - \left( \sigma _{31} \dfrac{\sigma _2}{r_2} u_1 + \sigma _{32}\right) wi \right] \\&u_3\left[ -d_1 \sigma _{32} u_1^2 +u_1 \sigma _1 \left( 1- \dfrac{r_2}{\sigma _2} \dfrac{\beta _0+u_2}{u_1} \right) \left( \sigma _{31} \dfrac{\sigma _2}{r_2} u_1 + \sigma _{32} \right) + \sigma _{32} \dfrac{r_2}{\sigma _2} wi\right] \\&-\sigma _2u_1u_2 \sigma _1 \left( 1- \dfrac{r_2}{\sigma _2} \dfrac{\beta _0+u_2}{u_1}\right) -w^2 + r_2u_2d_1 \dfrac{\sigma _2}{r_2} u_1^2 - \left( r_2u_2 + d_1 \dfrac{\sigma _2}{r_2} u_1^2 \right) wi\\ \end{aligned} \right) . \end{aligned}$$

The conjugate, ${{\bar{p}}}$, of p satisfying the normalization condition ${{\bar{p}}}\cdot q =1$ is

$$\begin{aligned} {{\bar{p}}}=\dfrac{1}{s_0 s_1} \left( \begin{aligned}&u_3\left[ \sigma _2 \sigma _{31} u_1u_2 + \left( \sigma _{31} \dfrac{\sigma _2}{r_2} u_1 + \sigma _{32}\right) wi \right] \\&u_3 \left[ -d_1 \sigma _{32} u_1^2 +u_1 \sigma _1 \left( 1- \dfrac{r_2}{\sigma _2} \dfrac{\beta _0+u_2}{u_1} \right) \left( \sigma _{31} \dfrac{\sigma _2}{r_2} u_1 + \sigma _{32} \right) - \sigma _{32} \dfrac{r_2}{\sigma _2} wi\right] \\&-\sigma _2u_1u_2 \sigma _1 \left( 1- \dfrac{r_2}{\sigma _2} \dfrac{\beta _0+u_2}{u_1}\right) -w^2 + r_2u_2d_1 \dfrac{\sigma _2}{r_2} u_1^2 + \left( r_2u_2 + d_1 \dfrac{\sigma _2}{r_2} u_1^2 \right) wi\\ \end{aligned} \right) , \end{aligned}$$

where

$$\begin{aligned}&\begin{aligned} \text {Re} s_1&= -u_3\sigma _2 \sigma _{31} u_1u_2 w^2 -\left( \sigma _{31} \dfrac{\sigma _2}{r_2} u_1 + \sigma _{32}\right) \left( r_2u_2 +\sigma _{32} \dfrac{r_2}{\sigma _2} u_3\right) w^2\\&\quad -\sigma _2u_2 \sigma _{32} u_1u_3^2 \left[ -d_1 \sigma _{32} u_1^2 +u_1 \sigma _1 \left( 1- \dfrac{r_2}{\sigma _2} \dfrac{\beta _0+u_2}{u_1} \right) \left( \sigma _{31} \dfrac{\sigma _2}{r_2} u_1 + \sigma _{32} \right) \right] \\&\quad + u_2 \sigma _{32} u_3 r_2 w^2\\&\quad +u_3 \sigma _2 \sigma _{31} u_1u_2 \left[ -\sigma _2u_1u_2 \sigma _1 \left( 1- \dfrac{r_2}{\sigma _2} \dfrac{\beta _0+u_2}{u_1}\right) -w^2 + u_2d_1 \sigma _2 u_1^2 \right] \\&\quad - \left( \sigma _{31} \dfrac{\sigma _2}{r_2} u_1 + \sigma _{32}\right) \left( r_2u_2 + d_1 \dfrac{\sigma _2}{r_2} u_1^2 \right) w^2,\\ \end{aligned} \\&\begin{aligned} \text {Im} s_1&= u_3\sigma _2 \sigma _{31} u_1u_2 \left( r_2u_2 +\sigma _{32} \dfrac{r_2}{\sigma _2} u_3\right) w -w^3\left( \sigma _{31} \dfrac{\sigma _2}{r_2} u_1 + \sigma _{32}\right) \\&\quad +u_3\left[ -d_1 \sigma _{32} u_1^2 +u_1 \sigma _1 \left( 1- \dfrac{r_2}{\sigma _2} \dfrac{\beta _0+u_2}{u_1} \right) \left( \sigma _{31} \dfrac{\sigma _2}{r_2} u_1 + \sigma _{32} \right) \right] \sigma _2u_2 w \\&\quad + u_2 \sigma _{32} u_1u_3^2 \sigma _{32} r_2 w\\&\quad +u_3 \sigma _2 \sigma _{31} u_1u_2 \left[ -\sigma _2u_1u_2 \sigma _1 \left( 1- \dfrac{r_2}{\sigma _2} \dfrac{\beta _0+u_2}{u_1}\right) -w^2 + u_2d_1 \sigma _2 u_1^2 \right] \\&\quad +u_3 \sigma _2 \sigma _{31} u_1u_2(r_2u_2 + d_1 \dfrac{\sigma _2}{r_2} u_1^2 ) w\\&\quad +u_3\left( \sigma _{31} \dfrac{\sigma _2}{r_2} u_1 + \sigma _{32}\right) w \left[ -\sigma _2u_1u_2 \sigma _1 \left( 1- \dfrac{r_2}{\sigma _2} \dfrac{\beta _0+u_2}{u_1}\right) -w^2 + u_2d_1 \sigma _2 u_1^2\right] .\\ \end{aligned} \end{aligned}$$

A long but straightforward computation shows that the first Lyapunov coefficient, $l_1(P^*)$, is

$$\begin{aligned} l_1(P^*)= -\frac{0.0008}{2w}<0, \end{aligned}$$

where $w = 102.2008$. Since $l_1(P^*)<0$, the bifurcation is supercritical and a unique stable limit cycle bifurcation from $P^*$ for $\sigma _2< \sigma _2^0$.