PEST CONTROL BY GENERALIST PARASITOIDS: A BIFURCATION THEORY APPROACH

. Magal et al. [13] studied both spatial and non-spatial host-parasitoid models motivated by biological control of horse-chestnut leaf miners that have spread through Europe. In the non-spatial model, they considered pest control by predation of leaf miners by a generalist parasitoid with a Holling type II functional (Monod) response. They showed that there can be at most six equilibrium points and discussed local stability. We revisit their model in the non-spatial case, identify cases missed in their investigation and discuss conse-quences for possible pest control strategies. Both the local stability of equilibria and global properties are considered. We use a bifurcation theoretical approach and provide analytical expressions for fold and Hopf bifurcations and for the criticality of the Hopf bifurcations. Our numerical results show very interesting dynamics resulting from codimension one bifurcations including: Hopf, fold, transcritical, cyclic-fold, and homoclinic bifurcations as well as codimension two bifurcations including: Bautin and Bogdanov-Takens bifurcations, and a codimension three Bogdanov-Takens bifurcation.


1.
Introduction. Leaf miners [6] are the larvae of various insects such as moths, beetles, and flies. They are harmful invasive species and live in and eat the leaf tissue of plants. Once they emerge, if it is on an edible plant such as spinach or kale, leaf miner feeding tunnels give the plants an unpleasant appearance so that the plants may be unmarketable resulting in financial loss for farmers.
To control this harmful pest, at least two options are available: (i) chemical control [17] and (ii) biological control [10]. For chemical control, pesticides (e.g., insecticidal soap, silica gel, Rotenone, Pyrethrins, etc.) are used. However, the use of pesticides can be toxic to the environment and risk the health of humans and other animals. Due to such risks, instead, biological control such as the introduction of natural enemies in the form of predators, parasitoids, and/or pathogens is often used for the management of such pests.
Magal et al. [13] investigated a spatial and non-spatial host-parasitoid model for the biological control of a leaf miner population, e.g., Cameraria ohridella (commonly known as micro-moths) that rapidly spread across Europe attacking horsechestnut trees. As a control strategy, they considered the predation of leaf miners by a generalist parasitoid population with a Holling type II (Monod) response function [7,8]. In this paper, we focus on their non-spatial model of leaf miner and parasitoid interaction: where N and P denote the densities of the host (leaf miners) and the parasitoid (generalist predator) populations at time t. All of the parameters r 1 , r 2 , K 1 , K 2 , E, h, and γ are positive constants. In the absence of the parasitoid the host is assumed to grow logistically. The parasitoids are assumed to be generalists and so they are also assumed to grow logistically in the absence of the host population. The intrinsic growth rates are denoted by r 1 and r 2 and the carrying capacities by K 1 and K 2 . E and h in response function Φ(N ) represent the encounter rate and harvesting time, respectively. The parasitoids are assumed to predate on the host population. Hence, the growth of the host population is limited by the parasitoid population. The parameter γ denotes the conversion efficiency of parasitized hosts into new parasitoids. For notational simplicity, we use the following mathematical form for the response function: Magal et al. [13] considered the number of coexistence equilibria where host and parasitoid populations coexist using host and parasitoid nullclines, carried out local stability analysis, and described the host-parasitoid dynamics using phase portraits. They also discussed possible bifurcations without detailed analysis. In this paper, we revisit their non-spatial model doing a more detailed bifurcation analysis and identify a missed case. We then consider the ramifications for pest control strategies. Throughout this paper, due to the significance of the host (leaf miners) and the parasitoid carrying capacities for pest control, we primarily focus on how the host-parasitoid dynamics change with respect to changes of the carrying capacities.
This paper has six sections. In Section 2, we begin by analyzing host and parasitoid nullclines for system (1.1). The nullcline analysis shows that the system has at most six equilibria; three boundary equilibria that always exist and up to three positive (or coexistence) equilibria. We also determine the existence and the number of coexistence equilibria in terms of the relative values of the parameters. This involves finding analytical expressions for fold bifurcations and a cusp point at which two fold curves intersect. In Section 3, by means of a linearized stability analysis, we study the local stability of all of the equilibria. Furthermore, we obtain 0 N P K 2 Figure 1. Graphs of the host nullcline, F (N ), drawn in a thick solid curve and the parasitoid nullcline, G(N ), drawn in a thin solid curve for model (2.1) with the Holling type II functional response.
necessary and sufficient conditions for the global stability of the pest-free equilibrium and sufficient conditions for the global stability of the coexistence equilibrium when it is unique. Section 4 contains the Hopf bifurcation analysis. In this section, we prove the existence of Hopf bifurcations and determine analytical expressions for the two-parameter curves of Hopf bifurcations in the K 1 K 2 -plane. We also determine an analytic expression for the criticality of the Hopf bifurcations and for codimension two Bautin bifurcations. We show that the system undergoes not only a codimension two Bogdanov-Takens (BT) bifurcation, but also a codimension three BT bifurcation that occurs when the BT bifurcation occurs at the cusp point. In Section 5, we provide one and two-parameter bifurcation diagrams, with the corresponding phase portraits. Our numerical work in this section illustrates homoclinic bifurcations and the existence of multiple limit cycles via cyclic fold bifurcations, as well as the other bifurcations that we determined analytically in the previous sections. Finally, in Section 6, we summarize our results and discuss the implications for pest control focussing on the effect of variations in K 1 and K 2 . We explore the missing case in [13] with this in mind. Some of the calculations were assisted by the use of Maple [14]. All figures were generated by using Ipe [9] or MATLAB [15] and/or XPPAUT [4].
2. Nullclines and equilibria. We consider the following system with functional response Φ(N ) given in (1.3): with Φ(N ) = aN/(1 + bN ). The host nullclines are N = 0, P = F (N ) (2.2) and the parasitoid nullclines are The host nullcline P = F (N ) (thick solid curve in Figure 1) has the following properties : Thus, the graph of the host nullcline has P -intercept at the point (0, r 1 /a H ) and N -intercept at (K 1 , 0). The graph is concave down with one local maximum at in the interior of the first quadrant, when K 1 > 1/b. When K 1 ≤ 1/b, the graph of F (N ) is a monotonically decreasing function for N ≥ 0. Seo and Wolkowicz [18] demonstrate that for the prey nullcline, F (N, where r > 0 and K > 0, the value of F (N, K) at a local maximum (if exists) moves up and to the right as K increases and at a local minimum (if exists) moves up and to the left under the assumptions Φ(N ) ∈ C 2 ([0, ∞), R), Φ(0) = 0, Φ (N ) > 0, and Φ (N ) < 0, for all nonnegative N . The parasitoid nullcline P = G(N ) (the thin solid curve in Figure 1) has the following properties: The parasitoid nullcline is a positive increasing function with horizontal asymptote P = K 2 (γa + r 2 b)/(r 2 b) in the interior of the first quadrant. It is concave down and has P -intercept at the point (0, K 2 ). Equilibria can be found at the points of intersection of the host and parasitoid nullclines where dN/dt and dP/dt are both zero. System (2.1) has at most six equilibria, Figure 2. The number of coexistence equilibria, denoted by E * , in different regions of the K 1 K 2 -plane. C 1 (thick solid curve) and C 2 (thin solid curve) indicate branches of saddle-node bifurcations and the dot-dashed line indicates a branch of transcritical bifurcations.
In the white regions, no coexistence equilibria exist and the preyextinction equilibrium E K2 is globally asymptotically stable. See (2.11) for θ 1,2 and κ. where , to obtain a cubic equation (2.8) There are at most three equilibria in the interior of the first quadrant where both the host population and the parasitoid population coexist. We call them coexistence equilibria. They exist provided that N * > 0 and F (N * ) > 0. This is equivalent to assuming 0 < N * < K 1 . Figure 2 illustrates the number of coexistence equilibria in different regions of the K 1 K 2 -plane. In the figure, C 1 (thick solid curve) and C 2 (thin solid curve), where two coexistence equilibria appear and separate or collide and disappear, indicate branches of saddle-node (or fold) bifurcations: where K 1 ≥ 8γa−r2b b(γa+r2b) and In particular, a cusp point occurs at (2.10) Note that the cusp point appears in the interior of the first quadrant if 8aγ −r 2 b > 0 and it is located below the straight line, K 2 = r 1 /a (see Figure 2). In Figure 2, θ 1,2 and κ are given by The reader may refer to Appendix A to see how Figure 2 is produced.
3. Stability analysis. In order to determine the local stability of the equilibria, we use the Jacobian matrix given by, By a linearized stability analysis, E 0 = (0, 0) is an unstable node and E K1 = (K 1 , 0) is a saddle point. The Jacobian at E K2 = (0, K 2 ), where only the parasitoids survive, is Thus, E K2 is a stable node if K 2 > r 1 /a and a saddle point if K 2 < r 1 /a. The Jacobian at the coexistence equilibrium E * = (N * , P * ) can be expressed as where P * = F (N * ) = G(N * ), and the characteristic equation is and a saddle point if F (N * ) > G (N * ). Table 1 summarizes the stability conditions for equilibria in model (2.1). From Table 1, we see that the model possesses a transcritical bifurcation at Proof. The N -and P -axes of (2.1) are invariant. Since the vector field is C 1 , solutions with positive initial conditions remain positive. G(N ) is an increasing function with finite horizontal asymptote equal to m. For any solution with positive initial conditions, by (2.1)(a), for any t ≥ 0, if N (t) ≥ K 1 , N (t) < 0 and by (2.1(b)), if P (t) > m, then P (t) < 0. Therefore, solutions are bounded and Γ is positively invariant.
Equilibrium It remains to show that all solutions eventually enter Γ. By the Poincaré-Bendixson Theorem [16], the ω-limit set of any orbit of this system must contain a periodic orbit or an equilibrium point. However, all equilibrium points lie inside Γ, the N coordinates on any periodic orbit must be less than K 1 , and the P coordinates on any periodic orbit must be less than or equal to m. Hence, any periodic orbit must lie entirely inside Γ, and so any orbit with positive initial conditions must enter Γ.
The ultimate goal of a pest control program is to eradicate the pest populations. Here, we discuss the global stability of the host-free (i.e., the pest-free) equilibrium E K2 = (0, K 2 ).
, then the pest-free equilibrium E K2 is globally asymptotically with respect to all solutions initiating in the set {(N, P ) N ≥ 0, P > 0}.
Proof. No coexistence equilibria exist, when E 0 = (0, 0) is an unstable node, and E K1 = (K 1 , 0) is a saddle point with stable manifold along the N -axis. In this case, E K2 is locally asymptotically stable and is the only stable equilibrium point. Since by Theorem 3.1, all orbits are bounded, using the Poincaré-Bendixson Theorem and phase portrait analysis, illustrated in Figure 3, all orbits with N (0) ≥ 0 and P (0) > 0 converge to the unique stable boundary equilibrium E K2 . Theorem 3.3. If F (N ) < 0 for all 0 < N < K 1 or either r 1 ≤ r 2 or (r 1 > r 2 and K 1 ≤ r1+r2 b(r1−r2) ), then system (2.1) has no nontrivial periodic orbits in the interior of R 2 + .
Proof. Write system (2.1) in the form, where Define the auxiliary function, .
Hence, if F (N ) < 0 for 0 < N < K 1 , by the Dulac criterion [3,16], system (2.1) has no nontrivial periodic orbits in the interior of R 2 + . Consider the direction field. If a periodic orbit exists, it must cross the curve G(N ) horizontally and the curve F (N ) vertically and hence must be entirely con- Hence, if r 1 ≤ r 2 or if r 1 > r 2 and K 1 ≤ r1+r2 b(r1−r2) , again by the Dulac criterion, system (2.1) has no nontrivial periodic orbits in the interior of R 2 + .
a , then there is a unique coexistence equilibrium, E * , and it is globally asymptotically stable.
Proof. (i) The result follows from standard phase plane analysis (see Figure 3(a)).
(ii) The result follows from the Poincaré-Bendixson Theorem, since the boundary equilibria are all unstable in this case, none of their ω-limit sets intersect int(R 2 + ), and by Theorem 3.3, there are no periodic orbits.
Corollary 3.2. If there exists a unique interior equilibrium E * and either r 1 ≤ r 2 or r 1 > r 2 and , then E * is globally asymptotically stable with respect to int(R 2 + ). Proof. The result follows from the Poincaré-Bendixson Theorem, since the boundary equilibria are all unstable in this case, none of their ω-limit sets intersect int(R 2 + ), and by Theorem 3.3, there are no periodic orbits.
, then E * is globally asymptotically stable with respect to the int(R 2 + ). Proof. Consider the Lyapunov function adapted from the one proposed by Harrison [5]: for all N > 0 and P > 0, with equality if and only if P = P * and (N = N * or N = 0). By the LaSalle Invariance Principle [12], all bounded orbits that remain in int(R 2 + ) converge to the largest invariant subset of the set {(N, P ) ∈ R 2 + V = 0}. By Theorem 3.1, all orbits are bounded, and hence E * is globally asymptotically stable with respect to the int(R 2 + ).
4. Poincaré-Andronov-Hopf bifurcation. So far, we have found two branches of fold bifurcations, C 1 and C 2 , and a cusp point at ( in Section 2 and a transcritical bifurcation that occurs at K 2 = r 1 /a in Section 3. In this section, we discuss Hopf bifurcations, also known as Poincaré-Andronov-Hopf bifurcations, as well as associated codimension two and three bifurcations (see Kuznetsov [11]).
Theorem 4.1. The Jacobian evaluated at E * has a pair of pure imaginary eigenvalues if If for some parameter ρ in the model, (i) and (ii) are satisfied at a value of the parameter, ρ c , and in addition, then there is a Hopf bifurcation as ρ passes through ρ c .
Proof. We showed in (3.4) that the characteristic equation of the Jacobian evaluated at the coexistence equilibrium E * = (N * , P * ) is given by By assumption (i), the constant term is positive. It follows that there is a pair of pure imaginary roots if and only if the coefficient of λ vanishes. Noting that at a coexistence equilibrium, P * = G(N * ) = K 2 1 + γ r2 Φ(N * ) and using (ii), this coefficient is equal to Hence, if (i) and (ii) hold at ρ = ρ c , there is a pair of pure imaginary eigenvalues.
As well, transversality holds as that parameter is varied if d dρ α(ρ c ) = 0 and so there is a Hopf bifurcation. Now, let us find an analytic expression for the Hopf bifurcation curves. The Ncoordinate of a Hopf bifurcation point must satisfy the condition G (N * ) > F (N * ), the two equations in (2.7), and Φ(N * )F (N * )−r 2 F (N * )/K 2 = 0, which is equivalent to Solving for N 3 in (2.7) and substituting the result for N 3 into (4.1), we have The two roots of the quadratic equation are For N ± to be positive real numbers, it requires B < 0 and B 2 − 4 A C ≥ 0. That is, and which is equivalent to Also, both N -and P -coordinates of a Hopf bifurcation point must be real and positive. The P -coordinate of the Hopf bifurcation, P ± = F (N ± ), is given by which is the same as (4.4). When , the terms in the numerator are positive, so that P ± are real and positive.
In order to find curves of Hopf bifurcations with respect to K 1 and K 2 , we substitute N ± for N into (4.1) and solving for K 2 and obtain and ∆ H in (4.6). Thus, we consider H 1 and H 2 where Recall that the Hopf bifurcation point must also satisfy the condition, G (N * ) > F (N * ) that is equivalent to Substituting N ± for N into (4.8) and solving for K 2 yields and ∆ H in (4.6). When all the results are combined, the analytic expressions of the Hopf bifurcation curves are given by and K + 1 is defined in (4.5), ∆ H in (4.6), 1 and 2 in (4.9), and h 1 , h 2 , h 3 , and h 4 in (4.7).
Corollary 4.1. A curve of neutral saddles for system (2.1) is given by where where ψ = (r 2 + γΦ)(ΦG − r 2 − γΦ) and Φ, F, G, and their derivatives are evaluated at N = N * . A Hopf bifurcation is sub-critical if q is positive and super-critical if q is negative.
Proof. See Appendix B.
Lemma 4.1. Let Γ be any periodic orbit of (2.1). Then Proof. (i) Since F (N ) < 0 on the entire periodic orbit in this case, ∆ Γ < 0 where ∆ Γ is defined in Lemma 4.1, and so by the Poincaré criterion (see Coppel [1]), the periodic orbit is orbitally asymptotically stable.
, then the equilibrium lies to the right of the maximum of F (N ) and hence, close to the critical value of the Hopf bifurcation, F (N ) < 0 on the entire periodic orbit. Therefore, again ∆ Γ < 0 where ∆ Γ is defined in Lemma 4.1, and so by the Poincaré criterion (see Coppel [1]), the periodic orbit is orbitally asymptotically stable.
Substituting γ ± for γ back into (2.10), gives the corresponding values of K 1 and K 2 at this codimension three bifurcation: The three parameters, γ (or γ ± ), K 1 , and K 2 , are all positive if and only if r 1 ≥ 9r 2 . In particular, if one selects a = 0.4, b = 8, r 1 = 0.36, and r 2 = 0.04, then there is a codimension three Bogdanov-Takens bifurcation at γ = 1, K 1 = 1/2, and K 2 = 5/6. 5. Numerical studies. In this section, we illustrate our analytical results using phase portraits and one-and two-parameter bifurcation diagrams. We demonstrate numerically that system (1.1) exhibits very interesting dynamics resulting from codimension one bifurcations (Hopf, fold, transcritical, cyclic-fold, and homoclinic bifurcations), codimension two bifurcations (Bautin and Bogdanov-Takens bifurcations), and a codimension three bifurcation in which a Bogdanov-Takens bifurcation occurs at a cusp point. We focus on changing the carrying capacities, K 1 and K 2 , and consider how this impacts the host-parasitoid interactions and we discuss the implications for pest control. Figure 4 is a two-parameter bifurcation diagram in the K 1 K 2 -plane. Because many of the bifurcations occur very close together, we also provide a caricature of this two-parameter bifurcation diagram in Figure 5, where certain regions are expanded. This caricature is surrounded by representative phase portraits for each region labelled 1 to 13 separated by the bifurcation curves, except regions 3 and 13. We separate these two region to stress two dramatically different outcomes as far as the application to pest control is concerned although the dynamics in the region are basically the same. In Figure 5, we display representative phase portraits of the system that involve a homoclinic orbit in Figure 6. We observe a curve of homoclinic bifurcations very close to the fold curve, C 1 defined in (2.9a). A portion of these curves appear to overlap. However, it might be the case that the fold bifurcation occurs just before the homoclinic bifurcation and not at the same time.
In the two-parameter bifurcation diagram, Figure 4, K 1 and K 2 are taken as the bifurcation parameters. The values of the other parameters are fixed at a = 0.4, b = 8, r 1 = 2, r 2 = 0.04, and γ = 10. This two-parameter bifurcation diagram is divided into regions denoted by R(i, j), where i = 3, 4, 5, 6, indicates the number of equilibria and j = 0, 1, 2, indicates the number of limit cycles in each region. We zoom in on two portions of the bifurcation diagram, indicated by rectangles. The upper zoomed in portion shown in (b) includes a Bautin bifurcation, and the lower zoomed in portion shown in (c) includes a cusp point and a Bogdanov-Takens bifurcation (BT). This is just above the other Bautin bifucation and the inset includes a portion of the cyclic fold bifurcation curve emanating from that Bautin bifurcation that is too close to the Hopf bifurcation curve to be visible in (a). In the figure legend for the curve types, H + , H − , FB, TB, and CFB indicate super-critical Hopf, subcritical Hopf, fold, transcritical, and cyclic-fold bifurcations, respectively.  In each region of the bifurcation diagram, the symbol R(i, j) shows the number of equilibria (i = 3, 4, 5, 6) and the number of limit cycles (j = 0, 1, 2). When there is a unique limit cycle, the stability of the limit cycle is indicated by a superscript: S for stable and U for unstable.
Using the analysis of the Hopf bifurcation given in Section 4, the values of the parameters K 1 and K 2 at the two Bautin bifurcation points were determined analytically. Appendix C contains our Maple calculations used to find these values of K 1 and K 2 : (K 1 , K 2 ) ≈ {(0.7619259522, 1.098082881), (3.260434020, 3.800173134)} Next, we consider two one-parameter bifurcation diagrams with K 1 as the bifurcation parameter for different fixed values of K 2 , in order to describe some possible implications for pest control as K 1 increases from 0. We keep all of the other parameters the same as in the two-parameter bifurcation diagram shown in Figure 4. In the first one-parameter bifurcation diagram, shown in Figure 7(e), with representative phase portraits shown in (a)-(d), we fixed K 2 = 1.14408169, and in the second In Figure 7, we consider the sequence of bifurcations as the carrying capacity of the leaf miner population increases from 0, with the carrying capacity of the predator (parasitoid) population fixed at K 2 = 1.14408169. This is a little bit above the lower Bautin bifurcation, but below the cusp point (see Figures 4 and 5). We showed analytically that those points occur at (K 1 , K 2 ) = (0.7619259522, 1.098082881) and (K 1 , K 2 ) = (9.16, 1.15909098381), respectively, (see (2.10) and Appendix C). As K 1 increases form 0, the dynamics change from those shown in region 3 to those shown in 7 to 12 to 13. Figure 7(a) demonstrates that when the carrying capacity of the leaf miner population is small, i.e., to the left of the super-critical Hopf bifurcation, all solutions converge to the asymptotically stable coexistence equilibrium, at which there is a relatively small population of leaf miners. Hence, in this case, the population of leaf miners is fairly well controlled by the parasitoids, although the infestation of the pest is not entirely eradicated. When K 1 = 0.6, a value of K 1 between the super-and sub-critical Hopf bifurcations (see Figure 7(b)), there is an orbitally asymptotically stable periodic orbit with a relatively large amplitude. In this case, the size of the pest outbreak fluctuates periodically from very large to very small. If one is aware that this is the case, other measures should be taken to try to eradicate the outbreak when the leaf miner population becomes small to try to avoid another inevitably large outbreak if nothing else is done. When K 1 = 0.871 is even larger as in Figure 7(c)), so that it is between the sub-critical Hopf bifurcation and the cyclic fold bifurcation, there is bistability. Besides the large amplitude orbitally asymptotically stable periodic orbit, there is an unstable periodic orbit and  Figure 4. Equilibria can be found at the intersection of the host and parasitoid nullclines. Filled circles indicate stable equilibria and open circles indicate unstable equilibria. SLC and ULC indicate stable and unstable limit cycles, respectively. (e) One-parameter bifurcation diagram with K 1 as the bifurcation parameter. Branches of stable equilibria are drawn using thick solid curves and unstable equilibria are drawn using thin solid curves. The largest and smallest values of the N -coordinate of the orbitally asymptotically stable limit cycles are shown using filled circles and of the unstable limit cycles using open circles. H + and H − indicate super-and sub-critical Hopf bifurcations, respectively. CFB indicates the cyclic-fold bifurcation (or saddle-node bifurcation of limit cycles).
an asymptotically stable equilibrium point inside the periodic orbits. If there are measures that might be successful at eradication of the leaf miners when the leaf miner population size gets very small, then efforts to reduce K 1 might be helpful as well, but care must be taken not to do anything to move the population size so that it is in the basin of attraction of the stable coexistence equilibrium. If not possible, then rather than having large outbreaks periodically, it might be better to try to increase K 1 slightly so that the two periodic orbits disappear in the cyclic fold bifurcation and the coexistence equilibrium is asymptotically stable, as depicted in Figure 7(d), where K 1 = 0.9 is to the right of the cyclic fold bifurcation. However, this is not ideal, since the size of the leaf miner population is still substantial.
Next we consider a different sequence of bifurcations that occurs when the carrying capacity of the predator (parasitoid) population K 2 is increased to K 2 = 1.3, a value of K 2 above the Bogdanov-Takens bifurcation. See Figure 8. All of the other parameters are the same as in Figures 4 and 7. As the carrying capacity of the leaf miner population, K 1 , increases from 0, as can be seen more easily in Figure 5, just above the Bogdanov-Takens bifurcation, the dynamics move from region 3 to region 7 to region 8 or possibly from region 3 to region 7 to very briefly in region 11 to region 8. In Figure 8, we see that as K 1 initially increases from 0, the only coexistence equilibrium is asymptotically stable. It then loses stability through a super-critical Hopf bifurcation giving birth to an orbitally asymptotically stable periodic orbit. This is followed by a saddle-node bifurcation of equilibria in which an unstable node and a saddle appear and a homoclinic bifurcation in which the periodic orbit disappears. They appear to occur for the same value of K 1 (but the saddle-node bifurcation might only be very close to, but just before the homoclinic bifurcation). This finally results in the stabilization of one of the three coexistence equilibria, the one with the largest pest population. The two unstable coexistence equilibria eventually disappear in a second saddle-node bifurcation. After the first saddle-node bifurcation, the pest population increases as K 1 increases. Representative phase portraits are given in Figure 8(a)-(d) showing the dynamics between the different bifurcations.
The implications for pest control in the second case with K 2 = 1.3 instead of K 2 = 1.14408, are very similar to those for the previous case, even though the sequence of bifurcations is quite different. In particular, when K 1 is very small, then there is an outbreak, but it is relatively small and would be difficult to eradicate completely. As K 1 increases, there is a periodic orbit that gets very close to the P -axis where the pest population is very small, and some other measure such as the application of pesticide or some effect due to climate change or stochasticity, might have a chance to eradicate the pest entirely. Increasing K 1 further, eventually results in a globally asymptotically stable coexistence equilibrium with a relatively large pest population. One difference in the dynamics is that in the former case, there is a range of K 1 , though small, in which there is bistability, i.e. a stable limit cycle and a stable coexistence equilibrium, that does not occur in the latter case. One significant difference with regard to pest control is that even though the difference in K 2 is quite small in these two cases, the range of the parameter K 1 in which there is a large amplitude periodic orbit is larger in the case that K 2 is larger. Therefore, if other measures can be taken to completely eradicate the pest when the pest population gets very small (due to the oscillatory nature of the dynamics), increasing the carrying capacity of the parasitoid would be advisable. But if these measures are not successful, then there is likely to be very larger outbreaks periodically, unless K 2 can be increased above its value where regions 3, 4, and 6 intersect (see Figure 5), and this might make the problem worse.
6. Discussion. In this paper, we revisit the non-spatial host-parasitoid model, system (1.1), proposed by Magal et al. [13] that was motivated by the invasion of Europe by horse-chestnut leaf miners. We give a more detailed analysis of the model and revise their criteria for control strategies of the leaf miner population.
System (1.1) always has three boundary equilibria and up to three coexistence (positive) equilibria. We obtained analytical conditions that divide the K 1 K 2 -plane into regions in which there are no coexistence equilibria, one, two, or three. We considered the local and global stability of these equilibria in Section 3. We then showed that the model displays interesting dynamical behavior using a bifurcation theory approach. In Section 4, criteria for the existence of Hopf bifurcations were obtained and an analytical expression for determining the criticality of the Hopf bifurcations and the two Bautin bifurcations were provided. Analytic expressions for the transcritical bifurcation, the two curves of fold bifurcations, and the cusp point where the fold curves meet in the K 1 K 2 -plane were also provided. In Section 5, a two-parameter bifurcation diagram was given providing numerical evidence for cyclic fold bifurcations, homoclinic bifurcations, and codimension two and three Bogdanov-Takens bifurcations. Different sequences of bifurcations and their implications for pest control were also discussed in Section 5, focussing on when the carrying capacity of the parasitoid population is kept fixed and the carrying capacity of the pest population is increased from 0. In [13], the authors assumeed that the outbreak can be controlled if there is no coexistence equilibrium. In this case the host-free equilibrium on the P -axis is globally asymptotically stable (see Figure 3). However, their criteria characterizing this case mistakenly includes a small region in which there are two coexistence equilibria (see Figure 9 and the hatched region in Figure 10). Rather than equating the nullclines of the host and parasitoid equations, F (N ) = G(N ), and determining the branches of saddle-node (or fold) and transcritical bifurcations, as we did, in [13] they considered three cases based on the relative values of K 2 and r 1 /a and of F (N ) and G(N ) at the maximum value of F (N ), which occurs at N M = (K 1 −1/b)/2 (see Figure 1). In the case that K 2 > r 1 /a and F (N M ) < G(N M ), they claimed that no coexistence equilibria are possible, and hence the parasitoid population eventually eradicates the leaf miner population. However, they overlooked the possibility shown in Figure 9 when two coexistence equilibria are possible in this case. In particular, they claimed that there are no coexistence equilibria if K 2 > r 1 /a, which is the value of K 2 at the transicritical bifurcation, and either K 1 < 1/b or K 1 > 1/b  Figure 10. Discrepancy (hatched region) between the graphs of A = B in [13] and the portion of the fold curve, C 1 , above the curve of transcritical bifurcations (TB) in the two-parameter bifurcation diagram. The parameter values are the same as the ones used in Figure 9. Magal et al. [13] claim that there are no interior equilibria in the hatched region. However, there are actually two interior equilibria in that region.
and A ≤ B, where A = r 1 (bK 1 + 1) 2 4abK 1 and However, as shown in Figure 9, two coexistence equilibria exist, even though K 2 > r 1 /a, A < B, and K 1 > 1/b. This demonstrates that the criteria based on A = B, or equivalently overestimates the region of initial condition independent eradication of the pest. It mistakenly includes the hatched region above the transcritical bifurcation shown in Figure 10. See Figure 5 for the representative dynamics in phase space for region 2. Hence, although eradication is possible, the outcome would depend on the initial conditions. Therefore, the curve A = B must be replaced by the portion of the fold curve C 1 above the transcritical bifurcation.
In [13], the authors only considered what they called control and possible control of the pest to be possible when the carrying capacity of the parasitoid is above its value at the transcritical bifurcation, i.e., K 2 ≥ r 1 /a. We also considered the implications for pest reduction, even when K 2 < r 1 /a, since it might not be possible to find a parasitoid population with K 2 ≥ r 1 /a or the parasitoid population that is already established might not satisfy this property. Besides trying to manipulate the environment to move the carrying capacities to region 1 in Figure 5 where the model predicts eventual initial condition independent eradication of the pest, or region 2 where there is potential for complete eradication based on initial conditions, our two-parameter bifurcation diagram given in Figures 4 and its caricature given in Figure 5 allowed us to consider other possible pest reduction strategies.
First we considered the case in which the host population had invaded and had established itself at or near its carrying capacity. Our analysis indicates that introduction of the parasitoid population always reduces the pest population below its carrying capacity and is therefore always somewhat effective, provided the parasitoid population itself is not a nuisance. However, the effect might be small, as in region 13, (as a opposed to region 3, where the eventual size of the host population would be relatively small, and this is why we labelled these regions differently even though the dynamics are the same, i.e., convergence to the unique coexistence equilibrium). In regions 2, 4, 6, and 8, where K 1 is fairly large, the outcome would be initial condition dependent, and since the pest is already established at its carrying capacity, even introducing a large number of parasitoids would probably still result in the host population surviving at a new steady state with a relatively large size. If, however, K 1 is small, as in regions 5, 7, 9, 10, 11, and 12, substantial outbreaks would be followed by periods of time in which the pest population gets very small, and the introduction of some other measures to try to totally eradicate the pest population at that time might work, especially in the light of possible stochastic effects.
Next we considered the scenario in which the parasitoid population is established at its carrying capacity, and a small number of the pest population tries to invade. They would likely be successful and able to establish themselves at levels close to their carrying capacity if the carrying capacities are in regions 6 and 8, but the dynamics would again be oscillatory if they invaded with carrying capacity parameters in regions 5, 7, 9, 10, 11, and 12 with the same potential for eradication using other measures.
Appendix A. The number of coexistence equilibria in different regions of the K 1 K 2 -plane shown in Figure 2. The number of positive equilibria is determined by the number of positive roots of the cubic equation, (A.2) In this section, we denote the corresponding positive equilibria by E * = (N * , P * ) where N * is a root of equation (A.1). Note that P * is positive for N * > 0, because the parasitoid nullcline G(N ) is always positive for all positive N . That implies that the number of E * is the same as the number of positive roots of the cubic equation. Our ultimate goal in this section is to generate the K 1 K 2 -plane divided by regions corresponding to the number of E * .
The host nullcline F (N ) is a strictly decreasing function for N > 0, when K 1 ≤ 1/b. In this case, there is a unique E * if K 2 ≤ r 1 /a and no coexistence equilibrium if K 2 > r 1 /a. In our following analysis, we restrict our attention to K 1 > 1/b.
Descartes' rule of signs [2] is a useful technique to determine the number of positive real zeros of a polynomial function. Regardless of the actual values of coefficients, the number of positive roots is the same as the number of the sign changes on the coefficients. There are no, exactly one, two or no, or three or one E * if there are no, one, two, or three sign changes, respectively.
The signs of the coefficients in (A.2) are Table A includes the conditions of the number of E * , by means of the Descartes' rule of signs.

(A.4)
Notice that κ is either undefined or negative if γ ≥ br 2 /a. The function Z 1 (K 1 ) is monotonically increasing and concave down. There exist a vertical asymptote at K 1 = 0 and horizontal asymptote at K 2 = 2br 1 r 2 /(a(br 2 + aγ)). Figure 11 shows regions that correspond to the number of E * on the K 1 K 2plane based on the results in Table A. No coexistence equilibria E * occur in the white region and exactly one E * occurs in the light grey region. Unfortunately, the Descartes' rule of signs does not tell us the regions of exactly two or exactly three E * : it only tells us regions of two or no E * in the dark grey regions and three or one E * in the hatched regions. For that reason, we need to further examine the dark grey and hatched regions.
First, we rewrite the cubic equation as and let the left-hand side of the equation be L(N ) and the right-hand side be C. The number of E * is equal to the number of points of intersection of the graphs of One E * Two or no E * Three or one E * Thus, there are one N -intercept at 0 if ∆ R < 0, equivalently, K 2 > Z 2 (K 1 ) where two intercepts (one of the two is a double root) if ∆ R = 0, or three intercepts K 2 < Z 2 (K 1 ). After a little bit of algebra, one can find that N + is positive if K 1 > 2/b or if K 1 < 2/b and K 2 > Z 1 (K 1 ) and N − is positive if K 1 > 2/b and K 2 > Z 1 (K 1 ). Other properties of the function, L(N ), are where the prime represents the derivative with respect to N and ∆ C = r 1 r 2 [r 1 r 2 (bK 1 + 1) 2 − 3aK 1 K 2 (aγ + br 2 )] where There are five possible graphs of L(N ) for N > 0 as shown in Figure 12. However, we neglect Figure 12(a) since we only care about two cases: (i) two or no E * and (ii) three or one E * . Figure 12(b)-(e) shows that a local minimum occurs at N c + := N min and local maximum occurs at N c − := N max if they exist. That implies that a fold bifurcation, at which either two equilibria collide and disappear or two equilibria appear out of the blue, arises when C is equal to L(N ) evaluated either at N min or N max . The expressions of the fold bifurcations as a function of K 1 are where K 1 ≥ 8γa−r2b b(γa+r2b) := θ 1 and s 2 = b 2 (r 2 b + γa) 2 , s 1 = 2b(r 2 b + γa)(r 2 b + 10γa), s 0 = r 2 b(r 2 b + 20γa) − 8γ 2 a 2 .
Note that C 1 is obtained from the local minimum and C 2 from the local maximum of L(N ). A cusp point, the point of intersection of the two fold curves, occurs at If γ ≤ br 2 /(8a), a cusp point does not appear in the interior of the first quadrant. In this case, the vertical asymptote exists on the K 2 -axis. The graph of C 1 (K 1 ) is convex with a local minimum at K 1 = κ. As we increase K 1 , C 2 (K 1 ) monotonically decreases and its graph crosses point (θ 2 , r 1 /a) and approaches the horizontal asymptote, K 2 = 4br 1 r 2 γ/(aγ + br 2 ) 2 (< r 1 /a). If br 2 /(8a) < γ ≤ br 2 /(2a), The graph of C 1 (K 1 ) is convex with a local minimum at the point, (κ, r 1 /a) . The graph of C 2 (K 1 ) decreases starting from the cusp point, and approaches the horizontal asymptote. The cusp point appears above K 2 = r 1 /a. If γ > br 2 /(2a), C 1 (K 1 ) is an increasing function for K 1 ≥ θ 1 . When br 2 /(2a) < γ ≤ br 2 /a, C 2 (K 1 ) is a concave function whose local maximum occurs at the point, (κ, r 1 /a). Its graph, again, approaches the horizontal asymptote as K 1 is increased. Finally, if γ > br 2 /a, C 2 (K 1 ) is monotonically increasing toward the horizontal asymptote.
In the following two subsections, we will find the regions of two and three E * . Following the Descartes' rule of signs, we have two ambiguous regions that correspond to two or zero and three or one E * (see dark grey and hatched regions in Figure 11). Once we determine the regions of exactly two or three E * , the remaining portions would correspond to no or one E * .
A.1. Regions with exactly two coexistence equilibria in the K 1 K 2 -plane. Let P min and P max be the local minimum and maximum values of L(N ), respectively. When we move the graph of C = −C from bottom to top, sweeping through the graph of L(N ) in Figure 12, two E * occur when P min < C < 0 in (b) or P min < C ≤ 0 in (e), which corresponds to r 1 a ≤ K 2 < C 1 (K 1 ) where C 1 (K 1 ), a function of K 1 , is the right-hand side of (A.8a). That is, confining the dark grey regions in Figure 11, we have two E * in the regions corresponding to r 1 /a ≤ K 2 < C 1 (K 1 ).
A.2. Regions with exactly three coexistence equilibria in the K 1 K 2 -plane.
In a similar way, three E * appear when P min < C < P max in Figure 12(c) or 0 < C < P max in Figure 12(d) and (e). That is, The regions of three E * , therefore, are portions of the hatched regions in Figure  12(b) and (c) that satisfy the above conditions. Figure 2 displays the resulting regions that correspond to the number of E * in the K 1 K 2 -plane. In particular, in Figure 2(a), when 0 < γ < br 2 /(2a), the graph of C 2 is not shown since it is irrelevant to separating the regions of the different number of E * .