Dynamics of a plant-herbivore model with a chemically-mediated numerical response

A system of two ordinary differential equations is proposed to model chemically-mediated interactions between plants and herbivores by incorporating a toxin-modifiednumerical response. This numerical response accounts for the reduction in the her-bivore's growth and reproduction due to chemical defenses from plants. It is shownthat the system exhibits very rich dynamics including saddle-node bifurcations, Hopfbifurcations, homoclinic bifurcations and co-dimension 2 bifurcations. Numerical sim-ulations are presented to illustrate the occurrence of multitype bistability, limit cycles,homoclinic orbits and heteroclinic orbits. We also discuss the ecological implicationsof the resulting dynamics.


Introduction
Ecological systems are fascinating and complex. So are interactions between plants and herbivores. On the one hand, herbivores feed on plants and on the other hand, it is frequently observed that plants can produce chemicals which either are toxic to the herbivores or can inhibit the herbivores' growth and reproduction, preventing the plants from becoming extinct [6]. It has been shown in many works (see, for example, [1,2,3,4,12,14,17]) that plant chemical defenses can have a great influence on interactions of plants and herbivores. Hence it is important to take plant chemical defenses into consideration in models of the dynamics of plant-herbivore interactions.
In general, plant chemical defense affects interactions between plants and herbivores from two aspects: (i) directly reducing the amount of plant biomass consumed by herbivores and (ii) inhibiting herbivores' growth and reproduction and increasing their mortality rate [11]. For example, related to the first aspect, tomatoes release jasmonic acid to inhibit the tobacco hornworm's digestion and stop them from getting the needed nitrogen for their own growth [7]. As another example, related to the second aspect, soybeans are resistant to attack by bruchid beetles [13]: bruchid beetle larvae are quickly killed by chemicals released by soybeans shortly after the larvae burrow beneath the seed coat.
When it comes to mathematical modeling, on the first aspect, some recent work has proposed several models with complex toxin-determined functional responses (see [5,9,10]). Rich dynamics such as oscillations from limit cycles have been observed.
We show that plant chemical defense enriches the dynamics of plant-herbivore interaction via promoting multiple types of saddle-node bifurcations, multiple types of bistability, and homoclinic bifurcations. We first formulate the model and non-dimensionalize it in Section 2. We then show the model is well posed in Section 3. We carry out the detailed mathematical analysis including existence and stability of equilibria in Section 4. Section 5 is devoted to numerically exploring various bifurcations and bistability. A summary and discussion is given in the last section.

Model formulation
Let P (t) and H(t) denote the density of plant biomass and the density of herbivore biomass at time t, respectively. Assume that the plant grows logistically when its consumer, the herbivore is absent. Plant consumption by the herbivore is modelled by a Holling type II functional response f , which takes the form f (x) = ex/(1 + hex), where e represents the rate herbivores encounter plants, and h represents the herbivore's handling time on its food. To account for the plant's chemical defenses, different from [5], we assume that a plant's chemical defense does not directly change the herbivore's consumption behavior, but rather affects the herbivore's growth and reproduction rates. More precisely, we assume that herbivore growth and reproduction rates are functions of the total grazing pressure on the plant, x = Hf (P ), and that the numerical response, denoted by g, has the specific functional form g(x) = x/(1 + (x/a) 2 ). Thus, the numerical response increases with total plant consumption, x, when x is low, but decreases if the plant consumption is high. The positive constant a is the scale parameter of the response and is the level of plant consumption for which the numerical response peaks.
Our model is then described by the following coupled differential equations Here R is the intrinsic growth rate and K is the carrying capacity of the plants. The constant B represents the conversion efficiency of consumed plants into new herbivore biomass. The constant d is the death rate for herbivore. For convenience, we first non-dimensionalize the system by setting and dropping the hats to obtain In System (2.2), the key parameter µ can be interpreted as a measure of the effectiveness of the plant at inhibiting the herbivore's digestion and growth. When µ = 0, the model reduces to the Rosenzweig-MacArthur model [15]. The Rosenzweig-MacArthur model has two boundary equilibria (0, 0) and (1, 0) and a unique positive equilibrium The equilibrium (0, 0) is always unstable, the equilibrium (1, 0) is stable if m < 1 + b and is unstable if m > 1 + b. The positive equilibrium can change from being stable to being unstable via a Hopf bifurcation resulting in the appearance of a limit cycle.
Proof. The existence and uniqueness of the solution to System (2.2) follows from standard results on ordinary differential equations [18,Theorem 1.1.8]. Note that if x(0) = 0, then x(t) = 0 for all t > 0 and if y(0) = 0, then y(t) = 0 for all t > 0. For any given initial condition (x(0), y(0)) ∈ R 2 + , we claim that the unique solution (x(t), y(t)) of (2.2) remains non-negative for t > 0. Suppose this is not true, then there must exist t 1 > 0 such that x(t) ≥ 0, y(t) ≥ 0 for t ∈ [0, t 1 ] and either x(t 1 ) = 0 with x (t 1 ) < 0 or y(t 1 ) = 0 with y (t 1 ) < 0. This is impossible since dx dt x=0 = dy dt y=0 = 0. This proves that the unique solution (x(t), y(t)) ∈ R 2 + for t > 0. Next we show the unique solution is bounded. Consider the equation with X(0) = x(0). By a comparison theorem [8, Theorem 1.4.1], we obtain Note that x(t) and y(t) are nonnegative. We then have Thus This implies that Thus, y(t) is bounded.

Stability and bifurcations
In this section, we study the dynamics of System (2.2) by considering the stability of the boundary equilibria, E 0 = (0, 0) and E 1 = (1, 0) and the existence, stability and bifurcations of the positive equilibria.  2), the equilibrium E 0 is a saddle point with a one dimensional stable manifold and a one dimensional unstable manifold; The herbivore-free equilibrium E 1 is locally asymptotically stable if and is a saddle point with a one dimensional stable manifold and a one dimensional unstable manifold if Proof. The Jacobian matrix at E 0 is given by which has two eigenvalues: λ 1 = r > 0 and λ 2 = −1 < 0. Hence, by [18,Theorem 1.1.3], E 0 is a saddle point with a one dimensional stable manifold and a one dimensional unstable manifold. Next we consider the stability of E 1 . The Jacobian matrix at E 1 is which has two eigenvalues: λ 1 = −r < 0 and λ 2 = m/(1 + b) − 1. If m < 1 + b holds, then λ 2 < 0 and thus E 1 is locally asymptotically stable, and if m > 1 + b holds, then λ 2 > 0. Again, by [ To explore the existence as well as the stability of positive equilibria, we first establish the following result on the existence of positive equilibria. Proof. A positive equilibrium, (x * , y * ), is a positive solution to the following system: 2) has an positive equilibrium (x * , y * ) with y * = r(1 − x * )(1 + bx * ). It follows from (4.3b) that This shows that m > 1 + b is a necessary condition for System (2.2) to have a positive equilibrium. In fact, it is also a sufficient condition ensuring the existence of at least one positive equilibrium for System (2.2). This follows from the fact that Q(0) = 1/µr 2 > 0 and Q This implies that Q(x) always has one negative zero and one positive zero which is larger than 1. Note that Q(x) can have at most 5 zeros, we can exclude a negative zero and a positive zero which is larger than 1. This shows that Q(x) can have at least one and at most three zeros on the interval (0, 1). That is, System (2.2) has at most three positive equilibria.
Remark 4.4. The herbivore-free equilibrium E 1 loses the stability at m = 1 + b and positive equilibria appear when m > 1 + b. This implies that a transcritical bifurcation occurs at m = 1 + b.
Next, we analyze the existence of positive equilibria for System (2.2) graphically. From (4.3), we can derive the following relation: An intersection point of the two curves given by w = F 1 (x) and w = F 2 (x) for x ∈ (0, 1) determines a positive zero of Q(x) in the interval (0, 1). As demonstrated in Figure 1 that Q(x) = 0 can have either one, or two or three positive roots in (0, 1). Figure 1 implies that System (2.2) may undergo saddle-node bifurcations. Since µ and m are the parameters associated with herbivore's growth, which may be influenced by plant's chemical defenses, we then focus on how µ and m affect the dynamics of System (2.2). Suppose all other parameter values (r and b) are fixed, then we can determine a curve in the µ-m space such that saddle-node bifurcations occur on this curve. We call this curve as a saddle-node bifurcation curve. As it is very difficult to express the x value of the positive equilibrium in terms of the parameters µ and m, we parameterize µ and m in terms of other parameters and x. On the bifurcation curve, we must have F 1 (x) = F 2 (x) and This yields is the unique positive zero of g 1 (x). The properties of the two functions m s (x) and µ s (x) are given in the following lemma.
Proof. It is straightforward to obtain (4.10), which shows that there exist 1 > 0 and 2 > 0 such that both µ s (x) and m s (x) are decreasing for x ∈ (0, 1 ) and increasing for x ∈ ( 2 ,x). A direct calculation shows that sgn(µ s (x)) = sgn(m s (x)) = sgn(G(x)) for x ∈ (0,x), where G(x) is given by Therefore, G(x) < 0 for x ∈ (0, 1 ) and G(x) > 0 for x ∈ ( 2 ,x). Note that lim x→−∞ G(x) > 0 and G(1) = −b − 1 < 0. This implies that G(x) has one negative zero and two positive zeros located in the intervals (0,x) and (x, 1), respectively. Hence, there exists a unique x c ∈ (0,x) such that The above lemma allows us to establish the following result. defines a curve in the µ−m space. Define two branches C u and C l as and then where x c is given in Lemma 4.5. Lemma 4.5 indicates that the two branches are monotone in the sense that µ s (x) and m s (x) are monotone in x for (µ, m) on C u and C l ,respectively. If (µ, m) ∈ C, then the graphs of F 1 (x) and F 2 (x) intersect tangentially at some point, say (x * , F 1 (x * )), which determines an additional positive equilibrium (x * , y * ) for System (2.2). That is, there are exactly two positive equilibria if (µ, m) is located on C. If (µ, m) is located between C u and C l (see Region III in Figure 2), then the curve given by w = F 1 (x) (locally) lies below the above mentioned point (x * , F 1 (x * )) and hence there are two intersection points for the graphs of F 1 (x) and F 2 (x) resulting in two additional positive equilibria of System (2.2). That is, there are exactly three positive equilibria when (µ, m) is in Region III. If (µ, m) is located on the other side of the curve C, then the curve given by w = F 1 (x) (locally) lies above the point (x * , F 1 (x * )) and hence there is no intersection point for the graphs of F 1 (x) and F 2 (x). Consequently, System (2.2) does not admit any positive equilibria near (x * , y * ) and there is exactly one positive equilibrium. This shows that there exists a saddle-node bifurcation curve which consists of C u and C l . Note that the two monotone branches C u and C l meet tangentially at CP with µ s (x c ) = m s (x c ) = 0. This indicates that CP is a cusp point, at which a codimension-2 bifurcation occurs [16] The corresponding characteristic equation is then given by and where g 1 (x) is defined in (4.8).
Proof. Since (x, y) is a positive equilibrium of System (2.2), then m > 1 + b and x ∈ ((m − b) −1 , 1). Note that T (x) can be rearranged as Note that ( Proof. If x ∈ (0,x), then g 1 (x) > 0 and we can rewrite D(x) as where m s (x) is defined in (4.7b). If 1 + b < m < m c , then m s (x) > m for x ∈ (0,x). Thus D(x) > 0 for x ∈ (0,x). Using (4.3), we can also write D(x) as Note that 4bm − 2b 2 > 0, g 2 (0) = 2 > 0, then the cubic function g 2 (x) must have a negative zero and two positive zeros x * 1,m and x * 2,m . By the property of cubic functions, we know that D(x) < 0 for x ∈ (x * 2,m , x * 1,m ) and Similarly, we can prove the following lemma.
is a saddle point, the other two equilibria can be stable or unstable nodes. In both cases, if at a positive equilibrium (x H , y H ) such that T (x H ) = 0 and D(x H ) > 0, then a Hopf bifurcation occurs at this positive equilibrium.
Proof. The proof follows directly from Lemmas 4.7-4.9.

Rich dynamics: numerical exploration
For each positive equilibrium (x * , y * ), since we cannot explicitly express x * in terms of system parameters, we are not able to determine the sign of T (x * ). In this section, we present some numerical simulations to explore possible dynamics that may emerge in System (2.2).
We take µ as the bifurcation parameter. According to Theorem 4.6, for each m > m c , there are two saddle-node bifurcation values, µ 1,m < µ 2,m , which give two saddle-node bifurcation points SN 1 = (x * 1,m , y * 1,m ) and SN 2 = (x * 2,m , y * 2,m ) satisfying x * 1,m > x * 2,m (See Theorem 4.6 (ii)). These two saddlenode bifurcation points divide the set of equilibria into three subsets: Moreover, by Theorem 4.10, every positive equilibrium (x * , y * ) ∈ E M is a saddle point. Since System (2.2) admits three positive equilibria when µ ∈ (µ 1,m , µ 2,m ), we see that a pair of saddle-node bifurcations occurs at SN 1 and SN 2 . Depending on the stability of the positive equilibria, we may have the following possible scenarios.
Scenario I. This corresponds to the case where (x * , y * ) is locally asymptotically stable for (x * , y * ) ∈ E U ∪ E L . That is, a stable node and a saddle point bifurcate from both saddle-node points SN 1 and SN 2 . Such a bifurcation diagram is shown in Figure 3, with a corresponding phase portrait given in Figure 4. In this case, System (2.2) exhibits Type-I bistability: a stable positive equilibrium coexists with another stable positive equilibrium (see Figure 4).
Scenario II. In this case, there exists a Hopf bifurcation point H P = (x H , y H ) ∈ E L and (x * , y * ) is locally asymptotically stable for (x * , y * ) ∈ E U ∪ E LS and (x * , y * ) is unstable for (x * , y * ) ∈ E LU . Here Note that a stable node and a saddle point bifurcate from both the saddle-node points SN 1 and SN 2 , but as µ decreases, the positive equilibrium on the node branch of SN 2 becomes unstable due to the occurrence of Hopf bifurcation. See Figure 5 for a typical bifurcation diagram. Consequently, besides Type-I bistability, Type-II bistability emerges: a stable positive equilibrium coexists with a stable limit cycle (See Figure 6).
Scenario III. In this case, the positive equilibrium (x * , y * ) is locally asymptotically stable for (x * , y * ) ∈ E U and (x * , y * ) is unstable for (x * , y * ) ∈ E L . A stable node and a saddle-point bifurcate from the saddle-node bifurcation point SN 1 , while an unstable node and a saddle-point bifurcate from the saddle-node bifurcation point SN 2 . A bifurcation diagram is presented in Figure 7. In the case that the system admits three positive equilibria, E 1 ∈ E L , E 2 ∈ E M and E 3 ∈ E U , multiple heteroclinic orbits exist: one from (1, 0) to E 3 , one from E 1 to E 3 and two from E 2 to E 3 . A phase portrait is shown in Figure 8. Scenario IV. In this case there exists a Hopf bifurcation point (x H , y H ) ∈ E U and (x * , y * ) is locally asymptotically stable for (x * , y * ) ∈ E U S and (x * , y * ) is unstable for ( As shown in Figure 9, locally, an unstable node and a saddle-point bifurcate from both saddle-node points SN 1 and SN 2 . Homoclinic bifurcation and type-III bistability. In Scenario II, there exists a Hopf bifurcation point H P on E L corresponding a critical value µ H ∈ (µ 1,m , µ 2,m ). This results in a stable limit cycle when µ is smaller than but nearby µ H . Depending on the values of the system parameters, the continuation of the limit cycle induced from Hopf bifurcation will be terminated via a homoclinic bifurcation, where the limit cycle branch touches a saddle point on E M . On the other hand, when µ < µ 1,m , System (2.2) has only one unstable positive equilibrium (x * , y * ) ∈ E LU and both boundary equilibria are unstable. Since it is a two dimensional system, and the solutions are nonnegative and bounded, there must exist at least one limit cycle around the positive equilibrium (x * , y * ). As µ increases, it is possible for the limit cycle branch to touch the saddle point resulting in a homoclinic bifurcation. This type of dynamics induces the type-III bistability: a stable positive equilibrium coexists with a half-stable homoclinic orbit. We take the parameter values as r = 3.5, b = 5.5 and m = 39.25. As shown in Figure 10, there exists a critical value µ H ≈ 16.52 at which a Hopf bifurcation occurs resulting a stable limit cycle. The continuation of the stable limit cycle terminates as µ is decreased to µ hom 2 ≈ 11.07 at which a homoclinic bifurcation occurs resulting a homoclinic orbit. Moreover, there exists another critical value µ hom 1 ≈ 7.02 at which a homoclinic bifurcation occurs resulting another homoclinic orbit. A representative phase portrait is depicted in Figure 11.

Codimension-2 bifurcations
Regarding both µ and m as bifurcation parameters, codimension−2 bifurcations including Bogdanov-Takens bifurcation and cusp point are also possible in our model. To verify this, we use the Matcont package to plot the bifurcation curves in the µ−m plane to get Figure 12.  As shown in Figure 12, between the two saddle-node bifurcation curves, there is a Hopf bifurcation curve. In addition, there are a cusp point (0.90, 13.11) and two Bogdanov-Takens points (1.05, 13.80) and (7.09, 26.29) on both the Hopf bifurcation curve and the saddle-node bifurcation curves.

Summary and discussion
In this paper we have proposed a new plant-herbivore model that includes a toxin-determined numerical response to account for the effect of plant's chemical defenses. We have shown that the toxindetermined numerical response can induce very rich dynamics including saddle-node bifurcations, Hopf bifurcation, homoclinic and Bogdanov-Takens bifurcation. This results in three types of bistability: Type-I bistability (one stable equilibrium coexists with another stable equilibrium); Type-II bistability (one stable equilibrium coexists with one stable limit cycle) and Type-III bistability (one stable equilibrium coexists with one half stable homoclinic orbit). The model's key parameter µ measures the effectiveness of the plant's inhibiting effect on the herbivore's growth and reproduction. In principle, µ increases the plant's ability to maintain a relatively high population density preventing it from extinction. Depending on other model parameters, various dynamic behaviours are possible. For instance, Figure 10 represents a bifurcation diagram for certain combination of parameters. When there is no chemical defense, that is, µ = 0, the system has a unique stable limit cycle induced by the Holling type II functional response. Both the plant population and the herbivore population oscillate about the unique positive equilibrium. If µ > 0, then there are 6 intervals in which dynamics are different: Case (i), µ ∈ (0, µ 1,m ). There is a unique positive equilibrium, which is unstable. The solutions tend to a stable limit cycle. As µ increases, the amplitude of the limit cycle decays. Case (ii), µ ∈ (µ 1,m , µ hom 1 ). There are three positive equilibria. Among them, one is a stable node, one is a saddle point, and the other is an unstable node. The solutions with initial conditions different from these equilibria approach either the stable node on E U or a stable limit cycle enclosing the unstable node on E LU . Case (iii), µ ∈ (µ hom 1 , µ hom 2 ). There are three positive equilibria and the solutions approach the stable node on E U . In addition, there are two heteroclinic orbits from the saddle point to the stable node and one heteroclinic orbit from the unstable node to the stable node. Case (iv), µ ∈ (µ hom 2 , µ H ). As in case (ii), the solutions approach either the stable node on E U or a stable limit cycle enclosing the unstable node on E LU . Case (v), µ ∈ (µ H , µ 2,m ). The solutions approach either the stable node on E U or the stable node on E LS . Case (vi), µ > µ 2,m . There is a unique positive equilibrium to which the solutions approach. Biologically, in case (i), the plant population undergoes periodic fluctuations with decaying amplitudes. In cases (ii) and (iv), the outcome is initial condition dependent: depending on the initial population, the plant population either maintains at a constant level or fluctuates about a lower level. In cases (iii) and (vi), the plant population stabilizes at a constant level, which is higher than the case with µ = 0. But case (iii) has a threshold (the stable manifold of the saddle point) for perturbations from equilibrium which lead to catastrophic plant-populations declines. In case (v), depending on the initial population, the plant population maintains at either a high level (corresponding to the equilibrium on the upper branch) or a low level (corresponding to the equilibrium on the lower branch).
Bistability exhibited in our model can be used to design suitable control strategies for resource management purpose. For instance, in the case of Type-II bistability, a stable positive equilibrium on the upper equilibrium branch coexists with a stable limit cycle that oscillates about an unstable equilibrium located on the lower equilibrium branch. If the plant population is in the basin of attraction of the limit cycle, and it is required to increase the plant population and avoid sustained oscillations, then a strategy that can shift the solution to the basin of the attraction of the stable equilibrium should be sufficient.