A comparative analysis of host–parasitoid models with density dependence preceding parasitism

ABSTRACT

We present a systematic comparison and analysis of four discrete-time, host–parasitoid models. For each model, we specify that density-dependent effects occur prior to parasitism in the life cycle of the host. We compare density-dependent growth functions arising from the Beverton–Holt and Ricker maps, as well as parasitism functions assuming either a Poisson or negative binomial distribution for parasitoid attacks. We show that overcompensatory density-dependence leads to period-doubling bifurcations, which may be supercritical or subcritical. Stronger parasitism from the Poisson distribution leads to loss of stability of the coexistence equilibrium through a Neimark–Sacker bifurcation, resulting in population cycles. Our analytic results also revealed dynamics for one of our models that were previously undetected by authors who conducted a numerical investigation. Finally, we emphasize the importance of clearly presenting biological assumptions that are inherent to the structure of a discrete-time model in order to promote communication and broader understanding.

KEYWORDS:

• Host–parasitoid models
• discrete-time models
• bifurcations
• Jury conditions
• stability

2010 MATHEMATICS SUBJECT CLASSIFICATIONS:

• 37N25
• 92D25
• 92D40
• 37G35
• 39A60

1. Introduction

The interactions between insect parasitoids and their hosts are of great interest to ecologists. Roughly 8.5% of insect species are parasitoids [12], and they play a significant role in regulating their hosts. Because parasitoid species are specialists on suitable prey, they are often used in biological control programs. This has fueled much interest in developing a better understanding of the dynamics of parasitoids and their hosts. Mathematical models of these host–parasitoid systems are also notable because of the simple and specific modelling assumptions that result from the direct connection between parasitized hosts and parasitoid offspring.

Nicholson and Bailey [33] laid the foundation for the study of discrete-time host–parasitoid models. Their basic model assumed that oviposition by parasitoids is limited by the number of encounters with hosts and not by parasitoid egg-supply. In addition, they assumed that the number of encounters with hosts is proportional to host abundance and that hosts are equally susceptible to randomly distributed encounters. Their model, however, yields unstable dynamics. As a result, much of the subsequent literature has sought to investigate factors that induce stability.

In a particularly influential paper, Beddington et al. [6] incorporated density-dependent host recruitment, resulting in the model (1a) $\begin{array}{rl}{N}_{t+1}& =N_t{\mathrm{e}}^{r\left(1-N_t/K\right)}{\mathrm{e}}^{-a{P}_t},\end{array}$(1a) (1b) $\begin{array}{rl}{P}_{t+1}& =c{N}_t\left(1-{\mathrm{e}}^{-a{P}_t}\right).\end{array}$(1b) Here, $N_t$ is the host density, $P_t$ is the parasitoid density, r is the intrinsic rate of growth, K is the host carrying capacity, a is the parasitoid searching efficiency or area of discovery, and c is the parasitoid clutch size. For a detailed explanation of searching efficiency, see [33].

Beddington et al. [6] did not specify the life-stage of the host species for which $N_t$ is the density. This is in contrast to Nicholson and Bailey [33], who provide extensive biological detail for their model. Beddington et al.'s model also fails to provide a coherent explanation of when the density dependence and parasitism occur during the life-cycle of the host. Specifying the order of events is critical when both density dependence and parasitism affect the host population.

Model (1(1a) $\begin{array}{rl}{N}_{t+1}& =N_t{\mathrm{e}}^{r\left(1-N_t/K\right)}{\mathrm{e}}^{-a{P}_t},\end{array}$(1a) ) is an example of the more generalized model (2a) $\begin{array}{rl}{N}_{t+1}& =N_{t}g\left(N_t\right)f\left(P_t\right),\end{array}$(2a) (2b) $\begin{array}{rl}{P}_{t+1}& =c{N}_t[1-f(P_t\left)\right].\end{array}$(2b) Model (2(2a) $\begin{array}{rl}{N}_{t+1}& =N_{t}g\left(N_t\right)f\left(P_t\right),\end{array}$(2a) ) assumes that parasitism affects the original $N_t$ hosts, so that a fraction of hosts, $f\left(P_t\right)$, survive parasitism. The survivors then produce offspring with a per capita recruitment, $g\left(N_t\right)$, that depends on the original number of hosts. The model also assumes that new parasitoids are produced in proportion to the number of parasitized hosts. Murdoch et al. [30] note that the host biology described above is unlikely, though May et al. [28] provide the example of the winter moth (Operophtera brumata) and a fly, Cyzenis albicans, that have this biology.

May et al. [28] evaluated model (2(2a) $\begin{array}{rl}{N}_{t+1}& =N_{t}g\left(N_t\right)f\left(P_t\right),\end{array}$(2a) ) along with two other models to investigate whether the temporal sequence of host density-dependence and parasitism can affect the dynamics of the populations. In a conclusion that is consistent with the earlier findings of Wang and Gutierrez [37], May et al. noted that the ‘sequence of density dependence and parasitism in the host life-cycle can have a significant effect on the population dynamics’ [28]. May et al. further recommended that model (2(2a) $\begin{array}{rl}{N}_{t+1}& =N_{t}g\left(N_t\right)f\left(P_t\right),\end{array}$(2a) ) be abandoned unless the biology of a particular system demands it.

Numerous investigators [4,9,10,13,15,16,21,22,25] have nevertheless cited Beddington et al. [6] and use the structure of model (2(2a) $\begin{array}{rl}{N}_{t+1}& =N_{t}g\left(N_t\right)f\left(P_t\right),\end{array}$(2a) ). These authors often derive their models from previous work, without a careful explanation of the underlying biology. Many books [1,8,11,31,36] also present some version of model (2(2a) $\begin{array}{rl}{N}_{t+1}& =N_{t}g\left(N_t\right)f\left(P_t\right),\end{array}$(2a) ). Mills et al. [29], Murdoch et al. [30], and Hassell [14], are among the few authors who recognize and discuss the biological assumptions inherent in model (2(2a) $\begin{array}{rl}{N}_{t+1}& =N_{t}g\left(N_t\right)f\left(P_t\right),\end{array}$(2a) ).

In this paper, we carefully develop, analyze, and compare four models that assume that density-dependent growth precedes parasitism. These models correspond to a biologically reasonable alternative system presented by May et al. [28]. We consider two functions for density dependence of the host and two functions for parasitism. For each combination of these nonlinear functions, we perform stability analyzes to determine dynamics and bifurcations. From these analyzes, we conclude that stronger nonlinearity in the density-dependence term produces different effects than stronger parasitism.

This paper has eight sections. In the second section, we present the biological assumptions underlying our models and the general form of our equations. In the third section, we outline our methods of analysis. In the following four sections, we present four models. The first two models use a fractional function for parasitism, while the next two models use an exponential form. The first and third models assume compensatory density-dependence, while the second and fourth models include overcompensatory density-dependence. The first model yields highly stable dynamics. The second model has a period-doubling route to chaos. A Neimark–Sacker bifurcation occurs in the third model. The fourth model has exponential functions for both density dependence and parasitism, leading to the greatest variability in dynamics. For certain parameter values, there are two interior equilibria, no more than one of which is stable. Both a Neimark–Sacker bifurcation and a subcritical period-doubling bifurcation occur in this model. We conclude with a discussion of the value of understanding the differences between these models.

2. Model formulation

We now consider the model (3a) $\begin{array}{rl}{N}_{t+1}& =N_{t}G\left(N_t\right)F\left(P_t\right),\end{array}$(3a) (3b) $\begin{array}{rl}{P}_{t+1}& =c{N}_{t}G\left(N_t\right)[1-F(P_t\left)\right].\end{array}$(3b) Although this model is consistent with more than one biological scenario, we make several specific choices here. Let $N_t$ be the density of reproducing host adults, and let $P_t$ be the density of adult female parasitoids. $G\left(N_t\right)$ is the host per-capita-recruitment. $F\left(P_t\right)$, in turn, is the fraction of host offspring that escape parasitism, while $1-F\left(P_t\right)$ is the fraction of host offspring that succumb to parasitism.

In order to analyze zero-growth isoclines more easily, we let $H\left(P_t\right)$ be the fraction of hosts that succumb to parasitism per adult female parasitoid, (4) $H\left(P_t\right)=\frac{1-F\left(P_t\right)}{P_t}.$(4) System (3(3a) $\begin{array}{rl}{N}_{t+1}& =N_{t}G\left(N_t\right)F\left(P_t\right),\end{array}$(3a) ) can now be written (5a) $\begin{array}{rl}{N}_{t+1}& =N_{t}G\left(N_t\right)[1-P_{t}H(P_t\left)\right],\end{array}$(5a) (5b) $\begin{array}{rl}{P}_{t+1}& =c{N}_{t}G\left(N_t\right)P_{t}H\left(P_t\right),\end{array}$(5b) where c is the clutch size. More precisely, c is the average number of female parasitoids laid on a single host that emerge and successfully become reproducing adults. This model is consistent with the second formulation discussed by May et al. [28].

Figure 1 illustrates a host life-cycle scenario that matches the biological assumptions of system (5(5a) $\begin{array}{rl}{N}_{t+1}& =N_{t}G\left(N_t\right)[1-P_{t}H(P_t\left)\right],\end{array}$(5a) ). As above, $N_t$ is the density of reproducing host adults. These adults lay eggs that hatch into larvae. The larvae compete for resources, and $N_{t}G\left(N_t\right)$ larvae survive to the end of larval development. The larvae become pupae, some of which are parasitized, leaving $N_{t+1}$ adults in the next generation.

Figure 1. A host life-cycle diagram that illustrates a set of biological assumptions that match the formulation of the model set-up with density-dependent competition preceding parasitism. $N_t$ is the density of viable adult hosts that reproduce, and $P_t$ is the density of adult female parasitoids.

Although the scenario we have described is that of a pupal parasitoid, we emphasize that this is not the only biological scenario described by systems (3(3a) $\begin{array}{rl}{N}_{t+1}& =N_{t}G\left(N_t\right)F\left(P_t\right),\end{array}$(3a) ) and (5(5a) $\begin{array}{rl}{N}_{t+1}& =N_{t}G\left(N_t\right)[1-P_{t}H(P_t\left)\right],\end{array}$(5a) ). The key point, emphasized by Murdoch [30] and Hassell [14], is that this formulation matches a host life-cycle in which density-dependent competition precedes parasitism.

We now return to the model. For host density-dependent recruitment, we compare Beverton–Holt growth, (6) $N_{t}G\left(N_t\right)=\frac{R_0{N}_t}{1+\frac{(R_0-1)}{K}{N}_t},$(6) and the Ricker curve, (7) $N_{t}G\left(N_t\right)=N_t{\mathrm{e}}^{r\left(1-N_t/K\right)}.$(7) Here $R_0=\exp \left(r\right)$ is the net reproductive rate, $r=\ln \left(R_0\right)$ is the intrinsic rate of growth, and K is the carrying capacity. Recall that the Beverton–Holt growth function is compensatory while the Ricker growth function is overcompensatory.

Early investigators [33,35] used the zero term of the Poisson distribution for $F\left(P_t\right)$, the fraction of hosts that escape parasitism. May [27] considered varying levels of aggregation and proposed the use of the zero term of the negative binomial distribution, (8) $F\left(P_t\right)={\left(1+\frac{a{P}_t}{\kappa }\right)}^{-\kappa }.$(8) May's use of this function influenced Livadiotis et al. [26], who studied system (2(2a) $\begin{array}{rl}{N}_{t+1}& =N_{t}g\left(N_t\right)f\left(P_t\right),\end{array}$(2a) ) with κ-parameterized functions for both parasitism and density-dependent intraspecific competition. The formulation used by Livadiotis et al. highlights the similarities in the exponential ($\kappa \to \mathrm{\infty }$) and rational ($\kappa =1$) functions most commonly used for $F\left(P_t\right)$ and $G\left(N_t\right)$.

In this paper, we focus on two forms of May's function, given by $\kappa =1$ and $\kappa \to \mathrm{\infty }$. From Equation (4(4) $H\left(P_t\right)=\frac{1-F\left(P_t\right)}{P_t}.$(4) ), these values of κ give the fraction of hosts that succumb to parasitism per adult female parasitoid as (9) $H\left(P_t\right)=\frac{a}{1+a{P}_t}$(9) and (10) $H\left(P_t\right)=\frac{1}{P_t}\left(1-{\mathrm{e}}^{-a{P}_t}\right)$(10) respectively.

Other than in May et al.'s paper [28], system (5(5a) $\begin{array}{rl}{N}_{t+1}& =N_{t}G\left(N_t\right)[1-P_{t}H(P_t\left)\right],\end{array}$(5a) ) has not been studied in a way that compares functional forms for modelling parasitism and density dependence. Using Equations (6(6) $N_{t}G\left(N_t\right)=\frac{R_0{N}_t}{1+\frac{(R_0-1)}{K}{N}_t},$(6) ), (7(7) $N_{t}G\left(N_t\right)=N_t{\mathrm{e}}^{r\left(1-N_t/K\right)}.$(7) ), (9(9) $H\left(P_t\right)=\frac{a}{1+a{P}_t}$(9) ), and (10(10) $H\left(P_t\right)=\frac{1}{P_t}\left(1-{\mathrm{e}}^{-a{P}_t}\right)$(10) ) with system (5(5a) $\begin{array}{rl}{N}_{t+1}& =N_{t}G\left(N_t\right)[1-P_{t}H(P_t\left)\right],\end{array}$(5a) ), we will formulate four possible models and compare their dynamics in Sections 4–7.

3. Methods of analysis

Each of our four models can be written in the general ‘density-dependence first’ form of system (5(5a) $\begin{array}{rl}{N}_{t+1}& =N_{t}G\left(N_t\right)[1-P_{t}H(P_t\left)\right],\end{array}$(5a) ). We now nondimensionalize. If we let $y_t=a{P}_t$, $x_t=N_t/K$, and b = acK, we obtain (11a) $\begin{array}{rl}{x}_{t+1}& =x_{t}u(x_t,y_t),\end{array}$(11a) (11b) $\begin{array}{rl}{y}_{t+1}& =y_{t}v(x_t,y_t),\end{array}$(11b) where (12) $\begin{array}{rl}u(x_t,y_t)& =g\left(x_t\right)[1-y_{t}h(y_t\left)\right],\end{array}$(12) (13) $\begin{array}{rl}v(x_t,y_t)& =b{x}_{t}g\left(x_t\right)h\left(y_t\right).\end{array}$(13) For Beverton-Holt growth, (14) $g\left(x_t\right)=\frac{R_0}{1+(R_0-1)x_t},$(14) while for Ricker growth, (15) $g\left(x_t\right)={\mathrm{e}}^{r(1-x_t)}.$(15) Both expressions of per-capita recruitment can be found from $G\left(N_t\right)$ by directly substituting $x_t=N_t/K$. We will call (14(14) $g\left(x_t\right)=\frac{R_0}{1+(R_0-1)x_t},$(14) ) fractional per-capita-recruitment, which produces compensatory density-dependence, and (15(15) $g\left(x_t\right)={\mathrm{e}}^{r(1-x_t)}.$(15) ) exponential per-capita-recruitment, which produces overcompensatory density-dependence.

Similarly, the fraction of hosts that succumb to parasitism, $yh\left(y\right)$, can be rewritten with (16) $h\left(y_t\right)=\frac{1}{1+y_t},$(16) for $\kappa =1$, and (17) $h\left(y_t\right)=\frac{1}{y_t}\left(1-{\mathrm{e}}^{-y_t}\right),$(17) for $\kappa \to \mathrm{\infty }$. For both cases, to get $h\left(y_t\right)$ from $H\left(P_t\right)$, it is necessary to divide by a and then substitute $y_t=a{P}_t$. We will refer to (16(16) $h\left(y_t\right)=\frac{1}{1+y_t},$(16) ) as fractional parasitism and (17(17) $h\left(y_t\right)=\frac{1}{y_t}\left(1-{\mathrm{e}}^{-y_t}\right),$(17) ) as exponential parasitism.

In all that follows, we assume $R_0\ge 1$ ($r\ge 0$), since we choose to consider cases where the host species can persist in the absence of the parasitoid species. For $R_0>1$ (r>0), the per-capita recruitment, $g\left(x_t\right)$, is a positive, monotonically decreasing function that starts from $R_0=\ln \left(r\right)$ at $x_t=0$ and crosses 1 at $x_t=1$. Similarly, $h\left(y_t\right)$ is positive and monotonically decreasing, with $h\left(0\right)=1$. Sample plots of $g\left(x\right)$, $xg\left(x\right)$, and $h\left(y\right)$ are shown in Figure 2.

Figure 2. These figures illustrate the behavior of $g\left(x\right)$, $xg\left(x\right)$, and $h\left(y\right)$ for functions used in our models. Fractional forms of $g\left(x\right)$ and $h\left(y\right)$ from Equations (14(14) $g\left(x_t\right)=\frac{R_0}{1+(R_0-1)x_t},$(14) ) and (16(16) $h\left(y_t\right)=\frac{1}{1+y_t},$(16) ) are shown with solid lines. Exponential forms of $g\left(x\right)$ and $h\left(y\right)$ from Equations (15(15) $g\left(x_t\right)={\mathrm{e}}^{r(1-x_t)}.$(15) ) and (17(17) $h\left(y_t\right)=\frac{1}{y_t}\left(1-{\mathrm{e}}^{-y_t}\right),$(17) ) are shown with dashed lines. Both functions for $g\left(x\right)$ are monotonically decreasing from $R_0$. Recruitment, $xg\left(x\right)$, is non-monotonic for the exponential form, while it is monotonic for the fractional form. Both fractional and exponential forms of $h\left(y\right)$ are positive and monotonically decreasing. The dashed curve remains above the solid curve as y increases. (a) Host per-capita-recruitment. (b) Host recruitment and (c) Parasitism.

To find the equilibria of system (11(11a) $\begin{array}{rl}{x}_{t+1}& =x_{t}u(x_t,y_t),\end{array}$(11a) ), we set $x_{t+1}=x_t=x^{\ast }$ and $y_{t+1}=y_t=y^{\ast }$. The equilibria occur at (0,0), (1,0), and at solutions of the system (18a) $\begin{array}{rl}1& =u(x^{\ast },y^{\ast })=g\left(x^{\ast }\right)[1-y^{\ast }h(y^{\ast }\left)\right],\end{array}$(18a) (18b) $\begin{array}{rl}1& =v(x^{\ast },y^{\ast })=b{x}^{\ast }g\left(x^{\ast }\right)h\left(y^{\ast }\right).\end{array}$(18b) For each of our models, it can be shown that b>1 is a necessary and sufficient condition for the existence of a unique positive solution to system (18(18a) $\begin{array}{rl}1& =u(x^{\ast },y^{\ast })=g\left(x^{\ast }\right)[1-y^{\ast }h(y^{\ast }\left)\right],\end{array}$(18a) ). For the fourth model, there is a region below b = 1 for which two positive solutions to system (18(18a) $\begin{array}{rl}1& =u(x^{\ast },y^{\ast })=g\left(x^{\ast }\right)[1-y^{\ast }h(y^{\ast }\left)\right],\end{array}$(18a) ) exist.

To determine the stability of the equilibria, we form the Jacobian matrix of partial derivatives for system (11(11a) $\begin{array}{rl}{x}_{t+1}& =x_{t}u(x_t,y_t),\end{array}$(11a) ), (19) $J(x,y)=\left(\begin{array}{cc}x{u}_x+u& x{u}_y\\ y{v}_x& y{v}_y+v\end{array}\right),$(19) where we drop the t subscripts for notational simplicity. After evaluating the partial derivatives, the Jacobian may be rewritten (20) $J(x,y)=\left(\begin{array}{cc}\left[x{g}^{\prime }\right(x)+g(x\left)\right][1-yh(y\left)\right]& -xg\left(x\right)\left[y{h}^{\prime }\right(y)+h(y\left)\right]\\ byh\left(y\right)\left[x{g}^{\prime }\right(x)+g(x\left)\right]& bxg\left(x\right)\left[y{h}^{\prime }\right(y)+h(y\left)\right]\end{array}\right),$(20) where we factor to separate the x and y dependencies. We will use the Jacobian evaluated at each of equilibria of system (11(11a) $\begin{array}{rl}{x}_{t+1}& =x_{t}u(x_t,y_t),\end{array}$(11a) ) to determine stability.

3.1. Host-only system summary

In the absence of the parasitoid, system (11(11a) $\begin{array}{rl}{x}_{t+1}& =x_{t}u(x_t,y_t),\end{array}$(11a) ) reduces to (21) $x_{t+1}=x_{t}g\left(x_t\right).$(21) For fractional recruitment, the host population goes extinct if $R_0<1$. All values of x are neutrally stable equilibria if $R_0=1$. For $R_0>1$, the host density will approach the asymptotically stable equilibrium, x = 1.

For exponential recruitment, the host population will die out if r<0. Similarly to the fractional recruitment case, all values of x are neutrally stable equilibria if r = 0. The host population will persist if r>0. If 0<r<2, the point x = 1 is an asymptotically stable equilibrium. At r = 2, a supercritical flip (period-doubling) bifurcation occurs. The attractor in the system becomes a two-cycle, which persists as r continues to increase, until another period-doubling occurs, and there is a stable four-cycle. Period doubling continues as r increases, until the system becomes chaotic.

The dynamics of this system are particularly relevant as we now consider the extinction and exclusion equilibria of system (11(11a) $\begin{array}{rl}{x}_{t+1}& =x_{t}u(x_t,y_t),\end{array}$(11a) ). Understanding the period-doubling route to chaos in the host-only system with exponential recruitment is also useful when interpreting the dynamics of the host–parasitoid system with exponential recruitment, especially in Section 5.

3.2. Extinction equilibrium

At the extinction point $(0,0)$, the Jacobian, (22) $\begin{array}{rl}J(0,0)& =\left(\begin{array}{cc}g\left(0\right)& 0\\ 0& 0\end{array}\right)=\left(\begin{array}{cc}{R}_0& 0\\ 0& 0\end{array}\right),\end{array}$(22) has eigenvalues $R_0$ and 0. Note that we used $g\left(0\right)=R_0$, which was mentioned previously. The extinction equilibrium is unstable for $R_0>1$. The zero eigenvalue indicates that for initial conditions with $x=0,y>0$, the system will collapse to the $(0,0)$ fixed point at the next generation due to the lack of hosts.

3.3. Exclusion equilibrium

The equilibrium point $(1,0)$ is known as an exclusion point [4,20,21]. Here, the host population persists at carrying capacity, while the parasitoid population is extinct. The Jacobian for this system is (23) $\begin{array}{rl}J(1,0)& =\left(\begin{array}{cc}{g}^{\prime }\left(1\right)+g\left(1\right)& -h\left(0\right)\\ 0& bh\left(0\right)\end{array}\right)=\left(\begin{array}{cc}{g}^{\prime }\left(1\right)+1& -1\\ 0& b\end{array}\right),\end{array}$(23) since $h\left(0\right)=1$ for Equation (16(16) $h\left(y_t\right)=\frac{1}{1+y_t},$(16) ) and $\underset{y\to 0}{lim}h\left(y\right)=1$ for Equation (17(17) $h\left(y_t\right)=\frac{1}{y_t}\left(1-{\mathrm{e}}^{-y_t}\right),$(17) ). The eigenvalues for this triangular system are thus (24) ${\lambda }_1=g^{\prime }\left(1\right)+1,{\lambda }_2=b.$(24)

Recall that $g^{\prime }\left(1\right)$ is negative since $g\left(x\right)$ is monotone decreasing for $R_0>1$ (r>0). Based on the eigenvalues in (24(24) ${\lambda }_1=g^{\prime }\left(1\right)+1,{\lambda }_2=b.$(24) ), we conclude that the exclusion equilibrium will be asymptotically stable if (25) $-2<g^{\prime }\left(1\right)<0$(25) and b<1. The second inequality in condition (25(25) $-2<g^{\prime }\left(1\right)<0$(25) ) is satisfied, so we will check the first inequality for both forms of the host per-capita-recruitment, $g\left(x\right)$.

For Equation (14(14) $g\left(x_t\right)=\frac{R_0}{1+(R_0-1)x_t},$(14) ), (26) $g^{\prime }\left(1\right)=\frac{(1-R_0)}{R_0},$(26) and the first inequality in (25(25) $-2<g^{\prime }\left(1\right)<0$(25) ) becomes (27) $-2R_0<1-R_0,$(27) which simplifies to $-1<R_0$. Since the net reproductive rate, $R_0$, is positive, this inequality is true, and the stability of the exclusion equilibrium point hinges on the value of b for our models that use fractional recruitment. For b<1, the equilibrium is asymptotically stable, and for b>1, the equilibrium is unstable.

For Equation (15(15) $g\left(x_t\right)={\mathrm{e}}^{r(1-x_t)}.$(15) ), (28) $g^{\prime }\left(1\right)=-r.$(28) Thus, $|{\lambda }_1|<1$ if and only if $-2<-r<0$. For our models that use exponential recruitment, b<1 is not sufficient for asymptotic stability of the exclusion equilibrium. If both b<1 and 0<r<2, then the exclusion equilibrium will be asymptotically stable.

Recall that $R_0\ge 1$ ($r\ge 0$) is necessary for the host to persist when alone. For the host-only system with fractional recruitment, the equilibrium x = 1 is asymptotically stable for $R_0>1$. For the host-only system with exponential recruitment, x = 1 is asymptotically stable if r>0 andr<2. For the host–parasitoid system with either fractional or exponential recruitment, when b<1, a small parasitoid population is unable to increase when the host is at x = 1.

3.4. Coexistence equilibria

The coexistence equilibria are the solutions to system (18(18a) $\begin{array}{rl}1& =u(x^{\ast },y^{\ast })=g\left(x^{\ast }\right)[1-y^{\ast }h(y^{\ast }\left)\right],\end{array}$(18a) ). Biologically, coexistence occurs when both $x^{\ast }$ and $y^{\ast }$ are positive. These equilibria can be explicitly determined for models using fractional parasitism, but not for exponential parasitism. Nevertheless, the coexistence equilibria can be approximated numerically for all cases.

We evaluate Jacobian matrix (19(19) $J(x,y)=\left(\begin{array}{cc}x{u}_x+u& x{u}_y\\ y{v}_x& y{v}_y+v\end{array}\right),$(19) ) at the coexistence equilibrium, $(x^{\ast },y^{\ast })$, and use Equations (18a(18a) $\begin{array}{rl}1& =u(x^{\ast },y^{\ast })=g\left(x^{\ast }\right)[1-y^{\ast }h(y^{\ast }\left)\right],\end{array}$(18a) ) and (18b(18b) $\begin{array}{rl}1& =v(x^{\ast },y^{\ast })=b{x}^{\ast }g\left(x^{\ast }\right)h\left(y^{\ast }\right).\end{array}$(18b) ), which are satisfied at coexistence equilibria, to simplify. We obtain (29) $J(x^{\ast },y^{\ast })=\left(\begin{array}{cc}{x}^{\ast }{u}_x+1& x^{\ast }{u}_y\\ y^{\ast }{v}_x& y^{\ast }{v}_y+1\end{array}\right).$(29) To avoid unnecessarily complicated algebra, we will not proceed from eigenvalues.

We will instead determine sufficient conditions for asymptotic stability of the coexistence equilibria by applying the Jury conditions [18] to each model. The Jury conditions are necessary and sufficient conditions for asymptotic stability in the linearized system. As such, they provide sufficient conditions for asymptotic stability in the nonlinear system. The Jury conditions can be written (30) $\begin{array}{rl}& 1-\tau +\mathrm{\Delta }>0,\end{array}$(30) (31) $\begin{array}{rl}& 1+\tau +\mathrm{\Delta }>0,\end{array}$(31) (32) $\begin{array}{rl}& \mathrm{\Delta }<1,\end{array}$(32) where τ is the trace and Δ is the determinant of the Jacobian matrix evaluated at the implicit or explicit coexistence equilibrium. In τ–Δ space, the intersection of the regions defined by the Jury conditions is a triangle, shown in [1,8] and described by [38], where both eigenvalues are less than one in modulus.

Recall that the stability of a hyperbolic equilibrium is accurately determined by analyzing its stability in the linearized system. For our system with $R_0$ (or r) and b parameters such that $(\tau ,\mathrm{\Delta })$ is inside or outside the triangle, the coexistence fixed point is hyperbolic. Outside the triangle, at least one eigenvalue exceeds one in magnitude, so the fixed point is unstable in the nonlinear system.

On the boundary of the triangle, we cannot use the linearized system to determine the stability of the coexistence fixed point in the nonlinear system. At least one eigenvalue is exactly one in magnitude on the triangle boundaries, meaning that the fixed point is non-hyperbolic, and nonlinear terms are needed to determine stability. However, we will not concern ourselves with stability of non-hyperbolic points since they generally coincide with bifurcations in our models.

Losing stability by violating the first Jury condition, inequality (30(30) $\begin{array}{rl}& 1-\tau +\mathrm{\Delta }>0,\end{array}$(30) ), corresponds to the dominant eigenvalue leaving the unit circle through positive one. If the system instead loses stability when inequality (31(31) $\begin{array}{rl}& 1+\tau +\mathrm{\Delta }>0,\end{array}$(31) ) is first violated, the dominant eigenvalue leaves the unit circle through negative one, causing a flip bifurcation. In both of these cases, the non-dominant eigenvalue remains within the unit circle when the bifurcation occurs. If both the first and second Jury conditions are satisfied and stability is lost through violating the third Jury condition, the eigenvalues are a complex-conjugate pair that leave the unit circle simultaneously.

For a system with τ and Δ that violate two of the three Jury conditions, both eigenvalues will be larger than one in magnitude (or equal to one in magnitude on the corners of the triangle). For further details on the Jury conditions and bifurcations, see Kuznetsov [23] and Whitley [38], though the Jury conditions are not referenced by name in these sources.

For matrix (29(29) $J(x^{\ast },y^{\ast })=\left(\begin{array}{cc}{x}^{\ast }{u}_x+1& x^{\ast }{u}_y\\ y^{\ast }{v}_x& y^{\ast }{v}_y+1\end{array}\right).$(29) ), (33) $\tau =2+x^{\ast }{u}_x+y^{\ast }{v}_y,$(33) and (34) $\mathrm{\Delta }=1+x^{\ast }{u}_x+y^{\ast }{v}_y+x^{\ast }{y}^{\ast }(u_x{v}_y-u_y{v}_x).$(34) Using these expressions, the first Jury condition, inequality (30(30) $\begin{array}{rl}& 1-\tau +\mathrm{\Delta }>0,\end{array}$(30) ), simplifies to (35) $x^{\ast }{y}^{\ast }(u_x{v}_y-u_y{v}_x)>0.$(35) The first Jury condition will be violated for parameter values such that $x^{\ast }=0$ or $y^{\ast }=0$. For a true coexistence equilibrium point with positive $x^{\ast }$ and $y^{\ast }$ values, inequality (35(35) $x^{\ast }{y}^{\ast }(u_x{v}_y-u_y{v}_x)>0.$(35) ) requires (36) $u_x{v}_y-u_y{v}_x>0.$(36)

We now consider the second Jury condition, inequality (31(31) $\begin{array}{rl}& 1+\tau +\mathrm{\Delta }>0,\end{array}$(31) ). After we write the inequality in terms of $u,v,x^{\ast },$ and $y^{\ast }$, the condition simplifies to (37) $4+2x^{\ast }{u}_x+2y^{\ast }{v}_y+x^{\ast }{y}^{\ast }(u_x{v}_y-u_y{v}_x)>0.$(37)

Finally, the third Jury condition, inequality (32(32) $\begin{array}{rl}& \mathrm{\Delta }<1,\end{array}$(32) ), can be expressed as (38) $1+x^{\ast }{u}_x+y^{\ast }{v}_y+x^{\ast }{y}^{\ast }(u_x{v}_y-u_y{v}_x)<1.$(38) These three Jury conditions, inequalities (35(35) $x^{\ast }{y}^{\ast }(u_x{v}_y-u_y{v}_x)>0.$(35) ), (37(37) $4+2x^{\ast }{u}_x+2y^{\ast }{v}_y+x^{\ast }{y}^{\ast }(u_x{v}_y-u_y{v}_x)>0.$(37) ), and (38(38) $1+x^{\ast }{u}_x+y^{\ast }{v}_y+x^{\ast }{y}^{\ast }(u_x{v}_y-u_y{v}_x)<1.$(38) ), will be used for each specific model to determine the requirements on parameters b and $R_0$ (or r) to ensure that the coexistence equilibrium is stable.

4. Model 1: compensatory host density-dependence and fractional parasitism

The first model we consider uses fractional per-capita-recruitment (14(14) $g\left(x_t\right)=\frac{R_0}{1+(R_0-1)x_t},$(14) ) and fractional parasitism (16(16) $h\left(y_t\right)=\frac{1}{1+y_t},$(16) ). The model is thus (39a) $\begin{array}{rl}{x}_{t+1}& =\left[\frac{R_0{x}_t}{1+(R_0-1)x_t}\right]\left(\frac{1}{1+y_t}\right),\end{array}$(39a) (39b) $\begin{array}{rl}{y}_{t+1}& =b\left[\frac{R_0{x}_t}{1+(R_0-1)x_t}\right]\left(\frac{y_t}{1+y_t}\right).\end{array}$(39b)

The coexistence equilibrium for this system is (40) $(x^{\ast },y^{\ast })=\left(\frac{1}{b},\frac{R_0}{1+\left(R_0-1\right)\left(\frac{1}{b}\right)}-1\right).$(40) As shown in Appendix A.1, for $R_0>1$, this equilibrium is in the interior of the first quadrant if and only if b>1. For b = 1, the equilibrium given by Equation (40(40) $(x^{\ast },y^{\ast })=\left(\frac{1}{b},\frac{R_0}{1+\left(R_0-1\right)\left(\frac{1}{b}\right)}-1\right).$(40) ) is the exclusion equilibrium, $(1,0)$. For $R_0=1$, system (39(39a) $\begin{array}{rl}{x}_{t+1}& =\left[\frac{R_0{x}_t}{1+(R_0-1)x_t}\right]\left(\frac{1}{1+y_t}\right),\end{array}$(39a) ) has a line of equilibria on the x-axis, and (40(40) $(x^{\ast },y^{\ast })=\left(\frac{1}{b},\frac{R_0}{1+\left(R_0-1\right)\left(\frac{1}{b}\right)}-1\right).$(40) ) reduces to $(1/b,0)$.

4.1. Stability region

Compensatory (fractional) host recruitment, $xg\left(x\right)$, and fractional parasitism are both rational functions, which correspond to a low κ index in the parameterized families of common recruitment and parasitism functions (see Livadiotis et al. [26]). When we use fractional per-capita-recruitment for $g\left(x\right)$ and fractional parasitism for $h\left(y\right)$, the model has a large stability region as seen in Figure 3(a). All three Jury conditions are satisfied for the region in parameter space defined by b>1, $R_0>1$. The first Jury condition is violated for b = 1. Crossing this line corresponds to a transcritical bifurcation. Both the first and third Jury conditions are violated for $R_0=1$. Details are given in Appendix A.2–A.4. Satisfying the three Jury conditions ensures that the coexistence equilibrium is asymptotically stable.

Figure 3. Stability regions for the positive coexistence equilibrium for Models 1–4. For $R_0>1$, the condition b>1 ensures that the coexistence equilibrium is in the interior of the first quadrant for Models 1–3. For Model 4, the system has a unique positive coexistence equilibrium for b>1. For $R_0>{\mathrm{e}}^2$, there is a region below b = 1 for which there are two positive coexistence equilibria. The coexistence equilibrium with the larger y coordinate is asymptotically stable in the shaded region below b = 1. The vertical solid lines are boundary curves where both the first and third Jury conditions are violated. The horizontal solid lines are boundary curves where the first Jury condition is violated. The dotted curves are the boundaries where the second Jury condition is first violated. The dashed curves are the boundaries where the third Jury condition is first violated. (a) Model 1. (b) Model 2. (c) Model 3 and (d) Model 4.

5. Model 2: overcompensatory host density-dependence and fractional parasitism

Our second model also uses fractional parasitism, but it incorporates the exponential per-capita-recruitment from Equation (15(15) $g\left(x_t\right)={\mathrm{e}}^{r(1-x_t)}.$(15) ). As seen in Figure 2(b), exponential recruitment is nonmonotonic, so we have introduced stronger nonlinearity in the density-dependence term. These choices yield the model (41a) $\begin{array}{rl}{x}_{t+1}& =x_t{\mathrm{e}}^{r(1-x_t)}\left(\frac{1}{1+y_t}\right),\end{array}$(41a) (41b) $\begin{array}{rl}{y}_{t+1}& =b{x}_t{\mathrm{e}}^{r(1-x_t)}\left(\frac{y_t}{1+y_t}\right).\end{array}$(41b)

The coexistence equilibrium for this system, (42) $(x^{\ast },y^{\ast })=\left(\frac{1}{b},{\mathrm{e}}^{r\left(1-1/b\right)}-1\right),$(42) is again in the interior of the first quadrant if r>0, b>1. This is shown in Section A.5. For b = 1, the equilibrium given by Equation (42(42) $(x^{\ast },y^{\ast })=\left(\frac{1}{b},{\mathrm{e}}^{r\left(1-1/b\right)}-1\right),$(42) ) is the exclusion equilibrium, $(1,0)$. For r = 0, system (41(41a) $\begin{array}{rl}{x}_{t+1}& =x_t{\mathrm{e}}^{r(1-x_t)}\left(\frac{1}{1+y_t}\right),\end{array}$(41a) ) has a line of equilibria on the x-axis, and (42(42) $(x^{\ast },y^{\ast })=\left(\frac{1}{b},{\mathrm{e}}^{r\left(1-1/b\right)}-1\right),$(42) ) reduces to $(1/b,0)$.

5.1. Stability region

As was true for Model 1, the first Jury condition is satisfied for b>1, r>0 ($R_0>1$), and the third Jury condition is satisfied for r>0 ($R_0>1$). Substituting the exponential form of density dependence in place of the fractional form from Model 1 introduces an additional stability criterion for the coexistence equilibrium for Model 2. The second Jury condition is now satisfied above the curve defined, for $3/2<\omega <2$, by (43) $r=\omega -\ln \left(2\omega -3\right),b=1-\frac{1}{\omega }\ln (2\omega -3).$(43) The derivation of these stability criteria is shown in Sections A.6–A.8.

The stability region for the coexistence equilibrium is shown in Figure 3(b). Crossing the line segment, b = 1, $1<R_0<{\mathrm{e}}^2$ violates the first Jury condition, resulting in a transcritical bifurcation. Crossing the line $R_0=1$ violates both the first and third Jury conditions. Crossing the dotted curve from the left in Figure 3(b) means that one of the real eigenvalues exceeds $-1$ in magnitude. This corresponds to a period-doubling or flip bifurcation. However, the stability analysis holds only in a neighborhood of the equilibrium point. We discuss below the existence of other stable phenomena for this model, including 2-cycles, 4-cycles, and invariant circles.

5.2. Bifurcations and attractors

For a fixed value of b as r increases, another attractor emerges. For certain values of b, there is a range of r values for which bistability is observed. We describe the behavior for various fixed b as r is increased in the specified range noting that $R_0=\exp \left(r\right)$. This is illustrated in Figure 4.

Figure 4. Bifurcation diagrams for Model 2 for fixed b as r increases. Figures on the left show x coordinates of stable (solid) and unstable (dashed) fixed points and cycles as r increases. Figures on the right show y coordinates of stable (solid) and unstable (dashed) fixed points and cycles as r increases. Detailed descriptions of the dynamics and bifurcations are given in the text. (a) $b=\mathrm{1.1.}$ (b) $b=\mathrm{1.1.}$ (c) $b=\mathrm{1.3.}$ (d) b = 1.3. (e) b = 1.5 and (f) b = 1.5.

• b = 1.004, 1.1, 1.2

  For r<2, the coexistence equilibrium is stable. As r increases beyond r = 2, the system initially has a stable interior equilibrium, an unstable equilibrium at $(1,0)$, and an unstable 2-cycle on the x-axis. (The alternating x values of the 2-cycle are shown for changing r values as dashed lines in Figure 4(a) with y = 0 as seen in Figure 4(b)). The extinction equilibrium and the two-cycle on the axis are unstable because y can increase when rare, such that the system goes to the interior equilibrium. As r increases further, the interior equilibrium undergoes a supercritical flip bifurcation giving rise to a stable 2-cycle. As r continues to increase, the 2-cycle moves towards the x-axis (the values of y go towards 0) before colliding with the unstable 2-cycle and exchanging stability as it passes into the fourth quadrant.

  To summarize, for r somewhat greater than 2, the host-only system has a cycle with period 2 that the parasitoid can invade, stabilizing host dynamics. Increasing r eventually causes the host and parasitoid dynamics to become periodic, but further increasing r decreases the parasitoid density until it cannot persist. After this, the system goes back to the host-only dynamics (a 2-cycle), which are now stable. In this case and the two below, the dashed lines at the left of the x plots that become solid when crossing into the gray shaded region represent the dynamics of the host when it is alone (see Figure 4(a,c,e)).

• b = 1.3, 1.4

  For r sufficiently high, the system has a stable interior equilibrium, an unstable equilibrium at $(1,0)$, and an unstable 2-cycle on the x-axis. As r increases further, a stable 2-cycle emerges with an accompanying unstable 2-cycle in a saddle-node bifurcation of the second iterate of the mapping. Shortly thereafter, the interior equilibrium undergoes a subcritical flip bifurcation when the unstable 2-cycle in the first quadrant collides with it, and the equilibrium loses stability. For further discussion of subcritical flip bifurcations, see [32,38]. As r continues to increase, the stable 2-cycle moves towards the x-axis before colliding with the unstable 2-cycle on the x-axis and exchanging stability as it passes into the fourth quadrant.

• b = 1.5

  For r sufficiently high, the system has a stable interior equilibrium, an unstable equilibrium at $(1,0)$, and an unstable 2-cycle on the x-axis. As r increases further, we first observe that the unstable two-cycle on the x-axis period doubles into a four-cycle. Then, a stable 2-cycle emerges in the interior of the first quadrant with an accompanying unstable 2-cycle in a saddle-node bifurcation of the second iterate of the mapping. Then, the stable 2-cycle undergoes a period doubling bifurcation such that a stable 4-cycle emerges. Shortly thereafter, the interior equilibrium undergoes a subcritical flip bifurcation when the unstable 2-cycle in the interior of the quadrant collides with it. After this bifurcation, the coexistence equilibrium is unstable. As r continues to increase, the stable 4-cycle moves towards the x-axis before colliding with the unstable 4-cycle and exchanging stability as it passes into the fourth quadrant.

It is evident that for higher values of b, the bifurcations associated with increasing r are more complicated. Indeed, for b = 1.9, r = 2.92, the system has an attractor with fractal dimension. A small increase in r to r = 2.9205 results in an attractor made up of four circles such that the union of the four circles is an invariant attractor. In both cases, the equilibrium point is locally stable with its own basin of attraction. These two cases are shown in Figure 5. Further increases in r result in a 4-cycle, which then period doubles into an 8-cycle.

Figure 5. Illustration of bistability between the equilibrium and another attractor for Model 2. The parameters for the left figure are b = 1.9, r = 2.92. For the right figure, $b=1.9,r=2.9205$. For the figure on the left, the attractor is a region with fractal dimension. The attractor on the right consists of four circles such that the union of the circles is an invariant attracting set. If we continue to increase r past r = 2.9205, we see a 4-cycle that then period doubles to an 8-cycle. Initial conditions for both figures were $(0.8,0.7)$ for the equilibrium and $(0.3,0.4)$ for the other attractor. For clarity of the attractors, we ran $\mathrm{100,000}$ iterations and plotted 30, 000 points for the attractors.

6. Model 3: compensatory host density-dependence and exponential parasitism

For the third model under consideration, we return to fractional recruitment, Equation (14(14) $g\left(x_t\right)=\frac{R_0}{1+(R_0-1)x_t},$(14) ), and now incorporate a stronger parasitism term. That is, we now take the limit as $\kappa \to \mathrm{\infty }$ in Equation (8(8) $F\left(P_t\right)={\left(1+\frac{a{P}_t}{\kappa }\right)}^{-\kappa }.$(8) ), which results in exponential parasitism seen in Equation (10(10) $H\left(P_t\right)=\frac{1}{P_t}\left(1-{\mathrm{e}}^{-a{P}_t}\right)$(10) ). Biologically, higher κ corresponds to higher parasitoid aggregation, detailed in [27].

The third model is (44a) $\begin{array}{rl}{x}_{t+1}& =\left[\frac{R_0{x}_t}{1+(R_0-1)x_t}\right]{\mathrm{e}}^{-y_t},\end{array}$(44a) (44b) $\begin{array}{rl}{y}_{t+1}& =b\left[\frac{R_0{x}_t}{1+(R_0-1)x_t}\right]\left(1-{\mathrm{e}}^{-y_t}\right).\end{array}$(44b) As with the other models, the coexistence equilibrium is in the interior of the first quadrant for $R_0>1$, b>1. This is shown in Section A.9. For b = 1, the coexistence equilibrium has collided with the exclusion equilibrium at $(1,0)$. For $R_0=1$, system (44(44a) $\begin{array}{rl}{x}_{t+1}& =\left[\frac{R_0{x}_t}{1+(R_0-1)x_t}\right]{\mathrm{e}}^{-y_t},\end{array}$(44a) ) has a line of equilibria on the x-axis. As mentioned in Section 3.4, we cannot derive an explicit expression for the coexistence equilibrium for models with exponential parasitism.

6.1. Stability region

Even without an explicit expression for the coexistence equilibrium, we can determine the stability criteria. The first Jury condition is satisfied for $R_0>1$, b>1. Satisfying the first Jury condition is a sufficient condition for satisfying the second Jury condition. The third Jury condition, in turn, is satisfied in the $R_0$–b plane below the curve (45) $R_0=\frac{y{\mathrm{e}}^{2y}}{{\mathrm{e}}^y-1},b=\frac{y^2{\mathrm{e}}^{2y}-y{\mathrm{e}}^y+y}{({\mathrm{e}}^y-1)(y{\mathrm{e}}^y-{\mathrm{e}}^y+1)},$(45) for positive y.

We determined the parametric boundary curve for the third Jury condition, inequality (32(32) $\begin{array}{rl}& \mathrm{\Delta }<1,\end{array}$(32) ), by solving the three equations (46) $\begin{array}{rl}\mathrm{\Delta }& =1,\end{array}$(46) (47) $\begin{array}{rl}u(x,y)=g\left(x\right)[1-yh(y\left)\right]& =1,\end{array}$(47) (48) $\begin{array}{rl}v(x,y)=bxg\left(x\right)h\left(y\right)& =1,\end{array}$(48) to eliminate x and write b and $R_0$ as functions of y. Equations (47(47) $\begin{array}{rl}u(x,y)=g\left(x\right)[1-yh(y\left)\right]& =1,\end{array}$(47) ) and (48(48) $\begin{array}{rl}v(x,y)=bxg\left(x\right)h\left(y\right)& =1,\end{array}$(48) ) are the equations for the host and parasitoid nullclines given in system (18(18a) $\begin{array}{rl}1& =u(x^{\ast },y^{\ast })=g\left(x^{\ast }\right)[1-y^{\ast }h(y^{\ast }\left)\right],\end{array}$(18a) ).

Points that satisfy Equations (47(47) $\begin{array}{rl}u(x,y)=g\left(x\right)[1-yh(y\left)\right]& =1,\end{array}$(47) ) and (48(48) $\begin{array}{rl}v(x,y)=bxg\left(x\right)h\left(y\right)& =1,\end{array}$(48) ) are coexistence equilibria and are technically denoted $(x^{\ast },y^{\ast })$. However, when we find parametric boundary curve (45(45) $R_0=\frac{y{\mathrm{e}}^{2y}}{{\mathrm{e}}^y-1},b=\frac{y^2{\mathrm{e}}^{2y}-y{\mathrm{e}}^y+y}{({\mathrm{e}}^y-1)(y{\mathrm{e}}^y-{\mathrm{e}}^y+1)},$(45) ), the variable y functions as a parameter such that the curve is traced out in $R_0$–b space for a range of y values greater than zero. The resulting curve does not require actually solving for any of the coexistence equilibrium points that correspond to the $(R_0,b)$ values on this curve. Hence, we do not use $x^{\ast }$ or $y^{\ast }$ notation here.

Caution is required with the method used to determine parametric curve (45(45) $R_0=\frac{y{\mathrm{e}}^{2y}}{{\mathrm{e}}^y-1},b=\frac{y^2{\mathrm{e}}^{2y}-y{\mathrm{e}}^y+y}{({\mathrm{e}}^y-1)(y{\mathrm{e}}^y-{\mathrm{e}}^y+1)},$(45) ). A curve that satisfies Equations (46(46) $\begin{array}{rl}\mathrm{\Delta }& =1,\end{array}$(46) )–(48(48) $\begin{array}{rl}v(x,y)=bxg\left(x\right)h\left(y\right)& =1,\end{array}$(48) ) is not necessarily a true boundary of the stability region; it could be the case that $\mathrm{\Delta }<1$ on both sides of the curve or $\mathrm{\Delta }>1$ on both sides of the curve. After finding the parametric expression for the candidate boundary curve, Equations (45(45) $R_0=\frac{y{\mathrm{e}}^{2y}}{{\mathrm{e}}^y-1},b=\frac{y^2{\mathrm{e}}^{2y}-y{\mathrm{e}}^y+y}{({\mathrm{e}}^y-1)(y{\mathrm{e}}^y-{\mathrm{e}}^y+1)},$(45) ), we verified numerically that the third Jury condition is satisfied below the curve and violated above the curve. Details for all three Jury conditions are given in Sections A.10–A.12.

As shown in the stability region in Figure 3(c), the Jury 3 condition is the interesting feature of the stability region for Model 3. Crossing the dashed curve in parameter space from below corresponds to violating the third Jury condition such that both eigenvalues leave the unit disc in the complex plane. For our model, this yields a supercritical Neimark–Sacker bifurcation [39], where the equilibrium loses stability and is replaced by a stable, quasiperiodic attractor that is topologically similar to a circle. These attractors are commonly referred to as invariant circles. The bifurcation itself is sometimes referred to as a discrete Hopf bifurcation. Crossing the line b = 1 violates the first Jury condition and results in a transcritical bifurcation. Crossing the line $R_0=1$ violates both the first and third Jury conditions.

6.2. Bifurcations and attractors

In order to illustrate the bifurcations and types of attractors for different parameter choices, we fix $R_0=2$ and increase b. Figure 6 shows the attractors for increasing b values. For b values below the Jury 3 curve in Figure 3(c), the coexistence point is stable. After the Neimark–Sacker bifurcation, the complex eigenvalues of the fixed point are larger than one in magnitude, and the attractor is either a quasiperiodic invariant circle or a periodic n-cycle, increasing in size as b increases.

Figure 6. To investigate the interior attractor for Model 3, we fix $R_0=2$ and increase b. The value of b that produces the stable equilibrium in these figures is b = 3. Moving outward from this equilibrium point, the invariant circles and phase-locked cycles correspond to b = 3.5, 4, 4.41, 5, 5.5, 6, 6.5, 7. Crossing the boundary of the stability region results in a Neimark–Sacker bifurcation such that an invariant circle becomes the stable attractor. The invariant circle grows and undergoes phase-locking alternating with invariant circles as b continues to increase. As b increases, the lower portion of the attractor approaches the x-axis.

To verify that our system exhibits a Neimark–Sacker bifurcation when the third Jury condition is violated, we numerically computed the arguments of the eigenvalues of the system for $R_0$ and b parameters on the Jury 3 curve for $R_0>1$. In doing so, we verified that the eigenvalues do not pass through the kth roots of unity for k = 1, 2, 3, 4. To guarantee that an invariant circle does appear when the third Jury condition is violated, eigenvalues cannot pass through any of the first four roots of unity [38,39].

We will now qualitatively describe the behavior of the system after the eigenvalues pass through the unit circle. For some values of b just above the Jury 3 curve, the system has a stable invariant circle. For other values of b also just above the Jury 3 curve, there is a pair of periodic orbits on the invariant circle, one stable and one unstable. When the rotation number of the periodic orbits is p/q, the system has a $p/q$ resonance, and the stable and unstable orbits have period q [3]. The numerator of the rotation number, p, is the number of revolutions around the invariant circle that it takes for the orbit to return to the starting point of the cycle. The set of parameter values for which the system has a periodic orbit with rational rotation number p/q is known as an Arnold tongue [39].

Alternately, Arnold tongues or resonance horns may describe a cusped region in the complex plane (rather than in parameter space) where eigenvalues within the horn correspond to the existence of a stable periodic orbit with rational rotation number [3,23,24]. The eigenvalue will typically intersect an infinite number of these resonance horns near the unit circle [38]. In our case, as b continues to increase, the eigenvalues of the coexistence point pass in and out of resonance horns or Arnold tongues. When the eigenvalues are within a resonance horn, the system is phase-locked, and we observe a stable n-cycle in the x–y plane.

As the eigenvalues continue to grow in magnitude, the Arnold tongues are wider, and there are broader windows of phase-locking in the bifurcation diagram as the parameter b increases. Within these windows, the system may undergo period-doubling or period-halving of the cycle as the eigenvalues enter and exit overlapping resonance horns where p and q simultaneously double or are halved. For other cases with eigenvalues in overlapping resonance horns, we may expect to observe coexistence of multiple n-cycles with differing rational rotation numbers. Note that certain n-cycles may not be visible in practice if a different n-cycle dominates the system [2,17]. Further, it is also possible in general for bistability to exist between a periodic cycle and some chaotic attractor as eigenvalues move further away from the unit circle. Aronson et al. [3] provide conditions for when a rotation number in a system is unique, as well as complications that may arise when there are multiple rotation numbers.

The behavior of system (44(44a) $\begin{array}{rl}{x}_{t+1}& =\left[\frac{R_0{x}_t}{1+(R_0-1)x_t}\right]{\mathrm{e}}^{-y_t},\end{array}$(44a) ) as eigenvalues increase in magnitude is visible in Figures 7 and 8. These figures show bifurcation diagrams of the x coordinates of the attractors to illustrate the changes in the system as b increases for fixed $R_0=2$.

Figure 7. Model 3 bifurcation diagram illustrating the x coordinates of the stable attractor for $R_0=2$ and varying b. Note the Neimark–Sacker bifurcation that occurs when the equilibrium loses stability and an invariant circle becomes the attractor, corresponding to multiple x values for a single value of b. (Bifurcation diagram for y not shown here.)

Figure 8. Model 3 bifurcation diagram with a narrower range of b values. As b increases, phase locking occurs along with period doubling and halving, interspersed with regions of invariant circles corresponding to a dense set of x coordinates of the attractor. These phenomena occur as the eigenvalues pass through the resonance horns corresponding to phase-locked n-cycles.

Returning to Figure 6, we note that as b increases, the stable attractor in the system grows, and the lower portion approaches the x-axis. While the interior of the first quadrant is invariant for system (44(44a) $\begin{array}{rl}{x}_{t+1}& =\left[\frac{R_0{x}_t}{1+(R_0-1)x_t}\right]{\mathrm{e}}^{-y_t},\end{array}$(44a) ), numerical simulations of the system for values of b much past b = 8.5 cause rounding of small positive values of y down to identically 0. Thus, numerical simulations are limited in their ability to demonstrate the behavior of the system for even finite parameter values. It is ecologically likely, however, that for sufficiently small values of y, stochastic events would wipe out the parasitoid population, and the host would go to the stable equilibrium of the host-only system. However, since the host-only equilibrium can be destabilized by the introduction of parasitoids in this parameter range, an immigrating parasitoid population could invade and return the system to the complicated dynamics until another stochastic extinction event occurs.

7. Model 4: overcompensatory host density-dependence and exponential parasitism

The fourth model uses an exponential form for both density dependence and parasitism. This corresponds to Equations (7(7) $N_{t}G\left(N_t\right)=N_t{\mathrm{e}}^{r\left(1-N_t/K\right)}.$(7) ) and (10(10) $H\left(P_t\right)=\frac{1}{P_t}\left(1-{\mathrm{e}}^{-a{P}_t}\right)$(10) ). The model is thus (49a) $\begin{array}{rl}{x}_{t+1}& =x_t{\mathrm{e}}^{r(1-x_t)}{\mathrm{e}}^{-y_t},\end{array}$(49a) (49b) $\begin{array}{rl}{y}_{t+1}& =b{x}_t{\mathrm{e}}^{r(1-x_t)}\left(1-{\mathrm{e}}^{-y_t}\right).\end{array}$(49b) Because of the exponential parasitism term, there is not an explicit expression for the coexistence equilibria solutions to system (49(49a) $\begin{array}{rl}{x}_{t+1}& =x_t{\mathrm{e}}^{r(1-x_t)}{\mathrm{e}}^{-y_t},\end{array}$(49a) ).

Unlike the previous models, b>1 is not necessary for the occurrence of a coexistence equilibrium in the interior of the first quadrant. As shown in Figure 3(d), there is a region above r = 2 and below b = 1 for which there are two coexistence equilibria, only one of which may be stable. For 0<r<2, when b = 1, the single coexistence equilibrium has collided with the exclusion equilibrium at $(1,0)$. For r = 0 ($R_0=1$), system (49(49a) $\begin{array}{rl}{x}_{t+1}& =x_t{\mathrm{e}}^{r(1-x_t)}{\mathrm{e}}^{-y_t},\end{array}$(49a) ) has a line of equilibria on the x-axis.

This model was previously studied by Kang et al. [19] in the context of a plant–herbivore system. However, our nondimensionalization and methods of analysis differ from theirs. In particular, Kang et al. [19] studied stability of the equilibria numerically, while we use analytic methods. This allows us to find a bifurcation that is missing from their analysis, discussed below.

7.1. Stability region

This model uses exponential forms for both density-dependent recruitment and parasitism. The stronger nonlinearity in density dependence and the stronger form of parasitism result in both the second and third Jury conditions functioning as interesting boundaries of the stability region, seen in Figure 3(d). The stability conditions are:

1. For 0<r<2, Jury condition 1 is satisfied above b = 1; for r>2, Jury condition 1 is satisfied above the curve (50) $r=\frac{y^2{\mathrm{e}}^y}{1+y{\mathrm{e}}^y-{\mathrm{e}}^y},b=\frac{y^2{\mathrm{e}}^y}{{\left({\mathrm{e}}^y-1\right)}^2},$(50)

2. Jury condition 2 is satisfied above the curve (51) $r=\frac{{\mathrm{e}}^y(y^2+2y+2)-2}{{\mathrm{e}}^y(y+1)-1},b=\frac{2y({\mathrm{e}}^y-1)+y^2{\mathrm{e}}^y(2+y)}{(2+y){\mathrm{e}}^{2y}-4{\mathrm{e}}^y-y+2},$(51) where $r\ge 2$,

3. Jury condition 3 is satisfied for r>0 ($R_0>1$) and below the curve (52) $r=\frac{{\mathrm{e}}^y(y^2+y-1)+1}{y{\mathrm{e}}^y},b=\frac{{\mathrm{e}}^y(y^3+y^2-y)+y}{\left({\mathrm{e}}^y-1\right)\left(y{\mathrm{e}}^y-{\mathrm{e}}^y+1\right)}.$(52)

The parametric curves (50(50) $r=\frac{y^2{\mathrm{e}}^y}{1+y{\mathrm{e}}^y-{\mathrm{e}}^y},b=\frac{y^2{\mathrm{e}}^y}{{\left({\mathrm{e}}^y-1\right)}^2},$(50) ), (51(51) $r=\frac{{\mathrm{e}}^y(y^2+2y+2)-2}{{\mathrm{e}}^y(y+1)-1},b=\frac{2y({\mathrm{e}}^y-1)+y^2{\mathrm{e}}^y(2+y)}{(2+y){\mathrm{e}}^{2y}-4{\mathrm{e}}^y-y+2},$(51) ), and (52(52) $r=\frac{{\mathrm{e}}^y(y^2+y-1)+1}{y{\mathrm{e}}^y},b=\frac{{\mathrm{e}}^y(y^3+y^2-y)+y}{\left({\mathrm{e}}^y-1\right)\left(y{\mathrm{e}}^y-{\mathrm{e}}^y+1\right)}.$(52) ) are all defined for positive y. When b = 1 and 0<r<2, the first Jury condition is violated. When r = 0 ($R_0=1$), the first and third Jury conditions are violated.

Details for determining all three conditions are given in Sections A.13–A.15. For each parametrically-defined curve, we used the equations for the host and parasitoid nullclines, Equations (18a(18a) $\begin{array}{rl}1& =u(x^{\ast },y^{\ast })=g\left(x^{\ast }\right)[1-y^{\ast }h(y^{\ast }\left)\right],\end{array}$(18a) ) and (18b(18b) $\begin{array}{rl}1& =v(x^{\ast },y^{\ast })=b{x}^{\ast }g\left(x^{\ast }\right)h\left(y^{\ast }\right).\end{array}$(18b) ), with either $1-\tau +\mathrm{\Delta }=0$, $1+\tau +\mathrm{\Delta }=0$, or $\mathrm{\Delta }=1$ to eliminate x and write b and r as functions of y. This process gave us candidate bifurcation curves. We then verified numerically that the relevant Jury condition is satisfied on one side of the curve and violated on the other side. Note also that we use y rather than $y^{\ast }$, since y functions as a parameter used to trace out a curve in r–b space.

Note that curve (50(50) $r=\frac{y^2{\mathrm{e}}^y}{1+y{\mathrm{e}}^y-{\mathrm{e}}^y},b=\frac{y^2{\mathrm{e}}^y}{{\left({\mathrm{e}}^y-1\right)}^2},$(50) ) is entirely below curve (51(51) $r=\frac{{\mathrm{e}}^y(y^2+2y+2)-2}{{\mathrm{e}}^y(y+1)-1},b=\frac{2y({\mathrm{e}}^y-1)+y^2{\mathrm{e}}^y(2+y)}{(2+y){\mathrm{e}}^{2y}-4{\mathrm{e}}^y-y+2},$(51) ) (not shown). Curves (50(50) $r=\frac{y^2{\mathrm{e}}^y}{1+y{\mathrm{e}}^y-{\mathrm{e}}^y},b=\frac{y^2{\mathrm{e}}^y}{{\left({\mathrm{e}}^y-1\right)}^2},$(50) ) and (51(51) $r=\frac{{\mathrm{e}}^y(y^2+2y+2)-2}{{\mathrm{e}}^y(y+1)-1},b=\frac{2y({\mathrm{e}}^y-1)+y^2{\mathrm{e}}^y(2+y)}{(2+y){\mathrm{e}}^{2y}-4{\mathrm{e}}^y-y+2},$(51) ) are visibly indistinguishable at the scale used in Figure 3(d). For a given r>2, the value of b must be above curve (51(51) $r=\frac{{\mathrm{e}}^y(y^2+2y+2)-2}{{\mathrm{e}}^y(y+1)-1},b=\frac{2y({\mathrm{e}}^y-1)+y^2{\mathrm{e}}^y(2+y)}{(2+y){\mathrm{e}}^{2y}-4{\mathrm{e}}^y-y+2},$(51) ) for stability. The dynamics of the system for parameters between curves (50(50) $r=\frac{y^2{\mathrm{e}}^y}{1+y{\mathrm{e}}^y-{\mathrm{e}}^y},b=\frac{y^2{\mathrm{e}}^y}{{\left({\mathrm{e}}^y-1\right)}^2},$(50) ) and (51(51) $r=\frac{{\mathrm{e}}^y(y^2+2y+2)-2}{{\mathrm{e}}^y(y+1)-1},b=\frac{2y({\mathrm{e}}^y-1)+y^2{\mathrm{e}}^y(2+y)}{(2+y){\mathrm{e}}^{2y}-4{\mathrm{e}}^y-y+2},$(51) ) are discussed in Section 7.2 below.

Returning to Figure 3(d), we consider the bifurcations that occur when the Jury conditions are violated. Crossing the dotted curve from above violates the second Jury condition, and the system undergoes a subcritical period-doubling bifurcation. Crossing the dashed curve from below corresponds to a supercritical Neimark–Sacker bifurcation [39], where the equilibrium loses stability and is replaced by a stable, invariant circle. Crossing the solid horizontal line, b = 1, 0<r<2 from above corresponds to a transcritical bifurcation where the unique coexistence equilibrium collides with the exclusion equilibrium on the x-axis.

While Kang et al. [19] use a different nondimensionalization, their parameter a is the same as our parameter b, the product of searching efficiency, parasitoid clutch size, and the host carrying-capacity. Thus, it holds from their work that for b>1, system (49(49a) $\begin{array}{rl}{x}_{t+1}& =x_t{\mathrm{e}}^{r(1-x_t)}{\mathrm{e}}^{-y_t},\end{array}$(49a) ) has a unique coexistence equilibrium. For r>2, our results indicate that above the first Jury condition curve, (50(50) $r=\frac{y^2{\mathrm{e}}^y}{1+y{\mathrm{e}}^y-{\mathrm{e}}^y},b=\frac{y^2{\mathrm{e}}^y}{{\left({\mathrm{e}}^y-1\right)}^2},$(50) ), and below b = 1, there are two coexistence equilibria. For the region above the second Jury condition curve, (51(51) $r=\frac{{\mathrm{e}}^y(y^2+2y+2)-2}{{\mathrm{e}}^y(y+1)-1},b=\frac{2y({\mathrm{e}}^y-1)+y^2{\mathrm{e}}^y(2+y)}{(2+y){\mathrm{e}}^{2y}-4{\mathrm{e}}^y-y+2},$(51) ), below both b = 1 and the third Jury condition curve, (52(52) $r=\frac{{\mathrm{e}}^y(y^2+y-1)+1}{y{\mathrm{e}}^y},b=\frac{{\mathrm{e}}^y(y^3+y^2-y)+y}{\left({\mathrm{e}}^y-1\right)\left(y{\mathrm{e}}^y-{\mathrm{e}}^y+1\right)}.$(52) ), the equilibrium point with the larger y value is stable. This region extends infinitely in the r direction since the third Jury condition curve, (52(52) $r=\frac{{\mathrm{e}}^y(y^2+y-1)+1}{y{\mathrm{e}}^y},b=\frac{{\mathrm{e}}^y(y^3+y^2-y)+y}{\left({\mathrm{e}}^y-1\right)\left(y{\mathrm{e}}^y-{\mathrm{e}}^y+1\right)}.$(52) ), remains above the second Jury condition curve, (51(51) $r=\frac{{\mathrm{e}}^y(y^2+2y+2)-2}{{\mathrm{e}}^y(y+1)-1},b=\frac{2y({\mathrm{e}}^y-1)+y^2{\mathrm{e}}^y(2+y)}{(2+y){\mathrm{e}}^{2y}-4{\mathrm{e}}^y-y+2},$(51) ), even after the curve (52(52) $r=\frac{{\mathrm{e}}^y(y^2+y-1)+1}{y{\mathrm{e}}^y},b=\frac{{\mathrm{e}}^y(y^3+y^2-y)+y}{\left({\mathrm{e}}^y-1\right)\left(y{\mathrm{e}}^y-{\mathrm{e}}^y+1\right)}.$(52) ) is below b = 1. See Figure 3(d) and details in Sections A.13–A.15. The second Jury condition curve, (51(51) $r=\frac{{\mathrm{e}}^y(y^2+2y+2)-2}{{\mathrm{e}}^y(y+1)-1},b=\frac{2y({\mathrm{e}}^y-1)+y^2{\mathrm{e}}^y(2+y)}{(2+y){\mathrm{e}}^{2y}-4{\mathrm{e}}^y-y+2},$(51) ), is the one that was missed by Kang et al. [19].

7.2. Bifurcations and attractors

In order to clearly illustrate the bifurcations and dynamics of the system, we fix r = 2.5 ($R_0\approx 12.18$) and increase b. We have chosen a value of r for which the second Jury condition curve is the lower boundary of the stability region, seen in Figure 3(d). The host and parasitoid nullclines are shown in Figure 9 with stable and unstable equilibria and other attractors for selected values of b. A bifurcation diagram for increasing b values is shown in Figure 10.

Figure 9. Phase portraits for r = 2.5 ($R_0\approx 12.18$) as b increases. The host nullcline is shown with the solid line. As b increases, the parasitoid nullcline, shown with a dotted line, changes shape. Attractors are shown with filled circles, while unstable equilibria are shown with open circles. (Note that the solid circle on the x-axis is one of two points that form a stable two-cycle on the axis.) At $b\approx 0.959$, a saddle-node bifurcation occurs. Note the tangency between the nullclines. This condition is equivalent to the condition for the first Jury condition to be satisfied. Both coexistence equilibria are initially unstable, but after the subcritical period-doubling bifurcation seen in Figure 11, the upper of the two equilibria is stable. The lower coexistence equilibrium moves towards the x-axis and collides with the exclusion equilibrium at b = 1. For b>1, the coexistence equilibrium is unique. When the third Jury condition is violated, the coexistence equilibrium undergoes a supercritical Neimark–Sacker bifurcation. A stable quasiperiodic invariant circle becomes the attractor. (a) b = 0.959. (b) b = 0.96. (c) b = 0.9615. (d) b = 0.98. (e) b = 1.05 and (f) b = 1.6.

Figure 10. Model 4 bifurcation diagram illustrating the y coordinates of the stable attractor for r = 2.5 ($R_0\approx 12.18$) and varying b. The fixed points emerge in a saddle-node bifurcation as b crosses the first Jury condition curve (50(50) $r=\frac{y^2{\mathrm{e}}^y}{1+y{\mathrm{e}}^y-{\mathrm{e}}^y},b=\frac{y^2{\mathrm{e}}^y}{{\left({\mathrm{e}}^y-1\right)}^2},$(50) ). The upper of the two equilibria undergoes a subcritical period-doubling bifurcation in which it becomes stable. The resulting unstable two cycle is shown with the dash-dot line. The dotted line is the unstable equilibrium, which crashes through exclusion equilibrium on the x-axis at b = 1. The unstable two-cycle also crashes through the x-axis. The stable equilibrium loses stability through a Neimark–Sacker bifurcation, and an invariant circle becomes the attractor, corresponding to multiple y values for a single value of b. (Bifurcation diagram for x not shown here.)

Figure 11. Model 4 bifurcation diagram with a much narrower range of b values, again for r = 2.5 ($R_0\approx 12.18$). Here, the saddle-node bifurcation is clearly visible such that both equilibria are unstable at the lowest b values for which they exist. The left figure shows the x coordinates and the right figure shows the y coordinates for the same range of b values. Dotted lines correspond to unstable equilibria. As b increases, the equilibrium with the larger y value undergoes a subcritical period-doubling bifurcation and gains stability as an unstable two cycle is born, shown with a dash-dot line. This behavior was missed in the numerical investigations by Kang et al. [19] but can be found analytically from Equations (50(50) $r=\frac{y^2{\mathrm{e}}^y}{1+y{\mathrm{e}}^y-{\mathrm{e}}^y},b=\frac{y^2{\mathrm{e}}^y}{{\left({\mathrm{e}}^y-1\right)}^2},$(50) ) and (51(51) $r=\frac{{\mathrm{e}}^y(y^2+2y+2)-2}{{\mathrm{e}}^y(y+1)-1},b=\frac{2y({\mathrm{e}}^y-1)+y^2{\mathrm{e}}^y(2+y)}{(2+y){\mathrm{e}}^{2y}-4{\mathrm{e}}^y-y+2},$(51) ). The solid line corresponds to where the equilibrium with the larger y value is stable.

For r = 2.5 and b below the first Jury condition curve, Equation (50(50) $r=\frac{y^2{\mathrm{e}}^y}{1+y{\mathrm{e}}^y-{\mathrm{e}}^y},b=\frac{y^2{\mathrm{e}}^y}{{\left({\mathrm{e}}^y-1\right)}^2},$(50) ), there are no coexistence equilibria. When we increase b to the first Jury condition curve, the host and parasitoid nullclines are tangent, seen in Figure 9(a). For slightly higher values of b, both coexistence equilibria are unstable. This differs from the claim made in Kang et al. [19] that one equilibrium is stable after the saddle-node bifurcation. However, the instability of both coexistence equilibria occurs for a tiny range of b values, from 0.959<b<0.961. The upper of the two equilibria undergoes a subcritical period-doubling bifurcation and gains stability as b crosses the dotted curve shown in Figure 3(d), the second Jury condition curve. The resulting unstable two-cycle was found numerically and is shown in Figures 10 and 11.

We continue with the bifurcations as b increases past b = 1. Returning to Figure 9, we see that at b = 1, the lower of the coexistence equilibria collides with the exclusion equilibrium as it passes into the fourth quadrant. As b continues to increase, the unstable two-cycle in the interior of the first quadrant eventually crashes through the x-axis, passing into the fourth quadrant. For sufficiently low values of b, the two-cycle on the axis is a competing stable attractor. When b crosses the third Jury condition curve, a Neimark–Sacker bifurcation results and a quasiperiodic stable invariant circle is born. As seen in Figure 10, the complex eigenvalues of the coexistence equilibrium point again pass in and out of Arnold tongues, resulting in phase-locking and stable n-cycles. A detailed discussion of this phenomena is in Section 6.2.

We note that for this model, we also see the development of a chaotic strange attractor. The collapse of the strange attractor in a crisis bifurcation is discussed by Kang et al. [19], as well as cases of more complicated bistability between boundary attractors and interior attractors. Hence, we do not discuss details here. One of the strange attractors is shown in Figure 12. Due to the use of a stronger nonlinearity in density dependence and stronger parasitism in Model 4, we see the greatest variability in dynamics and bifurcations in this system compared to Models 1, 2, and 3.

Figure 12. A strange attractor for Model 4 with r = 2.3, b = 2.2. Parameters chosen for aesthetic appeal of the strange attractor.

8. Discussion

We have developed a framework for investigating host–parasitoid systems where density dependence precedes parasitism in the life cycle of the host. Recall that these models have the form given in system (5(5a) $\begin{array}{rl}{N}_{t+1}& =N_{t}G\left(N_t\right)[1-P_{t}H(P_t\left)\right],\end{array}$(5a) ). Our analysis addresses all combinations of the most frequently used functions for host density-dependence and parasitism. The methods used in this paper can also be extended to models using other functional forms for recruitment and parasitism, including cases where there may not be an explicit expression for the coexistence equilibria. With our analytical approach, we were able to more fully categorize the dynamics of system (49(49a) $\begin{array}{rl}{x}_{t+1}& =x_t{\mathrm{e}}^{r(1-x_t)}{\mathrm{e}}^{-y_t},\end{array}$(49a) ), Model 4, which previously had been analyzed using numerical techniques [19]. Our systematic approach allows for direct comparison of four foundational models, each based on specific biological characteristics of host and parasitoid species.

Each model resulted in different dynamics. Through systematic comparison of the models, we identified the effects of stronger parasitism (corresponding to higher κ or parasitoid aggregation). We then contrasted these effects with the effects of stronger nonlinearity in the density-dependence term. As expected, fractional recruitment and parasitism yield stable dynamics. Stronger parasitism in the model leads to a restricted stability region for the coexistence equilibrium, seen in Figure 3(c,d). Both Models 3 and 4 include Neimark–Sacker bifurcations where the coexistence equilibrium is replaced with invariant circles. On the other hand, stronger nonlinearity in the density-dependence term produces period-doubling bifurcations and the potential for bistability. In the case of Model 2, the period doubling may be supercritical or subcritical, depending on the value of b. The period-doubling bifurcation observed in Model 4 is subcritical and only occurs for sufficiently large values of r.

For models with stronger parasitism resulting from higher parasitoid aggregation (Models 3 and 4), stability of the equilibrium is lost as b increases. Since b is proportional to host carrying-capacity, K, an increase in host carrying-capacity can result in loss of stability of the equilibrium for these models, consistent with the paradox of biological enrichment [34]. For the invariant circles and n-cycles that arise after the Neimark–Sacker bifurcation, the host population remains below the carrying capacity throughout the population cycles. On the other hand, the loss of stability through increased r in Model 2 yields drastic swings in host population size above and below carrying capacity, K, with relatively short period (2, 4, etc.). In these cases, the introduction of a parasitoid species could increase the host population size above its natural carrying capacity during some years of the population cycles. In agricultural scenarios, these host outbreaks could have devastating consequences.

Future work for the models presented in this paper requires comparison with data from host–parasitoid systems and consideration of what range of parameters are observed biologically. While we have provided a mathematical characterization of these systems, the biological implications need to be experimentally verified. As noted above, the period and amplitude of oscillations differ for the case of invariant circles arising in models with higher parasitism and the case of 2-cycles or 4-cycles arising in models with overcompensatory density-dependent effects. It would be beneficial to compare these models with data to determine if overcompensation does, in fact, lead to shorter-period, higher-amplitude oscillations in host population size in experimental systems.

In comparing with data, it is important to acknowledge environmental and demographic stochasticity, which will impact the ways that mathematically predicted n-cycles and quasi-periodic fluctuations in population size manifest in real populations. It is also important to consider whether the models presented here can be used for prediction in specific management scenarios or whether their use is more suited to development of biological control theory. Barlow [5] provides a survey of biological control models for specific real-world systems and emphasizes the value of models in understanding specific case studies, whether or not the models are used for practical management decisions.

As discussed in Section 1, the sequence of events in the host life-cycle also has important impacts on the population dynamics. The models investigated in this paper assume that density dependence precedes parasitism, which is an appropriate assumption for some species. For example, houseflies (Musca spp.) are attacked by pupal parasitoids such as Spalangia spp. and Muscidifurax spp. after significant density dependence in the early larval stages [28]. However, in other species, density dependence acts on the survivors of parasitism, which leads to the model (53a) $\begin{array}{rl}{N}_{t+1}& =N_t[1-P_{t}H(P_t\left)\right]G(N_t[1-P_{t}H(P_t\left)\right]),\end{array}$(53a) (53b) $\begin{array}{rl}{P}_{t+1}& =c{N}_t{P}_{t}H\left(P_t\right).\end{array}$(53b) Note that this model assumes not only that parasitism occurs first in the host life-cycle, but also that the parasitized hosts are functionally dead and unable to compete. For the fractional form for recruitment and the negative binomial form for parasitism, this model can lead to stable equilibria with hosts at a higher level than their carrying capacity in the absence of parasitoids [28,29].

Systems (53(53a) $\begin{array}{rl}{N}_{t+1}& =N_t[1-P_{t}H(P_t\left)\right]G(N_t[1-P_{t}H(P_t\left)\right]),\end{array}$(53a) ) and (5(5a) $\begin{array}{rl}{N}_{t+1}& =N_{t}G\left(N_t\right)[1-P_{t}H(P_t\left)\right],\end{array}$(5a) ) represent the scenarios where parasitism occurs either before or after density-dependent effects on the host, and May et al. [28] compared some specific models in these frameworks. However, more complicated parasitoid phenologies exist in nature. Cobbold et al. [7] explicitly consider koinobiont parasitoids, which do not kill their host immediately. This means that there is a period of time when parasitized hosts are competing with nonparasitized hosts, which cannot be accounted for with either system (53(53a) $\begin{array}{rl}{N}_{t+1}& =N_t[1-P_{t}H(P_t\left)\right]G(N_t[1-P_{t}H(P_t\left)\right]),\end{array}$(53a) ) or system (5(5a) $\begin{array}{rl}{N}_{t+1}& =N_{t}G\left(N_t\right)[1-P_{t}H(P_t\left)\right],\end{array}$(5a) ). Cobbold et al. [7] found that the delayed mortality of parasitized hosts may have implications for biological control. Differences in the timing of interaction between parasitoids and hosts lead to different predicted population dynamics. Thus, model formulation requires care and awareness of biological assumptions that are inherent to the structure of a model.

Host–parasitoid models have numerous avenues for the inclusion of additional biological complexities such as spatial heterogeneity, Allee effects, and multiple parasitoid species. In building towards these more biologically realistic models, it is important to understand the dynamics of simpler models, such as those analyzed and compared here. Hassell [13,14] has done excellent work in bridging the gap between simple mechanistic models for host–parasitoid systems and models for more complex and biologically realistic systems. Extending simple mechanistic models to investigate more complicated scenarios can only occur when the simple foundational models are well-understood and presented with explicit acknowledgment of biological assumptions.

Biological differences between models may be critical to communicate well with ecologists and experimentalists. We therefore urge researchers to exercise caution in formulation of models and underlying biological assumptions in order to promote communication and broader understanding of mathematical and theoretical findings.

Acknowledgments

The authors would like to express gratitude to the editor and reviewers for helpful comments and suggestions. Additionally, we thank Benjamin Liu for troubleshooting and assistance with formatting figures.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Appendices

Appendix A. Partial derivatives for Jury conditions

We begin by evaluating the partial derivatives that appear in the Jury conditions. From the definitions of $u(x,y)$ and $v(x,y)$, Equations (12(12) $\begin{array}{rl}u(x_t,y_t)& =g\left(x_t\right)[1-y_{t}h(y_t\left)\right],\end{array}$(12) ) and (13(13) $\begin{array}{rl}v(x_t,y_t)& =b{x}_{t}g\left(x_t\right)h\left(y_t\right).\end{array}$(13) ), we obtain the partial derivatives (A1) $\begin{array}{rl}{u}_x& =g^{\prime }\left(x\right)[1-yh(y\left)\right],\end{array}$(A1) (A2) $\begin{array}{rl}{u}_y& =g\left(x\right)[1-yh(y){]}^{\prime },\end{array}$(A2) (A3) $\begin{array}{rl}{v}_x& =b\left[g\right(x)+x{g}^{\prime }(x\left)\right]h\left(y\right),\end{array}$(A3) (A4) $\begin{array}{rl}{v}_y& =bxg\left(x\right)h^{\prime }\left(y\right).\end{array}$(A4) For Models 1 and 2, with fractional parasitism (16(16) $h\left(y_t\right)=\frac{1}{1+y_t},$(16) ), we use (A5) $h\left(y\right)=\frac{1}{1+y}=1-yh\left(y\right).$(A5) Recall that Models 3 and 4 have exponential parasitism, given by Equation (17(17) $h\left(y_t\right)=\frac{1}{y_t}\left(1-{\mathrm{e}}^{-y_t}\right),$(17) ). Furthermore, Models 1 and 3 use fractional per-capita-recruitment, with $g\left(x\right)$ defined in Equation (14(14) $g\left(x_t\right)=\frac{R_0}{1+(R_0-1)x_t},$(14) ), while Models 2 and 4 use exponential per-capita-recruitment, with $g\left(x\right)$ defined in Equation (15(15) $g\left(x_t\right)={\mathrm{e}}^{r(1-x_t)}.$(15) ). To simplify expressions for the partial derivatives for each model using the corresponding functions for $g\left(x\right)$ and $h\left(y\right)$, we also use the nullcline equations, $u(x^{\ast },y^{\ast })=1$ and $v(x^{\ast },y^{\ast })=1$. The resulting expressions for the partial derivatives are given in Table 1.

Table A1. Partial derivatives used to apply the Jury conditions to the coexistence equilibrium point(s) for each model.

	Model 1	Model 2	Model 3	Model 4
$u_x$	${\displaystyle \frac{(1-R_0)}{R_0}}g\left(x^{\ast }\right)$	$-r$	${\displaystyle \frac{(1-R_0)}{R_0}}g\left(x^{\ast }\right)$	$-r$
$u_y$	$-h\left(y^{\ast }\right)$	$-h\left(y^{\ast }\right)$	$-1$	$-1$
$v_x$	${\displaystyle \frac{1}{R_0{x}^{\ast }}}g\left(x^{\ast }\right)$	${\displaystyle \frac{1}{x^{\ast }}}-r$	${\displaystyle \frac{1}{R_0{x}^{\ast }}}g\left(x^{\ast }\right)$	${\displaystyle \frac{1}{x^{\ast }}}-r$
$v_y$	$-h\left(y^{\ast }\right)$	$-h\left(y^{\ast }\right)$	${\displaystyle \frac{1}{y^{\ast }h\left(y^{\ast }\right)}}[1-y^{\ast }h(y^{\ast })-h(y^{\ast }\left)\right]$	${\displaystyle \frac{1}{y^{\ast }h\left(y^{\ast }\right)}}[1-y^{\ast }h(y^{\ast })-h(y^{\ast }\left)\right]$

Appendix B. Model 1 stability calculations

B.1. Requirements for existence of coexistence equilibrium in the first quadrant

We now determine the conditions that ensure an equilibrium in the interior of the first quadrant. For this model, we can explicitly solve system (18(18a) $\begin{array}{rl}1& =u(x^{\ast },y^{\ast })=g\left(x^{\ast }\right)[1-y^{\ast }h(y^{\ast }\left)\right],\end{array}$(18a) ) for the coexistence equilibrium, (A6) $(x^{\ast },y^{\ast })=\left(\frac{1}{b},g\left(\frac{1}{b}\right)-1\right).$(A6) The x coordinate is positive for all positive b. The y coordinate is positive when (A7) $g\left(\frac{1}{b}\right)-1=\frac{R_0}{1+\left(R_0-1\right)\left(\frac{1}{b}\right)}-1>0,$(A7) which simplifies to b>1 since we assume $R_0>1$. Thus, the coexistence equilibrium exists and is in the first quadrant when $R_0>1,b>1$.

B.2. First Jury condition

For b>1, $R_0>1$, the x and y coordinates of the coexistence equilibrium are positive. We thus use partial derivatives from Table 1 to write inequality (36(36) $u_x{v}_y-u_y{v}_x>0.$(36) ), as (A8) $u_x{v}_y-u_y{v}_x=\left(\frac{R_0-1}{R_0}\right)g\left(x^{\ast }\right)h\left(y^{\ast }\right)+h\left(y^{\ast }\right)\left[\frac{1}{R_0{x}^{\ast }}g\left(x^{\ast }\right)\right]>0.$(A8) We now substitute in the expression for $g\left(x^{\ast }\right)$ from Equation (14(14) $g\left(x_t\right)=\frac{R_0}{1+(R_0-1)x_t},$(14) ) and simplify on the left-hand side to get (A9) $\frac{h\left(y^{\ast }\right)}{x^{\ast }}>0.$(A9) Since $h\left(y^{\ast }\right)$ is positive, the first Jury condition is satisfied whenever the coexistence equilibrium is in the first quadrant.

When b = 1, the y-coefficient from Equation (A6(A6) $(x^{\ast },y^{\ast })=\left(\frac{1}{b},g\left(\frac{1}{b}\right)-1\right).$(A6) ) is $y^{\ast }=0$, and the first Jury condition, inequality (35(35) $x^{\ast }{y}^{\ast }(u_x{v}_y-u_y{v}_x)>0.$(35) ) is violated. When $R_0=1$, Equation (A6(A6) $(x^{\ast },y^{\ast })=\left(\frac{1}{b},g\left(\frac{1}{b}\right)-1\right).$(A6) ) again gives us $y^{\ast }=0$, regardless of the value of b, such that the x-axis is a line of equilibrium points. For $R_0=1$, the first Jury condition, inequality (35(35) $x^{\ast }{y}^{\ast }(u_x{v}_y-u_y{v}_x)>0.$(35) ) is again violated.

B.3. Second Jury condition

Recall that the second Jury condition, inequality (31(31) $\begin{array}{rl}& 1+\tau +\mathrm{\Delta }>0,\end{array}$(31) ), is (A10) $1+\tau +\mathrm{\Delta }>0.$(A10) For this model, we will not show this directly. Instead, note that if $\tau >0$ and $1-\tau +\mathrm{\Delta }>0$, which is the first Jury condition, then $1+\tau +\mathrm{\Delta }>1-\tau +\mathrm{\Delta }>0$. This means $\tau >0$ and the satisfaction of the first Jury condition are sufficient criteria for the second Jury condition.

The first Jury condition is satisfied for $b>1,R_0>1$. We will show that in this case, the second Jury condition will also be satisfied. We proceed by showing that $\tau >0$ at the equilibrium. As seen in Equation (33(33) $\tau =2+x^{\ast }{u}_x+y^{\ast }{v}_y,$(33) ), $\tau =2+x^{\ast }{u}_x+y^{\ast }{v}_y$. We use the expressions for $u_x$ and $v_y$ from Table 1 and the definitions of $g\left(x\right)$ and $h\left(y\right)$ from Equations (14(14) $g\left(x_t\right)=\frac{R_0}{1+(R_0-1)x_t},$(14) ) and (16(16) $h\left(y_t\right)=\frac{1}{1+y_t},$(16) ) to express the trace, (A11) $\tau =2-\frac{R_0-1}{R_0}{x}^{\ast }g\left(x^{\ast }\right)-y^{\ast }h\left(y^{\ast }\right)=2-\left[\frac{(R_0-1)x^{\ast }}{1+(R_0-1)x^{\ast }}+\frac{y^{\ast }}{1+y^{\ast }}\right].$(A11) We thus seek to show that (A12) $2>\frac{(R_0-1)x^{\ast }}{1+(R_0-1)x^{\ast }}+\frac{y^{\ast }}{1+y^{\ast }}$(A12) for $x^{\ast },y^{\ast }>0$, $R_0>1$.

Both of the terms on the right-hand side of inequality (A12(A12) $2>\frac{(R_0-1)x^{\ast }}{1+(R_0-1)x^{\ast }}+\frac{y^{\ast }}{1+y^{\ast }}$(A12) ) are of the form $z(1+z{)}^{-1}$, where z is positive. Each term individually is less than one because z<1 + z, which indicates that $z(1+z{)}^{-1}<1$ for positive z. Therefore, (A13) $\tau =2-\left[\frac{(R_0-1)x^{\ast }}{1+(R_0-1)x^{\ast }}+\frac{y^{\ast }}{1+y^{\ast }}\right]>0.$(A13) It follows that the first Jury condition is a sufficient condition for the second Jury condition for Model 1.

For either b = 1 or $R_0=1$, we can directly calculate the terms in the second Jury condition, inequality (31(31) $\begin{array}{rl}& 1+\tau +\mathrm{\Delta }>0,\end{array}$(31) ). Direct calculation verifies that the second Jury condition is satisfied.

B.4. Third Jury condition

Recall that the third Jury condition is $\mathrm{\Delta }<1$. The determinant is given in terms of the partial derivatives in Equation (34(34) $\mathrm{\Delta }=1+x^{\ast }{u}_x+y^{\ast }{v}_y+x^{\ast }{y}^{\ast }(u_x{v}_y-u_y{v}_x).$(34) ). Using the expressions for $u_x,u_y,v_x$, and $v_y$ from Table 1, the third Jury condition simplifies to (A14) $1+\frac{(1-R_0)}{R_0}{x}^{\ast }g\left(x^{\ast }\right)<1.$(A14) When we use Equation (14(14) $g\left(x_t\right)=\frac{R_0}{1+(R_0-1)x_t},$(14) ) for $g\left(x^{\ast }\right)$, the condition can be expressed as (A15) $1-\frac{(R_0-1)x^{\ast }}{1+(R_0-1)x^{\ast }}<1,$(A15) which simplifies to (A16) $1+(R_0-1)x^{\ast }>1.$(A16) This is true for $R_0>1$ for the equilibrium in the interior of the first quadrant. Thus, the third Jury condition is satisfied for $R_0>1$. For $R_0=1$, the third Jury condition is violated.

Appendix C. Model 2 stability calculations

C.1. Requirements for existence of coexistence equilibrium in the first quadrant

For this model, we can again explicitly solve system (18(18a) $\begin{array}{rl}1& =u(x^{\ast },y^{\ast })=g\left(x^{\ast }\right)[1-y^{\ast }h(y^{\ast }\left)\right],\end{array}$(18a) ) for the coexistence equilibrium for Model 2, (A17) $(x^{\ast },y^{\ast })=\left(\frac{1}{b},{\mathrm{e}}^{r\left(1-1/b\right)}-1\right).$(A17) When b = 1, this equilibrium point is on the x-axis at $(x^{\ast },y^{\ast })=(1,0)$, which is the exclusion equilibrium. For the coexistence equilibrium to be in the interior of the first quadrant, it is necessary that (A18) $y^{\ast }={\mathrm{e}}^{r\left(1-1/b\right)}-1>0.$(A18) For r>0, this requires b>1. Note that we will not consider the case r<0, b<1 since we are interested in cases where the host species persists in the absence of the parasitoid.

C.2. First Jury condition: slopes of zero-growth isoclines

Using partial derivatives from Table 1, the first Jury condition, inequality (35(35) $x^{\ast }{y}^{\ast }(u_x{v}_y-u_y{v}_x)>0.$(35) ), is (A19) $x^{\ast }{y}^{\ast }\left[rh\left(y^{\ast }\right)+h\left(y^{\ast }\right)\left(\frac{1}{x^{\ast }}-r\right)\right]=y^{\ast }h\left(y^{\ast }\right)>0.$(A19) Because $y^{\ast }h\left(y^{\ast }\right)$ is positive for the coexistence equilibrium, this inequality holds for the equilibrium in the interior of the first quadrant. When b = 1, $y^{\ast }h\left(y^{\ast }\right)=0$, and the first Jury condition is violated. For r = 0, the x-axis is a line of equilibrium points, and the first Jury condition is again violated.

C.3. Second Jury condition

Again using partial derivatives from Table 1, the second Jury condition, inequality (37(37) $4+2x^{\ast }{u}_x+2y^{\ast }{v}_y+x^{\ast }{y}^{\ast }(u_x{v}_y-u_y{v}_x)>0.$(37) ) simplifies to (A20) $4-2x^{\ast }r-y^{\ast }h\left(y^{\ast }\right)=4-2xr-\frac{y^{\ast }}{1+y^{\ast }}>0.$(A20) The coordinates of the coexistence equilibrium point are given by Equation (A17(A17) $(x^{\ast },y^{\ast })=\left(\frac{1}{b},{\mathrm{e}}^{r\left(1-1/b\right)}-1\right).$(A17) ). Using these values, the stability condition is (A21) $3-\frac{2}{b}r+{\mathrm{e}}^{\left(r/b-r\right)}>0.$(A21) We now consider the transcendental equation, (A22) $3-\frac{2}{b}r+{\mathrm{e}}^{\left(r/b-r\right)}=0,$(A22) and introduce the parameter $\omega =r/b$ so that (A23) $3-2\omega +{\mathrm{e}}^{\omega -r}=0.$(A23) We solve for r as a function of ω, (A24) $r=\omega -\ln \left(2\omega -3\right),$(A24) and can then also write b as a function of ω, (A25) $b=\frac{r}{\omega }=1-\frac{1}{\omega }\ln \left(2\omega -3\right).$(A25) For $\omega >3/2$, Equations (A24(A24) $r=\omega -\ln \left(2\omega -3\right),$(A24) ) and (A25(A25) $b=\frac{r}{\omega }=1-\frac{1}{\omega }\ln \left(2\omega -3\right).$(A25) ) express the boundary of the region in parameter space where the coexistence equilibrium satisfies the second Jury condition. The point $(r,b)=(2,1)$ is where the Jury 2 curve intersects the b = 1 line, seen in Figure 3(b). The portion of the curve shown in Figure 3(b) corresponds to $3/2<\omega <2$. We numerically verified that the second Jury condition is satisfied above this curve and violated directly below it. (The portion of the curve defined by $\omega >2$ corresponds to b<1, such that there is no coexistence equilibrium in the interior of the first quadrant.) For b = 1, inequality (A21(A21) $3-\frac{2}{b}r+{\mathrm{e}}^{\left(r/b-r\right)}>0.$(A21) ) requires r<2.

C.4. Third Jury condition

The expression for the determinant from Equation (34(34) $\mathrm{\Delta }=1+x^{\ast }{u}_x+y^{\ast }{v}_y+x^{\ast }{y}^{\ast }(u_x{v}_y-u_y{v}_x).$(34) ) for this model simplies significantly to (A26) $\mathrm{\Delta }=1-r{x}^{\ast },$(A26) using the partial derivatives in Table 1. Since $x^{\ast }=1/b$, the third Jury condition is (A27) $1-\frac{r}{b}<1.$(A27) Since b>0 and we assumed $r\ge 0$, this inequality is satisfied for r>0. When r = 0, the third Jury condition is violated.

Appendix D. Model 3 stability calculations

D.1. Requirements for existence of coexistence equilibrium in first quadrant

As was true in Section A.1, we seek to determine the conditions that ensure that an equilibrium exists in the interior of the first quadrant, this time for Model 3, system (44(44a) $\begin{array}{rl}{x}_{t+1}& =\left[\frac{R_0{x}_t}{1+(R_0-1)x_t}\right]{\mathrm{e}}^{-y_t},\end{array}$(44a) ). The coexistence equilibrium cannot be solved for explicitly in this case, so we instead consider the nullclines.

Equation (18a(18a) $\begin{array}{rl}1& =u(x^{\ast },y^{\ast })=g\left(x^{\ast }\right)[1-y^{\ast }h(y^{\ast }\left)\right],\end{array}$(18a) ) is the host nullcline with intercepts $(0,\ln R_0)$ and $(1,0)$. To obtain the slope of this nullcline in the x–y plane, we first differentiate $u(x,y)=1$ with respect to x and get (A28) $u_x+u_y\frac{\mathrm{d}y}{\mathrm{d}x}=0.$(A28) The slope of the host nullcline is (A29) $\frac{\mathrm{d}y}{\mathrm{d}x}=-\frac{u_x}{u_y}=\frac{1-R_0}{R_0}g\left(x\right),$(A29) using the expressions for $u_x$ and $u_y$ from Table 1. Since $g\left(x\right)>0$ and we assume $R_0>1$, the host nullcline is monotone decreasing in the first quadrant from $(0,\ln R_0)$ to $(1,0)$.

We now consider the parasitoid nullcline, Equation (18b(18b) $\begin{array}{rl}1& =v(x^{\ast },y^{\ast })=b{x}^{\ast }g\left(x^{\ast }\right)h\left(y^{\ast }\right).\end{array}$(18b) ). To find the slope in the x-y plane, we differentiate $v(x,y)=1$ with respect to x to get (A30) $v_x+v_y\frac{\mathrm{d}y}{\mathrm{d}x}=0.$(A30) The slope for the parasitoid nullcline is thus (A31) $\frac{\mathrm{d}y}{\mathrm{d}x}=\frac{-v_x}{v_y}=-\frac{g\left(x\right)}{R_{0}x{v}_y}.$(A31) Since $g\left(x\right)>0$, the sign of $v_y$ will determine the sign of the slope of the parasitoid nullcline. Negative $v_y$ will indicate that the slope of the nullcline is positive.

We substitute $h\left(y\right)$ from Equation (17(17) $h\left(y_t\right)=\frac{1}{y_t}\left(1-{\mathrm{e}}^{-y_t}\right),$(17) ) into $v_y$ for Model 3, such that (A32) $v_y=\frac{1}{1-{\mathrm{e}}^{-y}}\left[{\mathrm{e}}^{-y}-\frac{1}{y}\left(1-{\mathrm{e}}^{-y}\right)\right]=\frac{1}{y\left(1-{\mathrm{e}}^{-y}\right)}\left(y{\mathrm{e}}^{-y}-1+{\mathrm{e}}^{-y}\right).$(A32) The denominator is positive for y>0, so we consider the numerator. For y>0, (A33a) $\begin{array}{rl}1+y& <{\mathrm{e}}^y,\end{array}$(A33a) (A33b) $\begin{array}{rl}(1+y){\mathrm{e}}^{-y}& <1,\end{array}$(A33b) (A33c) $\begin{array}{rl}{\mathrm{e}}^{-y}+y{\mathrm{e}}^{-y}-1& <0.\end{array}$(A33c) Thus, we conclude that $v_y<0$ for y>0. This means that the slope of the parasitoid nullcline is positive in the first quadrant. If there is an intersection of the host and parasitoid nullclines in the first-quadrant, it is unique.

To determine existence of the equilibrium, we examine the the x- and y-intercepts of the parasitoid nullcline, (A34) $1=\frac{bx{R}_{0}h\left(y\right)}{1+(R_0-1)x}.$(A34) After solving for x, we obtain (A35) $x=\frac{1}{b{R}_{0}h\left(y\right)+(1-R_0)},$(A35) where (A36) $h\left(y\right)=\frac{1}{y}\left(1-{\mathrm{e}}^{-y}\right).$(A36) To examine Equation (A35(A35) $x=\frac{1}{b{R}_{0}h\left(y\right)+(1-R_0)},$(A35) ), we consider the limiting behavior of $h\left(y\right)$ as $y\to -\mathrm{\infty }$, (A37) $\underset{y\to -\mathrm{\infty }}{lim}h\left(y\right)=\underset{z\to \mathrm{\infty }}{lim}h(-z)=\underset{z\to \mathrm{\infty }}{lim}\frac{-1}{z}\left(1-{\mathrm{e}}^z\right)=\mathrm{\infty }.$(A37) Thus, if we consider the limit as $y\to -\mathrm{\infty }$ in Equation (A35(A35) $x=\frac{1}{b{R}_{0}h\left(y\right)+(1-R_0)},$(A35) ), $x\to 0^{+}$. This nullcline does not have a y-intercept because as $x\to 0^{+}$, $y\to -\mathrm{\infty }$.

Since we know that in the first quadrant, the parasitoid nullcline is monotone increasing and the host nullcline has x-intercept at x = 1, we need to find the conditions for which the parasitoid nullcline's x-intercept lies between 0 and 1. For these conditions, there exists exactly one intersection of the parasitoid and host nullclines in the interior of the first quadrant. The x-intercept of the parasitoid nullcline is the solution to (A38) $1=bxg\left(x\right)h\left(0\right)=\frac{bx{R}_0}{1+(R_0-1)x},$(A38) which is (A39) $x_{\mathrm{i}\mathrm{n}\mathrm{t}}=\frac{1}{R_0(b-1)+1}.$(A39) We seek the conditions for which (A40) $0<\frac{1}{R_0(b-1)+1}<1.$(A40) This translates into the following two criteria, (A41) $R_0(b-1)+1>0,$(A41) and (A42) $R_0(b-1)+1>1,$(A42) which can be consolidated as (A43) $\begin{array}{rl}{R}_0(b-1)+1& >1,\\ R_0(b-1)& >0.\end{array}$(A43) This is true when b>1. So the x-intercept of the parasitoid nullcline occurs between 0 and 1 if and only if b>1.

We conclude that there is exactly one equilibrium point in the interior of the first quadrant if and only if b>1. When b>1, we can then determine if the coexistence equilibrium is stable. For b = 1, the equilibrium point is on the boundary of the first quadrant, at $(1,0)$. For $R_0>1$, b<1, there are no equilibria points in the interior of the first quadrant.

D.2. First Jury condition

For b>1, $R_0>1$, the x and y coordinates of the coexistence equilibrium are positive. We thus use partial derivatives from Table 1 to write the first Jury condition, inequality (36(36) $u_x{v}_y-u_y{v}_x>0.$(36) ), as (A44) $x^{\ast }{y}^{\ast }(u_x{v}_y-u_y{v}_x)=x^{\ast }{y}^{\ast }\left[\frac{1-R_0}{R_0}g\left(x^{\ast }\right)v_y+\frac{1}{R_0{x}^{\ast }}g\left(x^{\ast }\right)\right]>0,$(A44) which simplifies to (A45) $\frac{1}{R_0{x}^{\ast }}g\left(x^{\ast }\right)>\frac{R_0-1}{R_0}g\left(x^{\ast }\right)v_y.$(A45) The left-hand side of inequality (A45(A45) $\frac{1}{R_0{x}^{\ast }}g\left(x^{\ast }\right)>\frac{R_0-1}{R_0}g\left(x^{\ast }\right)v_y.$(A45) ) is positive, while the right-hand side is negative for $R_0>1$, since $v_y$ was shown to be negative in Section A.9. Thus, this inequality holds for the positive coexistence equilibrium.

When b = 1, the coexistence equilibrium has collided with the exclusion equilibrium at $(1,0)$. Since the y coordinate is 0, the first Jury condition, inequality (35(35) $x^{\ast }{y}^{\ast }(u_x{v}_y-u_y{v}_x)>0.$(35) ), is violated. When $R_0=1$, any point on the x-axis is a solution to system (44(44a) $\begin{array}{rl}{x}_{t+1}& =\left[\frac{R_0{x}_t}{1+(R_0-1)x_t}\right]{\mathrm{e}}^{-y_t},\end{array}$(44a) ). Since these equilibria points have y = 0, the first Jury condition is violated for this line of equilibrium points.

D.3. Second Jury condition

We will use the technique from Section A.3 to show that the first Jury condition is a sufficient condition for the second Jury condition. To do this, we must show that $\tau >0$ at the interior equilibrium. As seen in Equation (33(33) $\tau =2+x^{\ast }{u}_x+y^{\ast }{v}_y,$(33) ), $\tau =2+x^{\ast }{u}_x+y^{\ast }{v}_y$. We use the expressions for $u_x$ and $v_y$ from Table 1 and the definitions of $g\left(x\right)$ and $h\left(y\right)$ from Equations (14(14) $g\left(x_t\right)=\frac{R_0}{1+(R_0-1)x_t},$(14) ) and (17(17) $h\left(y_t\right)=\frac{1}{y_t}\left(1-{\mathrm{e}}^{-y_t}\right),$(17) ) to get (A46) $\tau =1-\left[\frac{(R_0-1)x^{\ast }}{1+(R_0-1)x^{\ast }}\right]+\frac{y^{\ast }{\mathrm{e}}^{-y^{\ast }}}{1-{\mathrm{e}}^{-y^{\ast }}},$(A46) after simplification.

For positive $y^{\ast }$, the last term is positive. Similarly to Section A.3, the middle term is of the form $z(1+z{)}^{-1}$, where z is positive. For $R_0>1$, this term individually is less than one because z<1 + z, which indicates that $z(1+z{)}^{-1}<1$ for positive z. Since the coordinates of the coexistence equilibrium are positive for b>1, $R_0>1$, we conclude that $\tau >0$ for b>1 and $R_0>1$. It follows that the first Jury condition is a sufficient condition for the second Jury condition.

For either b = 1 or $R_0=1$, we can directly calculate the terms in the second Jury condition, inequality (31(31) $\begin{array}{rl}& 1+\tau +\mathrm{\Delta }>0,\end{array}$(31) ). Direct calculation verifies that the second Jury condition is satisfied in these cases.

D.4. Third Jury condition

The third Jury condition is $\mathrm{\Delta }<1$. The determinant is given in terms of the partial derivatives in Equation (34(34) $\mathrm{\Delta }=1+x^{\ast }{u}_x+y^{\ast }{v}_y+x^{\ast }{y}^{\ast }(u_x{v}_y-u_y{v}_x).$(34) ). Using the expressions for $u_x,u_y,v_x$, and $v_y$ from Table 1 and much algebraic simplification, the third Jury condition is (A47) $\mathrm{\Delta }=\frac{g\left(x^{\ast }\right)}{R_{0}h\left(y^{\ast }\right)}<1.$(A47) We want to write the condition solely in terms of the parameters, $R_0$ and b, and find the curve in the $R_0$–b plane where stability changes. Because of the transcendental nature of the inequality, we will express this curve parametrically with $R_0$ and b as functions of y. (As explained in Section 6.1, we do not use $x^{\ast }$ and $y^{\ast }$ here since y will function as a parameter, rather than the y-coordinate value of a particular coexistence equilibrium.)

We first consider (A48) $1=\frac{g\left(x\right)}{R_{0}h\left(y\right)}=\frac{1}{h\left(y\right)\left[1+(R_0-1)x\right]},$(A48) and solve for x, (A49) $x=\left(\frac{1}{R_0-1}\right)\left[\frac{1-h\left(y\right)}{h\left(y\right)}\right].$(A49) We now incorporate the host nullcline, Equation (18a(18a) $\begin{array}{rl}1& =u(x^{\ast },y^{\ast })=g\left(x^{\ast }\right)[1-y^{\ast }h(y^{\ast }\left)\right],\end{array}$(18a) ), which is valid at the equilibrium point. Using $g\left(x\right)=R_{0}h\left(y\right)$ from Equation (A48(A48) $1=\frac{g\left(x\right)}{R_{0}h\left(y\right)}=\frac{1}{h\left(y\right)\left[1+(R_0-1)x\right]},$(A48) ), we get (A50) $1=g\left(x\right)[1-yh(y\left)\right]=R_{0}h\left(y\right)[1-yh(y\left)\right].$(A50) Solving for $R_0$ as a function of y yields (A51) $R_0=\frac{1}{h\left(y\right)[1-yh(y\left)\right]}.$(A51) Next, we need an expression for b as a function of y. To do this, we incorporate the parasitoid nullcline, Equation (18b(18b) $\begin{array}{rl}1& =v(x^{\ast },y^{\ast })=b{x}^{\ast }g\left(x^{\ast }\right)h\left(y^{\ast }\right).\end{array}$(18b) ), which is valid at the equilibrium point. Starting with Equation (18b(18b) $\begin{array}{rl}1& =v(x^{\ast },y^{\ast })=b{x}^{\ast }g\left(x^{\ast }\right)h\left(y^{\ast }\right).\end{array}$(18b) ), we replace x with the expression from (A49(A49) $x=\left(\frac{1}{R_0-1}\right)\left[\frac{1-h\left(y\right)}{h\left(y\right)}\right].$(A49) ) and also substitute $R_{0}h\left(y\right)$ for $g\left(x\right)$, using the determinant condition (A48(A48) $1=\frac{g\left(x\right)}{R_{0}h\left(y\right)}=\frac{1}{h\left(y\right)\left[1+(R_0-1)x\right]},$(A48) ). This gives us the equation, (A52) $1=bxg\left(x\right)h\left(y\right)=\frac{b}{R_0-1}\left[\frac{1-h\left(y\right)}{h\left(y\right)}\right]R_{0}h\left(y\right)h\left(y\right).$(A52) We solve for b to get (A53) $b=\left(\frac{R_0-1}{R_0}\right)\frac{1}{h\left(y\right)\left[1-h\left(y\right)\right]}.$(A53) We then eliminate the dependence on $R_0$ from the equation for b. This gives us b as a function of y, (A54) $b=\frac{1-h\left(y\right)[1-yh(y\left)\right]}{h\left(y\right)[1-h(y\left)\right]},$(A54) which does not simplify in a meaningful way. This equation combined with Equation (A51(A51) $R_0=\frac{1}{h\left(y\right)[1-yh(y\left)\right]}.$(A51) ) expresses the boundary of the region in parameter space where the coexistence equilibrium satisfies the third Jury condition. We numerically verified that the third Jury condition is satisfied above this curve and violated below the curve.

When $R_0=1$, the x-axis is a line of equilibrium points, as stated in Section A.10. Under these conditions, the expression for the determinant simplifies to $\mathrm{\Delta }=1$, and so the third Jury condition is also violated for $R_0=1$.

Appendix E. Model 4 stability calculations

From Section 7, recall that there may be one or two coexistence equilibria, depending on the parameter values. In the case of two coexistence equilibria, only the point with the larger y value may be stable, as discussed in Section 7.1. The analysis here pertains to the stability of the single unique coexistence equilibrium or the coexistence equilibrium point with the larger y value.

As Kang et al. [19] proved, for b>1, system (49(49a) $\begin{array}{rl}{x}_{t+1}& =x_t{\mathrm{e}}^{r(1-x_t)}{\mathrm{e}}^{-y_t},\end{array}$(49a) ) has a unique positive equilibrium. For b = 1, 0<r<2, the point $(1,0)$ is an equilibrium point, and there is no coexistence equilibrium in the interior of the first quadrant. For r>2 and b just less than 1, the system has both a stable coexistence equilibrium point and an unstable coexistence equilibrium point, as seen in Figure 9(d). For r>2 and b = 1, the unstable coexistence point collides with the exclusion equilibrium, $(1,0)$. This is all consistent with the analysis in Kang et al. [19]. However, for r>2, there is also a narrow range of b for which the system has two unstable coexistence equilibria, as in Figure 9(b).

E.1. First Jury condition

For b>1, r>0, the x and y coordinates of the coexistence equilibrium are positive. Using partial derivatives from Table 1, the first Jury condition, inequality (35(35) $x^{\ast }{y}^{\ast }(u_x{v}_y-u_y{v}_x)>0.$(35) ), is (A55) $x^{\ast }{y}^{\ast }\left\{-r\left[\frac{1}{y^{\ast }h\left(y^{\ast }\right)}\right]\left[1-y^{\ast }h\left(y^{\ast }\right)-h\left(y^{\ast }\right)\right]+\frac{1}{x^{\ast }}-r\right\}>0,$(A55) which simplifies to (A56) $r\left(\frac{{\mathrm{e}}^{y^{\ast }}-1-y^{\ast }{\mathrm{e}}^{y^{\ast }}}{y^{\ast }{\mathrm{e}}^{y^{\ast }}-y^{\ast }}\right)+\frac{1}{x^{\ast }}>0,$(A56) using Equation (17(17) $h\left(y_t\right)=\frac{1}{y_t}\left(1-{\mathrm{e}}^{-y_t}\right),$(17) ) for $h\left(y^{\ast }\right)$. We want to write the condition solely in terms of the parameters, r and b, and find the curve in the r-b plane where stability changes. Because of the transcendental nature of the inequality, we will write this curve parametrically with r and b as functions of y. To do so, we first consider (A57) $r\left(\frac{{\mathrm{e}}^y-1-y{\mathrm{e}}^y}{y{\mathrm{e}}^y-y}\right)+\frac{1}{x}=0,$(A57) where we use x and y rather than $x^{\ast }$ and $y^{\ast }$ as in Appendix A.12.

We now incorporate the host and parasitoid nullclines, Equations (18a(18a) $\begin{array}{rl}1& =u(x^{\ast },y^{\ast })=g\left(x^{\ast }\right)[1-y^{\ast }h(y^{\ast }\left)\right],\end{array}$(18a) ) and (18b(18b) $\begin{array}{rl}1& =v(x^{\ast },y^{\ast })=b{x}^{\ast }g\left(x^{\ast }\right)h\left(y^{\ast }\right).\end{array}$(18b) ). For Model 4, these equations simplify to (A58) $r-rx-y=0,$(A58) and (A59) $bx{\mathrm{e}}^{r-rx}\frac{1}{y}\left(1-{\mathrm{e}}^{-y}\right)=1.$(A59) Using Equation (A58(A58) $r-rx-y=0,$(A58) ), we now eliminate x from Equation (A57(A57) $r\left(\frac{{\mathrm{e}}^y-1-y{\mathrm{e}}^y}{y{\mathrm{e}}^y-y}\right)+\frac{1}{x}=0,$(A57) ) and write r as a function of y, which simplifies to (A60) $r=\frac{y^2{\mathrm{e}}^y}{1+y{\mathrm{e}}^y-{\mathrm{e}}^y}.$(A60) We now return to Equation (A59(A59) $bx{\mathrm{e}}^{r-rx}\frac{1}{y}\left(1-{\mathrm{e}}^{-y}\right)=1.$(A59) ) and again eliminate x. We can then write b as a function of y, using Equation (A60(A60) $r=\frac{y^2{\mathrm{e}}^y}{1+y{\mathrm{e}}^y-{\mathrm{e}}^y}.$(A60) ) to eliminate r. After algebraic simplification, we obtain (A61) $b=\frac{y^2{\mathrm{e}}^y}{({\mathrm{e}}^y-1{)}^2}.$(A61) Equations (A60(A60) $r=\frac{y^2{\mathrm{e}}^y}{1+y{\mathrm{e}}^y-{\mathrm{e}}^y}.$(A60) ) and (A61(A61) $b=\frac{y^2{\mathrm{e}}^y}{({\mathrm{e}}^y-1{)}^2}.$(A61) ) give the boundary of the region in parameter space where the coexistence equilibrium satisfies the first Jury condition. We numerically verified that the first Jury condition is satisfied above this curve and violated below it. Since this curve is just barely below the curve for the second Jury condition found in Section A.14, this curve does not contribute to the stability region shown in Figure 3(d).

Now consider b = 1 with 0<r<2. As discussed previously, there is no equilibrium in the interior of the first quadrant for this case because the coexistence equilibrium has collided with the exclusion equilibrium point at $(1,0)$. Since y = 0, the first Jury condition is violated for b = 1, 0<r<2. Note also that when we take the limit as $y\to 0$ for Equations (A60(A60) $r=\frac{y^2{\mathrm{e}}^y}{1+y{\mathrm{e}}^y-{\mathrm{e}}^y}.$(A60) ) and (A61(A61) $b=\frac{y^2{\mathrm{e}}^y}{({\mathrm{e}}^y-1{)}^2}.$(A61) ), we obtain $(r,b)=(2,1)$. This means that the first Jury condition curve described by Equations (A60(A60) $r=\frac{y^2{\mathrm{e}}^y}{1+y{\mathrm{e}}^y-{\mathrm{e}}^y}.$(A60) ) and (A61(A61) $b=\frac{y^2{\mathrm{e}}^y}{({\mathrm{e}}^y-1{)}^2}.$(A61) ) connects to the first Jury condition curve given by b = 1, 0<r<2. Finally, when r = 0, system (49(49a) $\begin{array}{rl}{x}_{t+1}& =x_t{\mathrm{e}}^{r(1-x_t)}{\mathrm{e}}^{-y_t},\end{array}$(49a) ) has a line of equilibria on the x-axis. Since each of these equilibrium points has y = 0, the first Jury curve is violated for r = 0.

E.2. Second Jury condition

We use partial derivatives from Table 1 and Equation (17(17) $h\left(y_t\right)=\frac{1}{y_t}\left(1-{\mathrm{e}}^{-y_t}\right),$(17) ) to write the second Jury condition, inequality (37(37) $4+2x^{\ast }{u}_x+2y^{\ast }{v}_y+x^{\ast }{y}^{\ast }(u_x{v}_y-u_y{v}_x)>0.$(37) ), as (A62) $\begin{array}{r}4-2r{x}^{\ast }+2y^{\ast }\left[\frac{y^{\ast }-{\mathrm{e}}^{y^{\ast }}+1}{y^{\ast }({\mathrm{e}}^{y^{\ast }}-1)}\right]+x^{\ast }{y}^{\ast }\left[\frac{-r(y^{\ast }-{\mathrm{e}}^{y^{\ast }}+1)}{y^{\ast }({\mathrm{e}}^{y^{\ast }}-1)}\right]-x^{\ast }{y}^{\ast }\left[-1\left(\frac{1}{x^{\ast }}-r\right)\right]>0.\end{array}$(A62) Our goal is now to determine a curve in parameter space where stability changes. We re-write inequality (A62(A62) $\begin{array}{r}4-2r{x}^{\ast }+2y^{\ast }\left[\frac{y^{\ast }-{\mathrm{e}}^{y^{\ast }}+1}{y^{\ast }({\mathrm{e}}^{y^{\ast }}-1)}\right]+x^{\ast }{y}^{\ast }\left[\frac{-r(y^{\ast }-{\mathrm{e}}^{y^{\ast }}+1)}{y^{\ast }({\mathrm{e}}^{y^{\ast }}-1)}\right]-x^{\ast }{y}^{\ast }\left[-1\left(\frac{1}{x^{\ast }}-r\right)\right]>0.\end{array}$(A62) ) as a equality and eliminate x using the expression from the host nullcline given in Equation (A58(A58) $r-rx-y=0,$(A58) ). When we simplify and solve for r, we obtain (A63) $r=\frac{2{\mathrm{e}}^y-2+2y{\mathrm{e}}^y+y^2{\mathrm{e}}^y}{{\mathrm{e}}^y-1+y{\mathrm{e}}^y},$(A63) where y>0 functions as a parameter.

To get b as a function of y, we first use Equation (A58(A58) $r-rx-y=0,$(A58) ) to eliminate x from the parasitoid nullcline, Equation (A59(A59) $bx{\mathrm{e}}^{r-rx}\frac{1}{y}\left(1-{\mathrm{e}}^{-y}\right)=1.$(A59) ). Then, we use Equation (A63(A63) $r=\frac{2{\mathrm{e}}^y-2+2y{\mathrm{e}}^y+y^2{\mathrm{e}}^y}{{\mathrm{e}}^y-1+y{\mathrm{e}}^y},$(A63) ) to eliminate r. We solve for b and obtain (A64) $b=\frac{2y({\mathrm{e}}^y-1)+y^2{\mathrm{e}}^y(2+y)}{(2+y){\mathrm{e}}^{2y}-4{\mathrm{e}}^y+2-y},$(A64) where y>0.

Equations (A63(A63) $r=\frac{2{\mathrm{e}}^y-2+2y{\mathrm{e}}^y+y^2{\mathrm{e}}^y}{{\mathrm{e}}^y-1+y{\mathrm{e}}^y},$(A63) ) and (A64(A64) $b=\frac{2y({\mathrm{e}}^y-1)+y^2{\mathrm{e}}^y(2+y)}{(2+y){\mathrm{e}}^{2y}-4{\mathrm{e}}^y+2-y},$(A64) ) give the boundary of the region in parameter space where the coexistence equilibrium satisfies the second Jury condition. Note that when we take the limit as $y\to 0$ for Equations (A63(A63) $r=\frac{2{\mathrm{e}}^y-2+2y{\mathrm{e}}^y+y^2{\mathrm{e}}^y}{{\mathrm{e}}^y-1+y{\mathrm{e}}^y},$(A63) ) and (A64(A64) $b=\frac{2y({\mathrm{e}}^y-1)+y^2{\mathrm{e}}^y(2+y)}{(2+y){\mathrm{e}}^{2y}-4{\mathrm{e}}^y+2-y},$(A64) ), we obtain $(r,b)=(2,1)$. For b = 1, inequality (A62(A62) $\begin{array}{r}4-2r{x}^{\ast }+2y^{\ast }\left[\frac{y^{\ast }-{\mathrm{e}}^{y^{\ast }}+1}{y^{\ast }({\mathrm{e}}^{y^{\ast }}-1)}\right]+x^{\ast }{y}^{\ast }\left[\frac{-r(y^{\ast }-{\mathrm{e}}^{y^{\ast }}+1)}{y^{\ast }({\mathrm{e}}^{y^{\ast }}-1)}\right]-x^{\ast }{y}^{\ast }\left[-1\left(\frac{1}{x^{\ast }}-r\right)\right]>0.\end{array}$(A62) ) requires r<2. The point $(r,b)=(2,1)$ is where the Jury 2 curve intersects the b = 1 line. Thus, at the point $(r,b)=(2,1)$, the second Jury curve given parametrically by Equations (A63(A63) $r=\frac{2{\mathrm{e}}^y-2+2y{\mathrm{e}}^y+y^2{\mathrm{e}}^y}{{\mathrm{e}}^y-1+y{\mathrm{e}}^y},$(A63) ) and (A64(A64) $b=\frac{2y({\mathrm{e}}^y-1)+y^2{\mathrm{e}}^y(2+y)}{(2+y){\mathrm{e}}^{2y}-4{\mathrm{e}}^y+2-y},$(A64) ) intersects the first Jury curve, which is described in Section A.13. We again used numerical tests to verify that the second Jury condition is satisfied above the curve given by Equations (A63(A63) $r=\frac{2{\mathrm{e}}^y-2+2y{\mathrm{e}}^y+y^2{\mathrm{e}}^y}{{\mathrm{e}}^y-1+y{\mathrm{e}}^y},$(A63) ) and (A64(A64) $b=\frac{2y({\mathrm{e}}^y-1)+y^2{\mathrm{e}}^y(2+y)}{(2+y){\mathrm{e}}^{2y}-4{\mathrm{e}}^y+2-y},$(A64) ) and violated below it.

E.3. Third Jury condition

We use partial derivatives from Table 1 to write the third Jury condition, inequality (38(38) $1+x^{\ast }{u}_x+y^{\ast }{v}_y+x^{\ast }{y}^{\ast }(u_x{v}_y-u_y{v}_x)<1.$(38) ), as (A65) $\begin{array}{r}1-r{x}^{\ast }+\frac{y^{\ast }}{y^{\ast }h\left(y^{\ast }\right)}[1-y^{\ast }h(y^{\ast })-h(y^{\ast }\left)\right]-\frac{r{x}^{\ast }{y}^{\ast }}{y^{\ast }h\left(y^{\ast }\right)}[1-y^{\ast }h(y^{\ast })-h(y^{\ast }\left)\right]+x^{\ast }{y}^{\ast }\left(\frac{1}{x^{\ast }}-r\right)<1.\end{array}$(A65) This simplifies to (A66) $\frac{y^{\ast }{\mathrm{e}}^{y^{\ast }}}{{\mathrm{e}}^{y^{\ast }}-1}(1-x^{\ast }r)<1,$(A66) where we use Equation (17(17) $h\left(y_t\right)=\frac{1}{y_t}\left(1-{\mathrm{e}}^{-y_t}\right),$(17) ) for $h\left(y\right)$.

To find the curve where stability of the equilibrium changes, we again consider an equation instead of the inequality. After eliminating x using Equation (A58(A58) $r-rx-y=0,$(A58) ) from the host isocline, we solve for r as a function of y, (A67) $r=\frac{y{\mathrm{e}}^y+y^2{\mathrm{e}}^y-{\mathrm{e}}^y+1}{y{\mathrm{e}}^y}.$(A67) As was done in Section A.14, we use the parasitoid nullcline, Equation (A59(A59) $bx{\mathrm{e}}^{r-rx}\frac{1}{y}\left(1-{\mathrm{e}}^{-y}\right)=1.$(A59) ), with Equation (A67(A67) $r=\frac{y{\mathrm{e}}^y+y^2{\mathrm{e}}^y-{\mathrm{e}}^y+1}{y{\mathrm{e}}^y}.$(A67) ) to write b as a function of y, (A68) $b=\frac{y^3{\mathrm{e}}^y+y^2{\mathrm{e}}^y-y{\mathrm{e}}^y+y}{({\mathrm{e}}^y-1)(y{\mathrm{e}}^y-{\mathrm{e}}^y+1)},$(A68) noting again that y functions as a parameter used to trace out a curve in r–b space.

Equations (A67(A67) $r=\frac{y{\mathrm{e}}^y+y^2{\mathrm{e}}^y-{\mathrm{e}}^y+1}{y{\mathrm{e}}^y}.$(A67) ) and (A68(A68) $b=\frac{y^3{\mathrm{e}}^y+y^2{\mathrm{e}}^y-y{\mathrm{e}}^y+y}{({\mathrm{e}}^y-1)(y{\mathrm{e}}^y-{\mathrm{e}}^y+1)},$(A68) ) give the boundary of the region in parameter space where the coexistence equilibrium satisfies the third Jury condition. We numerically verified that the third Jury condition is satisfied above the curve and violated below it.

Additionally, when r = 0, the x-axis is a line of equilibrium points, as stated in Section A.13. Under these conditions, the expression for the determinant simplifies to $\mathrm{\Delta }=1$, and so the third Jury condition is also violated for r = 0 ($R_0=1$).