Dynamics of a discrete-time pioneer–climax model

Abstract

We present a simple mathematical model for the dynamics of a successional pioneer–climax system using difference equations. Each population is subject to inter- and intraspecific competition; population growth is dependent on the combined densities of both species. Nine different geometric cases, corresponding to different orientations of the zero-growth isoclines, are possible for this system. We fully characterize the long-term dynamics of the model for each of the nine cases, uncovering diverse sets of potential behaviors. Competitive exclusion of the pioneer species and of the climax species are both possible depending on the relative strength of competition. Stable coexistence of both species may also occur; in two cases, a coexistence state is destabilized through a Neimark–Sacker bifurcation, and an attracting invariant circle is born. The invariant circle eventually disappears into thin air in a heteroclinic or homoclinic bifurcation, leading to the sudden transition of the system to an exclusion state. Neither global bifurcation has been observed in a discrete-time pioneer–climax model before. The homoclinic bifurcation is novel to all pioneer–climax models. We conclude by discussing the ecological implications of our results.

Introduction

In ecological succession, newly established habitat is initially composed of pioneer plants, hardy species that do best at low population densities and that can colonize unpopulated environments (Dalling 2019; Harper 1977; Turner 2001). Later in succession come climax species, which are strong competitors but poor colonizers (Harper 1977; Turner 2001). At low population densities, climax plants typically cannot persist on their own, and pioneer species act as nurse plants in a commensal interaction (Harper 1977; Selgrade and Namkoong 1990). As population density of the climax species increases, interaction dynamics become more competitive and the climax species tends to thrive (Pandolfi 2019; Selgrade and Namkoong 1990). In nature, there are many pairs of plants that interact in this manner. White bursage and creosote bush from the Sonoran Desert, a pioneer species and a climax species, are one such pair (McAuliffe 1988).

Modeling pioneer–climax dynamics with systems of both differential and difference equations has a rich history (Buchanan and Selgrade 1995; Franke and Yakubu 1994; Franke and Yakubu 1994; Kim and Marlin 1999; Selgrade and Namkoong 1990; Selgrade and Roberds 1996; Selgrade and Roberds 1997; Sumner 1996). In these models, the long-term behavior of both species depends on inter- and intraspecific competition, and each species reacts differently to these competitive effects. Most pioneer–climax models distinguish density-dependent effects from competition using a total-density variable, a linear combination of the population densities of both pioneer and climax species weighted by competitive effects. Competitive effects are thus dependent only on individual population densities. The models then assume that the per-capita replacement rate, or fitness, of each species is a function of the total-density variable, so that the species’ response to density is modeled by this fitness function.

A pioneer species’ fitness function is a monotonically decreasing function of total density (Selgrade and Namkoong 1990). A climax species, on the other hand, does poorly at low densities. As total density increases, its fitness increases, until overcrowding at high densities results in a decrease in fitness (Selgrade and Namkoong 1990), forming a one-humped function. Such fitness functions match ecological observations, incorporating the dependence of the climax species on the pioneer at low densities for resources such as shade or shelter (Harper 1977; McAuliffe 1988). This formulation also includes the pioneer species being generally outcompeted by the climax species as density increases (McAuliffe 1988).

Initially, most researchers assumed that the eventual result of any pioneer–climax dynamical model would be pioneer exclusion, an assumption supported by model results and ecological theory (Franke and Yakubu 1994; Harper 1977; Selgrade and Namkoong 1990). Habitat disturbances that prevent a stable state with only climax species are likely to occur (Dalling 2019; Tilman 1988), but scientists believed an undisturbed system would always result in climax species dominance. Selgrade and Namkoong (1990), however, demonstrated the possibility of long-term dynamics other than pioneer exclusion by finding a Hopf bifurcation that led to stable periodic behavior and coexistence.

After this seminal work, more evidence for a variety of steady-state behaviors emerged (Buchanan 1999; Franke and Yakubu 1994; Franke and Yakubu 1994; Kim and Marlin 1999). Much of this research used differential equations (Buchanan and Selgrade 1995; Cornejo-Pérez and Barradas 2015; Kim and Marlin 1999; Selgrade 1994; Sumner 1996). Studies that employed difference equations generally focused on analysis of the possible Neimark–Sacker bifurcation, the discrete-time equivalent of a Hopf bifurcation, or theoretical proof of exclusion and coexistence principles (Franke and Yakubu 1994; Franke and Yakubu 1994; Selgrade and Namkoong 1990; Selgrade and Namkoong 1991; Selgrade and Roberds 1996). From these initial studies, pioneer–climax research turned to studying how a forcing term, representing stocking or harvesting, might be able to restore stability in the system after a Hopf or period-doubling bifurcation has occurred (Buchanan and Selgrade 1995; Selgrade 1994; Selgrade 1998; Selgrade and Roberds 1996; Selgrade and Roberds 1998; Sumner 1996; Sumner 1998). To this point, no one has undertaken a more general overview of the behavior of a simple pioneer–climax difference-equation model.

We aim to provide such an overview, in order to understand the full long-term dynamics of a discrete-time model for the first time. With this comprehensive analysis, we gain further insight into how difference-equation models may contrast with differential-equation models, and we shed light on the unique behaviors a discrete-time model may reveal.

We present a basic pioneer–climax difference-equation model using similar assumptions to prior studies. In Sect. 2 (Model), we describe the model and its development. In this section, we also nondimensionalize our model. Using a total-density variable, we separate inter- and intraspecific competition effects and model the fitness of each species as a function of the total-density variable. We use simple rational forms, for both fitness functions, that are constructed in keeping with earlier model assumptions and ecological theory.

In Sect. 3 (Zero-growth isoclines and equilibria), we characterize the zero-growth isoclines that provide the framework for our stability analysis. We also present the fixed points, or equilibria, of our model, which occur at the intersection of the zero-growth isoclines. Nine different geometric configurations of the zero-growth isoclines occur, as in Buchanan's (1999) analysis using differential equations. These nine cases present scenarios of both competitive exclusion and coexistence.

We perform a steady-state analysis of the long-term behavior of the model in Sect. 4 (Stability analysis) and summarize the results of the stability analysis for all nine geometric cases. The full stability analysis is contained in Appendix A (Preliminary organization and derivatives for the Jacobian), Appendix B (Stability of boundary fixed points), and Appendix C (Stability of interior fixed points).

In two of the above cases, Neimark–Sacker bifurcations occur. We examine these bifurcations in Sect. 5 (Neimark-Sacker bifurcation). In Sect. 6 (Global bifurcations), we highlight a global bifurcation similar to the heteroclinic cycle found in the differential-equation model of Kim and Marlin (1999), as well as another global bifurcation, a homoclinic bifurcation, not previously found in pioneer–climax models. We also present numerical simulations that show the changing stabilities of the fixed points or limit cycles as the bifurcations occur.

In Sect. 7 (Discussion), we describe both the implications of our findings and how our work motivates future research. Our results provide insight into the persistence of pioneer and climax plant species and highlight previously unobserved global behavior of discrete-time pioneer–climax models. The investigation presented here also opens the door to new research, setting the stage for examining pioneer–climax models in a spatial context. In particular, this framework allows us to examine the invasion ecology of a spatial pioneer–climax model with flexible dispersal kernels and distinct growth and dispersal stages. As climate change alters ecosystems and induces range shifts in many species, we may see more successional communities, making such analyses highly valuable.

Model

Model development

We consider a model where sufficient resources exist that the per-capita replacement rates of the pioneer and climax species depend only on total-density variables, linear combinations of pioneer and climax densities weighted by competition coefficients. This construction allows us to distinguish between a species’ response to density and its response to competitive effects. Let ${P_t}$ and ${C_t}$ be the population densities of the pioneer and climax species at time t. The total-density variables for each species are

$${M_t} = {w_{11}} \ {P_t} + {w_{12}} \ {C_t},$$

(1)

$${N_t} = {w_{21} \ {P_t}} + {w_{22}} \ {C_t},$$

(2)

where $w_{ij}$ denotes the competitive effect of species j on species i. We take species 1 to be the pioneer species and species 2 to be the climax species.

The interaction between species is modeled by the system of difference equations

$$P_{t+1} = f(M_t) \ {P_t},$$

(3)

$$C_{t+1} = g(N_t) \ {C_t},$$

(4)

where $f({M_t})$ and $g({N_t})$ are the per-capita replacement rates, or fitness functions. We omit consideration of survivorship from one generation to the next (but see Sect. 7, Discussion). Note that for population growth of either species in this discrete-time model its fitness function $f({M_t})$ or $g({N_t})$ must be greater than 1.

We assume that the pioneer species’ fitness function $f({M_t})$ is a monotonically decreasing function with respect to density that has a single root $f(M_t) = 1$. To model the pioneer species’ dynamics under this assumption, consider the fitness function

$$f(M_t) = \frac{R_1}{1 + aM_t}.$$

(5)

The net reproductive rate of the pioneer species, ${R_1}$, is assumed to be greater than one to reflect the pioneer species’ ability to survive on its own at low densities. We group deleterious effects such as overcrowding into the positive parameter a, so that the pioneer species’ fitness always decreases with increasing total density, as shown in Fig. 1a.

The climax species’ fitness function $g(N_t)$ is a one-humped function of density with two roots $g(N_t) = 1$. Here, we take

$$\begin{aligned} g(N_t) = \frac{{R_2} (1 + {b_1} {N_t})}{1 + {b_2} {N_t}^2}. \end{aligned}$$

(6)

As with the pioneer species, ${{R_2}}$ is the climax species’ net reproductive rate. As the climax species require the pioneer species to facilitate growth at low densities, ${{R_2}}$ is assumed to be less than one. Beneficial effects from the pioneer species acting as nurse plant for the climax species at low densities are incorporated in the positive parameter ${{b_1}}$. The positive parameter ${{b_2}}$, multiplied by the square of the total density, measures deleterious effects on the population at high densities due to effects such as overcrowding. This arrangement gives the desired form of the climax fitness function, increasing at low densities before hitting a maximum and decreasing as total density grows (Fig. 1b). Thus, the fitness functions presented here match ecological theory on the per-capita replacement rates of pioneer and climax species (Harper 1977; McAuliffe 1988) as well as mathematical convention (Kim and Marlin 1999; Selgrade and Namkoong 1990) while maintaining a desirable simplicity.

Fig. 1

Fitness functions for both species from the dimensional model. a, the pioneer species’ fitness function $f(M_t)$. b, the climax species’ fitness function $g(N_t)$. The horizontal dashed gray lines indicate where the fitness function is equal to one, where the root(s) of the functions occur. Parameter values are ${R_1} = 1.2, {{R_2}} = 0.8, a = 0.15, {{b_1}} = 0.6$, ${{b_2}} = 0.15, {w_{11}} = 1, {w_{12}} = 2, {w_{21}} = 1.2, {w_{22}} = 1$. Note that these parameter values correspond to the nondimensionalized model presented below with ${w_{12}} = {{v_2}}$ and ${w_{21}} = {{v_1}}$

Nondimensionalization

We scale the above system with respect to intraspecific competition, letting

$$\begin{aligned} {x_t} = {w_{11}} \ {P_t},\end{aligned}$$

(7)

$$\begin{aligned} {y_t} = {w_{22}} \ {C_t}. \end{aligned}$$

(8)

This rescaling results in the reformulated total-density variables

$$\begin{aligned} {m_t} = {x_t} + \left( \frac{w_{12}}{w_{22}}\right) {y_t} = {x_t} + {{v_2}} {y_t},\end{aligned}$$

(9)

$$\begin{aligned} {n_t} = \left( \frac{w_{21}}{w_{11}}\right) {x_t} + {y_t} = {{v_1}} {x_t} + {y_t}, \end{aligned}$$

(10)

where ${{v_1}}$ and ${{v_2}}$ are now relative competition coefficients measuring the strength of interspecific competition proportionate to intraspecific competition.

We now have the nondimensional model

$$\begin{aligned} x_{t+1} = f(m_t) \ {x_t} = \frac{{R_1} {x_t}}{1 + {am_t}},\end{aligned}$$

(11)

$$\begin{aligned} y_{t+1} = g(n_t) \ {y_t} = \frac{{{R_2}} (1 + {{b_1}} {n_t}) {y_t}}{1 + {{b_2}} {n_t}^2}. \end{aligned}$$

(12)

This scaling creates a model with fewer parameters, simplifying analysis significantly.

Zero-growth isoclines and equilibria

To organize our analysis of the model given by Eqs. 11 and 12, we use zero-growth isoclines as our framework. A zero-growth isocline is the set of points at which the population density of a species is not changing. Thus, in our discrete-time system, the zero-growth isoclines of the pioneer or climax species occur when $x_{t+1} = {x_t} = x$ or $y_{t+1} = {y_t} = y$, respectively. Under these conditions, our total density variables are ${m_t} = m = x + {{v_2}} y$ and ${n_t} = n = {{v_1}} x + y$. This substitution gives us the equations

$$\begin{aligned} x = f(m) \ x,\end{aligned}$$

(13)

$$\begin{aligned} y = g(n) \ y. \end{aligned}$$

(14)

We see that Eq. 13 is satisfied either when $x = 0$ or when $f(m) = 1$. By construction, there is one root of the pioneer species’ fitness function at $f(m) = 1$. Thus, in total, there are two pioneer species zero-growth isoclines, the trivial one $x = 0$ and a nontrivial isocline where $f(m) = 1$ in the (x, y) plane. Both isoclines are linear.

Similarly, Eq. 14 has the solutions $y = 0$ or $g(n) = 1$. As the climax species’ fitness function has two roots, we therefore have three zero-growth isoclines for the climax species. These isoclines are the trivial isocline $y = 0$ and two nontrivial isoclines from the solutions to $g(n) = 1$ in the (x, y) plane. The climax species’ zero-growth isoclines are also linear.

Explicitly solving for the pioneer species’ zero-growth isocline given by $f(m) = 1$ in terms of x and y yields the line

$$\begin{aligned} x = {{m^*}} - {{v_2}} y, \end{aligned}$$

(15)

where

$$\begin{aligned} {{m^*}} = \frac{{R_1} - 1}{a}. \end{aligned}$$

(16)

Doing the same for the climax species’ zero-growth isoclines given by $g(n) = 1$ gives us the two lines

$$\begin{aligned} y = {{n_1}}^* - {v_1} x,\end{aligned}$$

(17)

$$\begin{aligned} y = {{n_2}}^* - {{v_1}} x, \end{aligned}$$

(18)

where

$$\begin{aligned} n_{1,2}^* = \frac{{{R_2}} {{b_1}} \pm \sqrt{{{R_2}}^2 {{b_1}}^2 - {4{b_2}}\left( 1 - {{R_2}}\right) }}{2 {{b_2}}}, \end{aligned}$$

(19)

with ${{n_1}}^*$ corresponding to the negative square root and ${{n_2}}^*$ the positive square root.

A fixed point of Eqs. 11 and 12 occurs at the intersection of two zero-growth isoclines, one isocline from each species, such that both populations are not changing over time. Given two zero-growth isoclines for the pioneer species and three for the climax species, there are six possible intersection points, and therefore six possible fixed points. We will label these fixed points as ${E_i} = \left( {x_i}, {y_i}\right)$ where the subscript denotes the number of the fixed point instead of the time step.

Four of these six, including the origin, are boundary equilibria and always exist. These boundary, or exclusion, equilibria are

$$\begin{aligned} {{E_1}} = \left( {x_1}, {y_1} \right) = \left( 0,0\right) , \end{aligned}$$

(20)

the equilibrium at the origin where neither species persists,

$$\begin{aligned} {{E_2}} = \left( {x_2}, {y_2} \right) = \left( {{m^*}}, \ 0\right) , \end{aligned}$$

(21)

the equilibrium where the pioneer species persists and excludes the climax species,

$$\begin{aligned} {{E_3}} = \left( {x_3}, {{y_3}} \right) = \left( 0, \ {{n_1}}^*\right) , \end{aligned}$$

(22)

the lower equilibrium on the y-axis where the climax species persists but the pioneer species does not, and

$$\begin{aligned} {{E_4}} = \left( {x_4}, {{y_4}} \right) = \left( 0, \ {{{n_2}}^*}\right) , \end{aligned}$$

(23)

the upper equilibrium on the y-axis where the climax species persists. The boundary equilibria result from the intersection of a trivial isocline ($x = 0$ or $y = 0$) with a zero-growth isocline of the other species.

Up to two interior equilibria may also exist. An interior coexistence state occurs when the nontrivial pioneer species’ isocline crosses one of the nontrivial climax species’ isoclines, i.e., when both $f(m) = 1$ and $g(n) = 1$. These equilibria are

$$\begin{aligned} {{E_5}} = ({{x_5}}, {{y_5}}) = \left( \frac{{{m^*}} - {{v_2}} {{n_1}}^*}{1 - {{v_1}} {{v_2}}}, \frac{{{n_1}}^* - {{v_1}} {{m^*}}}{1 - {{v_1}} {{v_2}}}\right) , \end{aligned}$$

(24)

which results from the pioneer species’ nontrivial zero-growth isocline crossing the lower climax species’ nontrivial zero-growth isocline, and

$${{E_6}} = ({{x_6}}, {{y_6}}) = \left( \frac{{{m^*}} - {{v_2}} {{{n_2}}^*}}{1 - {{v_1}} {{v_2}}}, \frac{{{{n_2}}^*} - {{v_1}} {{m^*}}}{1 - {{v_1}} {{v_2}}}\right),$$

(25)

where the pioneer species’ nontrivial zero-growth isocline crosses the upper climax species’ nontrivial zero-growth isocline.

In order to analyze the equilibria and understand the long-term dynamics of our model (Eqs. 11 and 12), we must consider the orientations of the isocline intersections that give rise to each equilibrium point. There are nine distinct geometric configurations of the zero-growth isoclines, some with interior equilibria and some without (Figs. 2 and 3). The cases are distinguished from each other by restrictions on the parameter values. In particular, we fix all parameters but the relative competition coefficients ${{v_1}}$ and ${{v_2}}$, and obtain limits on ${{v_1}}$ and ${{v_2}}$ in terms of the other parameters for each case.

In three cases, the nontrivial isoclines do not cross each other and so only the four boundary equilibria ${E_{1-4}}$ exist. These cases are shown in Fig. 2. In the six remaining cases, the nontrivial pioneer species’ isocline crosses one or both of the nontrivial climax species’ isoclines. The six configurations with interior equilibria are presented in Fig. 3.

Fig. 2

Zero-growth isocline orientations and stabilities of fixed points for cases where only boundary fixed points exist. The dashed-dotted lines are the pioneer species’ zero-growth isoclines; solid lines (including $y = 0$) are the climax species’ zero growth isoclines. Solid black circles indicate asymptotically stable fixed points. Circles with crosses are saddle points, while open circles signify other unstable fixed points. For all cases, parameter values are ${R_1} = 1.2, {{R_2}} = 0.8, a = 0.3, {{b_1}} = 0.6, {{b_2}} = 0.2$. a, ${{v_1}} = 0.35, {{v_2}} = 2.05$. b, ${{v_1}} = 1.2, {{v_2}} = 0.5$. c, ${{v_1}} = 4.5, {{v_2}} = 0.275$

Fig. 3

Zero-growth isocline configurations and stabilities of equilibria for cases where interior equilibria occur. The dashed-dotted lines are the pioneer species’ zero-growth isoclines; solid lines (including $y = 0$) are the climax species’ zero growth isoclines. Parameter values for all cases are ${R_1} = 1.2, {R_2} = 0.8, a = 0.3, {b_1} = 0.6, {b_2} = 0.2$. The same classification system as in Fig. 2 applies; gray shaded circles represent equilibria that undergo bifurcations. a and b, one interior equilibrium results from the crossing of the pioneer zero-growth isocline and the lower climax zero-growth isocline. a, ${{v_1}} = 1.2, {{v_2}} = 1.75$; b, ${{v_1}} = 0.5, {{v_2}} = 0.8$. c and d, one interior equilibrium results from crossing the upper climax zero-growth isocline. c, ${{v_1}} = 4.2, {{v_2}} = 0.6$; d, ${{v_1}} = 2, {{v_2}} = 0.3$. e and f, two interior equilibria occur. e, ${{v_1}} = 4, {{v_2}} = 2$; f, ${{v_1}} = 0.5, {{v_2}} = 0.23$

Stability analysis

For all nine cases, we now determine the stabilities of the four to six equilibria that occur in each case. To assess stability of an equilibrium, we can examine the eigenvalues of the Jacobian of Eqs. 11 and 12 evaluated at that equilibrium. If both eigenvalues are less than one in magnitude, i.e., $|{{\lambda _1}}| < 1$ and $|{{\lambda _2}}| < 1$, we have asymptotic stability. The general form of the Jacobian matrix for Eqs. 11 and 12 is

$$\begin{aligned} J &{}= \begin{bmatrix} f(m_t) + {x_t} \ f'(m_t) \frac{\partial m}{\partial x} &{} {x_t} \ f'(m_t) \frac{\partial m}{\partial y} \\ {y_t} \ g'(n_t) \frac{\partial n}{\partial x} &{} g(n_t) + {y_t} \ g'(n_t) \frac{\partial n}{\partial y} \end{bmatrix} \\ &{}= \begin{bmatrix} f(m_t) + {x_t} \ f'(m_t) &{} {{v_2}} \ {x_t} \ f'(m_t) \\ {{v_1}} \ {y_t} \ g'(n_t) &{} g(n_t) + {y_t} \ g'(n_t) \end{bmatrix}. \end{aligned}$$

(26)

Boundary equilibria

The eigenvalues of the Jacobian evaluated at each of the four boundary equilibria are given in Table 1. When $x = 0$ or ${y = 0}$, as in the boundary equilibria, an off-diagonal element of the Jacobian J is zero. This property allows us to simply read the eigenvalues off the diagonal of J. Additionally, for ${{E_2}}$ we may make the substitution $f(m) = 1$, as ${{E_2}}$ arises from the intersection of the isoclines $y = 0$ and $f(m) = 1$. Similarly, for ${E_{3,4}}$ we may say $g(n) = 1$, as these equilibria result from the intersections of $x = 0$ and the solutions to $g(n) = 1$.

The eigenvalues of all four boundary equilibria are real for any parameter values. We can also prove that the eigenvalues for all boundary equilibria are positive in all nine cases (see Appendix B, Stability of boundary fixed points, for details).

The origin is always a saddle point, as we assume ${R_1} > 1$ and ${{R_2}} < 1$. The stabilities of the other boundary equilibria are simple to determine in all cases by analyzing the eigenvalues. For ${E_{2-4}}$, we can establish whether ${{\lambda _1}}$ is greater or less than $+1$ by comparing where the non-zero component of the fixed point (i.e., ${x_2}$, ${{y_3}}$, or ${{y_4}}$) lands in relation to the root(s) where $g\left( {{v_1}} x\right) = 1$ or $f\left( {{v_2}} y\right) = 1$. This relation will be based on zero-growth isocline configuration for a given case. Similarly for ${E_{2-4}}$, ${{\lambda _2}}$ can be analyzed by determining the sign of the derivative of the fitness function evaluated at the given fixed point. If the derivative is positive, ${{\lambda _2}} > 1$, and if it is negative, ${{\lambda _2}} < 1$.

For the three cases in which only boundary equilibria exist, this section completes our stability analysis. These three cases and the stabilities of their equilibria are shown in Fig. 2. The stabilities of the boundary equilibria in the six cases where at least one interior fixed point exists are shown in Fig. 3. A full stability analysis for the boundary equilibria in all nine cases is presented in Appendix B (Stability of boundary fixed points).

Table 1 Boundary equilibria and their eigenvalues. Each equilibrium is generally represented by ${E_i} = \left( {x_i}, {y_i}\right)$

Interior equilibria

The remaining six cases, where at least one interior fixed point occurs, are of significant interest given the potential for stable coexistence states. We classify the stabilities of the boundary fixed points in these cases as above, by analyzing the magnitude of their eigenvalues. For the interior equilibria, however, analyzing the magnitude of the eigenvalues is difficult, as the form of the eigenvalues is complicated and the eigenvalues may be complex. Thus, to determine the stabilities of the interior equilibria, we use the Jury conditions (Jury 1964) applied to the Jacobian of Eqs. 11 and 12 evaluated at each equilibrium point.

For this system, three Jury conditions together prove stability of an equilibrium, given by

$$\begin{aligned} 1 - \tau + \Delta > 0,\end{aligned}$$

(27)

$$\begin{aligned} 1 + \tau + \Delta > 0,\end{aligned}$$

(28)

$$\begin{aligned} \Delta < 1, \end{aligned}$$

(29)

where $\tau$ and $\Delta$ are the trace and determinant of the Jacobian J. The first condition ensures there are no eigenvalues greater than $+1$, the second that there are no eigenvalues less than $-1$. The third condition guarantees that no complex eigenvalues lie outside the unit circle.

Violating any of the Jury conditions as parameter values change will result in the equilibrium point changing stability via a local bifurcation. Violating the first Jury condition leads to a transcritical bifurcation. Violating the second Jury condition results in a flip, or period-doubling, bifurcation. If the third condition is violated, a Neimark–Sacker bifurcation will occur. This bifurcation is the discrete-time analog of a Hopf bifurcation. By testing the Jury conditions on the Jacobian (Eq. 26) evaluated at the interior equilibrium point(s), for all six cases with coexistence states, we classify the stabilities of the interior equilibria as shown in Fig. 3.

For every case, the first Jury condition is either always satisfied or never satisfied for any parameter values within the case’s bounds. Thus, there are no transcritical bifurcations for either coexistence fixed point while that fixed point remains within the interior of the first quadrant. However, the interior fixed point will shift toward the x-axis if ${{v_1}}$ is varied and the y-axis if ${{v_2}}$ is varied. There are critical ${{v_1}}$ and ${{v_2}}$ values that distinguish each geometric case from the others, and precisely at these critical values the interior fixed point hits one of the axes, colliding with one of the boundary equilibria in the process. Therefore, transcritical bifurcations do occur as the interior fixed point passes through the boundary equilibria, corresponding to a change in geometric case. This bifurcation results in the coexistence equilibria exiting the first quadrant, leading to negative population densities.

The results of the transcritical bifurcations that occur are evident from Figs. 2 and 3. Compare, for example, Figs. 2c and 3c. As ${{v_2}}$ decreases, the pioneer species’ isocline (dashed-dotted line) in Fig. 3c shifts, its y-intercept increasing, and the interior equilibrium ${{E_5}}$ moves toward the y-axis. Eventually, ${{E_5}}$ collides with the y-axis at ${{E_4}}$, and the equilibria exchange stabilities. The previously stable equilibrium ${{E_4}}$ becomes a saddle point, and ${{E_5}}$ passes out of the first quadrant as an asymptotically stable fixed point. As ${{v_2}}$ decreases further, the pioneer species’ isocline lies fully above the two climax species’ isoclines, and we are now in the case shown in Fig. 2c, where we see that ${{E_4}}$ is indeed a saddle point and there is no interior fixed point in the first quadrant.

The second Jury condition is never violated, and so there are no period-doubling bifurcations. More precisely, we can prove, through use of Descartes’ Rule of Signs (Descartes 1886), that all real eigenvalues of the Jacobian (Eq. 26) are positive for both interior equilibria in all cases. This property rules out violation of the second Jury condition. This result is in contrast to several other studies that demonstrated period-doubling bifurcations in their discrete-time pioneer–climax models (Joshi and Blackmore 2010; Selgrade 1998; Selgrade and Roberds 1997; Selgrade and Roberds 1998).

For all four cases in which equilibrium ${{E_6}}$ exists (Figs. 3c−f), we can prove that all eigenvalues of the Jacobian evaluated at ${{E_6}}$ are strictly real. We can prove the same for two of the cases with equilibrium ${{E_5}}$ (Fig. 3b, f). Thus, we do not analyze the third Jury condition at all for ${{E_6}}$, nor for ${{E_5}}$ in two cases, as this condition applies to complex eigenvalues. We determine stabilities of the relevant interior equilibria in these cases by again using Descartes’ Rule of Signs (Descartes 1886), to find the number of eigenvalues that are above and below $+1$. As we also know that our real roots are positive, this information is all we require to classify the stabilities of the interior fixed point(s) in these cases.

In the remaining two cases in which ${{E_5}}$ occurs, one with only ${{E_5}}$ (Fig. 3a) and one with both interior equilibria (Fig. 3e), we cannot show that the eigenvalues at ${{E_5}}$ are always real. Numerically, we find that the eigenvalues may be real or complex for different parameter values. Thus, for these two cases we fix ${R_1}, {{R_2}}, a, {{b_1}}$, and ${{b_2}}$ and analyze the third Jury condition to find a region in the $({{v_1}}, {{v_2}})$ plane for which equilibrium ${{E_5}}$ is asymptotically stable. The parameters ${{v_1}}$ and ${{v_2}}$ were selected to highlight how relative competition coefficients affect stability. Competition is the main underlying factor influencing the density of both pioneer and climax species, as seen in the construction of our model. Therefore, the competition parameters were a natural choice to analyze.

Fig. 4

Stability region for equilibrium ${{E_5}}$, defined by the curves in $({{v_1}}, {{v_2}})$ parameter space where bifurcations occur. In both cases, ${{E_5}}$ is stable for parameter values above and to the left of the solid curve given by $\Delta = 1$; crossing this curve results in a loss of stability via a Neimark–Sacker bifurcation and the birth of an invariant circle. The lower dashed-dotted curves are the boundaries at which the global bifurcations occur; crossing the dashed-dotted curve results in the abrupt disappearance of the invariant circle. a, the first case corresponding to Fig. 3a with a heteroclinic bifurcation. The global bifurcation curve was found by using the bisection method to calculate, for a given ${{v_1}}$, the ${{v_2}}$ value for which a trajectory starting on the unstable manifold of ${{E_2}}$ ended at, or very near to, ${{E_3}}$, such that the unstable manifold of ${{E_2}}$ connected with the stable manifold of ${{E_3}}$. b, the second case corresponding to Fig. 3e with a homoclinic bifurcation. The global bifurcation curve for this case was found by using the bisection method to determine, for a given ${{v_1}}$, the ${{v_2}}$ value for which a trajectory starting on the unstable manifold of ${{E_6}}$ returned to ${{E_6}}$, such that the unstable and stable manifolds of ${{E_6}}$ connected to form a homoclinic connection. The dashed gray lines in each panel indicate the boundaries on ${{v_1}}$ and ${{v_2}}$ such that we are in the appropriate geometric case, i.e., that the isoclines intersect in the appropriate configuration. Parameter values for both cases are ${R_1} = 1.2, {{R_2}} = 0.8, a = 0.3, {{b_1}} = 0.6, {{b_2}} = 0.2$

Fig. 5

Growth of the invariant circle as ${{v_2}}$ decreases from 2.23 to just above 1.985 in the first case with a Neimark–Sacker bifurcation. The point in the middle of the invariant circles is ${{E_5}}$ for ${{v_2}} = 2.23$, when ${{E_5}}$ is asymptotically stable. We ran the system for 50,000 iterations; the final 15,000 iterations were plotted for ten ${{v_2}}$ values. Parameter values are ${R_1} = 1.2, {{R_2}} = 0.8, a = 0.3, {{b_1}} = 0.6, {{b_2}} = 0.2$

Fig. 6

Growth of the invariant circle in the second case as ${{v_2}}$ decreases from 5.2 to 2.9. The single point in the middle of the invariant circles is ${{E_5}}$ for ${{v_2}} = 5.2$, when ${E_5}$ is asymptotically stable. The system was run for 1,000,000 iterations; the final 200,000 iterations were plotted for ten ${{v_2}}$ values. This number of iterations was necessary to fill in the invariant circle for some ${{v_2}}$ values for which the system had nearly rational rotation numbers. Other parameter values are ${R_1} = 1.2, {{R_2}} = 0.8, a = 0.3, {{b_1}} = 0.6, {{b_2}} = 0.2$

The determinant of the Jacobian evaluated at ${{E_5}}$ is

$$\begin{aligned} \Delta = \left( 1 + {{x_5}} f'(m)\right) \left( 1 + {{y_5}} g'(n) \right) - {{v_1}} {{v_2}} {{x_5}} {{y_5}} f'(m) g'(n), \end{aligned}$$

(30)

with $m = {{x_5}} + {{v_2}} {{y_5}}$ and $n = {{v_1}} {{x_5}} + {{y_5}}$. After substituting in the explicit forms of the equilibrium and derivatives, it is possible to solve the equation $\Delta = 1$ for ${{v_2}}$ as a function of ${{v_1}}$ (or vice versa). Thus, for both cases we can get a curve for the third Jury condition that defines the region of stability for ${{E_5}}$ in the $({{v_1}}, {{v_2}})$ plane. These curves are shown in Fig. 4. Above and to the left of the solid curve, ${{E_5}}$ is asymptotically stable. Crossing the solid curve violates the third Jury condition and destabilizes ${{E_5}}$ via a Neimark–Sacker bifurcation. For the explicit form of the curves for the third Jury condition, as well as a full stability analysis of the interior fixed points, see Appendix C (Stability of interior fixed points).

In all, four of the six cases with interior equilibria have an asymptotically stable coexistence state for at least some parameter values. In two of these cases, equilibrium ${E_6}$ is always asymptotically stable. The other two have limited parameter regions where equilibrium ${{E_5}}$ is asymptotically stable (Fig. 4). Even after ${{E_5}}$ has become unstable in these two cases, an attracting invariant circle is born from the Neimark–Sacker bifurcation that destabilized ${{E_5}}$. We will now examine this local bifurcation in more detail.

Neimark–Sacker bifurcation

In two geometric cases we violate the third Jury condition for equilibrium ${{E_5}}$ by having a pair of complex conjugate eigenvalues exit the unit circle. These cases correspond to Fig. 3a, e. For the remainder of this section and the next (Sect. 6, Global bifurcations), we will refer to them as the first case and the second case.

In both cases, for parameter values above and to the left of the curve $\Delta = 1$, equilibrium ${{E_5}}$ is asymptotically stable (Fig. 4). As ${{v_2}}$ decreases and crosses this boundary, the equilibrium point loses stability via a Neimark–Sacker bifurcation and an attracting quasiperiodic solution, commonly known as an invariant circle, emerges. As ${{v_2}}$ continues to decrease, the invariant circle grows larger. The attracting invariant circle for a variety of ${{v_2}}$ values can be seen in Fig. 5 for the first case and in Fig. 6 for the second case.

As we have noted, loss of stability via a Neimark–Sacker bifurcation for fixed point ${{E_5}}$ occurs as the relative competition coefficient ${{v_2}}$ decreases for fixed ${{v_1}}$. Thus, as intraspecific competition in the climax species grows stronger relative to the competitive effect of the climax species on the pioneer species, the equilibrium destabilizes. Equivalently, if ${{v_2}}$ were fixed, the equilibrium loses stability as ${{v_1}}$ increases. Biologically, this situation corresponds to the competitive effect of the pioneer species on the climax species increasing relative to the competitive effect of the pioneer species on itself. In either context, ${{v_1}}$ increasing or ${{v_2}}$ decreasing, the fixed point ${{E_5}}$ is destabilized when a species’ competitive effect on the climax species grows stronger relative to its competitive effect on the pioneer species. In place of the fixed point, an attracting quasiperiodic solution arises, resulting in stable population fluctuations.

To fully characterize the Neimark–Sacker bifurcation, we must consider the possibility of resonance. Resonance cases occur when the complex conjugate eigenvalues $\lambda$, $\overline{\lambda }$ of our Jacobian exit the unit circle through a root of unity such that $\lambda = e^{2\pi i p/q}$ and $\overline{\lambda } = e^{-2\pi i p/q}$, where p and q are both positive integers and p/q is known as the rotation number (Aronson et al. 1982; Whitley 1983). This phenomenon is called strong resonance if $q \le 4$ and weak resonance for $q \ge 5$ (Lauwerier 1986; Whitley 1983). The bifurcation that occurs in a resonance case does not lead to the birth of an invariant circle; phase-locked periodic solutions or more complicated behaviors arise instead (Arnol’d 1977; Kuznetsov 2004; Lauwerier 1986).

In both cases with the Neimark–Sacker bifurcation, the eigenvalues of the Jacobian are precisely on the unit circle for any $({{v_1}}, {{v_2}})$ pair that lies on the curve $\Delta = 1$ (Fig. 4). Thus, if the eigenvalues are located at roots of unity for any such $({{v_1}}, {{v_2}})$ pair, resonance occurs for those parameter values. To check for cases of resonance, we numerically calculate the eigenvalues of the Jacobian for all $({{v_1}}, {{v_2}})$ pairs on the curve $\Delta = 1$, i.e., parameters for which the eigenvalues lie on the unit circle, and isolate the rotation number p/q. If the rotation number is rational, we have a p/q resonance (Aronson et al. 1982).

In both cases with a Neimark–Sacker bifurcation, the rotation numbers as our eigenvalues pass through the unit circle are very small. This small size precludes the possibility of strong resonance, as the smallest rotation number that would lead to strong resonance is $p/q = 1/4$.

There is no solid evidence for weak resonance either, although this absence could be due to numerical inaccuracy in calculating the rotation numbers. For some parameter values, particularly in the second case with a Neimark–Sacker bifurcation, it takes many thousands of iterations to densely fill the invariant circle. This situation indicates we may have rotation numbers that are very close to, but not quite, rational. Weak resonance cases give rise to phase-locked periodic orbits, but there is no evidence of periodic solutions in a bifurcation diagram. Instead, the bifurcation diagram shows only densely filled invariant circles or fixed points (Figs. 7 and 8). Thus, if phase-locking occurs, it is likely for an orbit of extremely high period that is unnoticeable in numerical calculations.

Fig. 7

Bifurcation diagram for the first case with bifurcations, showing the x coordinates of the attractors and other fixed points. Dashed lines represent saddle points or other unstable fixed points; solid lines represent asymptotically stable fixed points. The black dashed line on top of the gray solid line at $x = 0$ is included to show there are both stable and unstable fixed points with coordinate $x = 0$. As ${{v_2}}$ decreases, we see the birth of an invariant circle due to a Neimark–Sacker bifurcation. All x values on the invariant circle are mapped, such that for a single ${{v_2}}$ value there are many values of x. The invariant circle grows as ${v_2}$ decreases further until its sudden disappearance in a global bifurcation as the invariant circle runs into a heteroclinic cycle between the three saddle points of this case. Parameter values are ${R_1} = 1.2, {{R_2}} = 0.8, a = 0.3, {{b_1}} = 0.6, {{b_2}} = 0.2, {{v_1}} = 1.5$

Fig. 8

Bifurcation diagram of the second case showing the x coordinates of the attractors and other fixed points. Dashed lines show saddle points or other unstable fixed points; solid lines represent asymptotically stable fixed points. The black dashed line on top of the solid gray line at $x = 0$ shows there are both stable and unstable fixed points at $x = 0$. As with the first case, as ${{v_2}}$ decreases an invariant circle is born via a Neimark–Sacker bifurcation, eventually abruptly vanishing in a global bifurcation. Parameter values are $R_1 = 1.2, {R_2} = 0.8, a = 0.3, {b_1} = 0.6, {b_2} = 0.2, {v_1} = 2.85$ 

In a discrete-time system, it is not uncommon after a Neimark–Sacker bifurcation for the invariant circle to eventually deform or for further bifurcations to occur, leading to a strange attractor (Aronson et al. 1982; Curry and Yorke 1978; Neubert and Kot 1992; Selgrade and Namkoong 1991; Selgrade and Roberds 1996). Such behavior has been shown specifically in pioneer–climax discrete-time models (Selgrade and Namkoong 1991; Selgrade and Roberds 1996), as well as in a continuous-time pioneer–climax model (Selgrade 1994). We do not observe such behavior here. Instead, in both cases with a Neimark–Sacker bifurcation, a global bifurcation occurs as the invariant circle grows larger. We now turn to investigate these global bifurcations.

Global bifurcations

As can be seen in Figs. 5 and 6, in both cases with a Neimark–Sacker bifurcation the invariant circle grows closer and closer to the axes as ${{v_2}}$ decreases. Additionally, the invariant circle appears to approach some limiting curve with negative slope within the first quadrant. For a critical ${{v_2}}$ value, the invariant circle abruptly vanishes, leaving no attractor in the interior at all, in what is known as a dangerous global bifurcation (Thompson et al. 1994). The mechanism by which this disappearance of the invariant circle occurs is different in the two bifurcation cases, and we have two types of global bifurcation resulting in the blue-sky disappearance of an invariant circle.

Heteroclinic bifurcation

Fig. 9

Behavior of the system before, at, and after the heteroclinic bifurcation in the first case with bifurcations. Solid gray lines are stable manifolds, dashed gray lines are unstable manifolds. The axes are also invariant manifolds and form heteroclinic connections between ${{E_1}}$ and ${{E_2}}$ on the x-axis and ${{E_1}}$ and ${{E_3}}$ on the y-axis. Large black points are equilibria, trajectories of the system with an initial condition near ${E_5}$ are shown with smaller black points. The arrows indicate the direction of travel along the manifold. a, prior to the global bifurcation, the unstable manifold from ${{E_2}}$ passes below the stable manifold of ${{E_3}}$ and winds onto the attracting invariant circle. b, at the point of global bifurcation, a heteroclinic cycle forms and the unstable and stable manifolds of ${{E_2}}$ and ${{E_3}}$ link together to complete the cycle. c, after the global bifurcation, the heteroclinic connection is broken and no stable coexistence state exists. Almost all trajectories go to ${{E_4}}$ on the y-axis, save for solutions starting on the stable manifold of ${{E_3}}$ or on the axes between ${{E_1}}$ and ${{E_2}}$ or between ${{E_1}}$ and ${{E_3}}$. Parameter values are ${R_1} = 1.2, {{R_2}} = 0.8, a = 0.3, {{b_1}} = 0.6, {{b_2}} = 0.2, {{v_1}} = 1.5$

Fig. 10

Behavior of the system before, at, and after the homoclinic bifurcation in the second case with bifurcations. Solid gray lines are stable manifolds, dashed gray lines are unstable manifolds. The axes are again invariant manifolds with heteroclinic connections between ${{E_1}}$ and ${{E_2}}$ on the x-axis and ${{E_1}}$ and ${{E_3}}$ on the y-axis. Large black points represent equilibria, trajectories of the system starting near ${{E_5}}$ are shown with smaller black points. Arrows indicate direction of travel along the manifolds. a, prior to the global bifurcation, the unstable manifold from interior equilibrium ${{E_6}}$ passes below the stable manifold from ${{E_3}}$ on the y-axis and winds onto the attracting invariant circle. b, at the point of global bifurcation, a homoclinic connection forms between the stable and unstable manifolds of ${{E_6}}$, passing very close to the invariant manifolds on the axes. c, after the global bifurcation, no stable coexistence state exists and the homoclinic connection has been broken. Almost all trajectories go either to ${{E_4}}$ on the y-axis (shown here) or ${{E_2}}$ on the x-axis. Parameter values are ${R_1} = 1.2, {{R_2}} = 0.8, a = 0.3, {{b_1}} = 0.6, {{b_2}} = 0.2, {{v_1}} = 2.85$

Fig. 11

A magnified view of the lower part of the homoclinic connection in the second case with bifurcations. Here, the loop between the stable and unstable manifolds of ${{E_6}}$ is clear, and we more easily see that the invariant manifolds of ${E_6}$ do not intersect the manifolds on the axes of ${{E_1}}$ or ${{E_2}}$. Solid gray lines are stable manifolds, dashed gray lines are unstable manifolds, and points are equilibria. Arrows indicate direction of travel along the manifolds. Parameter values are ${R_1} = 1.2, {{R_2}} = 0.8, a = 0.3, {{b_1}} = 0.6, {{b_2}} = 0.2, {{v_1}} = 2.85,$ and ${{v_2}} = 2.878800266186866$

Recall that in the first case (Figs. 3a and 5), three saddle points are on the boundaries. These saddle points occur at ${{E_1}}$ at the origin, ${{E_2}}$ on the x-axis, and ${{E_3}}$, the lower fixed point on the y-axis. In this case, for a critical ${{v_2}}$ value, a heteroclinic cycle forms between the three saddle points. The invariant circle runs into the heteroclinic cycle and instantaneously vanishes. Thus, we have the blue-sky disappearance of an invariant circle due to a heteroclinic bifurcation.

The heteroclinic cycle that forms is relatively simple and is shown in Fig. 9b. The y-axis is the unstable manifold of ${{E_3}}$ and the stable manifold of ${{E_1}}$, so a heteroclinic connection is formed moving from ${{E_3}}$ into ${{E_1}}$ along the y-axis. Similarly, the x-axis is the unstable manifold of ${{E_1}}$ and the stable manifold of ${{E_2}}$; so the portion of the x-axis between ${{E_1}}$ and ${{E_2}}$ forms a heteroclinic connection as well. Finally, the unstable manifold of ${{E_2}}$ connects with the stable manifold of ${{E_3}}$, finishing off a triangular heteroclinic cycle between ${{E_1}}, {{E_2}}$ and ${{E_3}}$. Prior to and after the global bifurcation, the invariant manifolds of ${{E_2}}$ and ${{E_3}}$ do not intersect (Fig. 9a, c).

Numerically, this last heteroclinic connection between ${{E_2}}$ and ${{E_3}}$ appears to be a linear or near-linear connection (Fig. 9), suggesting these manifolds may coincide, or at the very least are nearly indistinguishable from one another. Kim and Marlin (1999) found in a differential-equation pioneer–climax model that, for the same geometric zero-growth isocline configuration, there was a critical parameter value for which the unstable manifold of ${{E_2}}$ and the stable manifold of ${{E_3}}$ coincided. While their continuous-time model was different from our discrete-time model, it offers support to the idea that these manifolds may indeed coincide; however, we leave an analytic proof for the future.

Homoclinic bifurcation

Three saddle points are also in the second case (Figs. 3e and 6), two on the boundaries at ${E_1}$ (the origin) and ${E_3}$ (the y-axis) and one just off of the x-axis at ${E_6}$. However, a heteroclinic cycle does not form between these three saddle points as in the first case; a homoclinic connection forms between the stable and unstable manifolds of the coexistence equilibrium ${E_6}$ instead (Figs. 10b and 11). The invariant circle runs into the homoclinic connection for a particular ${v_2}$ value, resulting in the blue sky disappearance of an invariant circle due to a homoclinic bifurcation.

The stable manifold of the interior equilibrium ${E_6}$ winds very close to the x and y axes but never intersects either axis, precluding a heteroclinic connection between ${E_6}$ and the origin (Fig. 11). Nevertheless, due to this closeness, we see the invariant circle grow very near to both axes. The invariant circle does not come close to the lower equilibrium on the y-axis ${E_3}$, however, even very close to the point at which the global bifurcation occurs. This distance can be seen by noting that the maximum height, or y coordinate, of the largest invariant circle is significantly lower in Fig. 6 than in Fig. 5. In addition to the fact that the manifolds of ${E_6}$ do not intersect the manifolds of ${E_1}$ or ${E_3}$, this behavior is further indicative of the differences between the heteroclinic cycle of the first case and the homoclinic connection in this second case.

Behavior near the global bifurcations

In both cases, as ${v_2}$ decreases the invariant circle grows larger until it is pushed up against the heteroclinic cycle or homoclinic connection (Figs. 5 and 6). At a critical ${v_2}$ value, the circle can grow no more, constrained by the invariant manifolds that have formed a heteroclinic cycle or homoclinic connection (Figs. 9b and 10b), and the invariant circle vanishes into thin air. The disappearance of the invariant circle corresponds to crossing the dashed-dotted line from above in Fig. 4. We see in Fig. 4 that for fixed ${v_1}$ an invariant circle exists only in the space between the solid line, when the invariant circle appears due to a Neimark–Sacker bifurcation, and the dashed-dotted line, which shows the point of global bifurcation when the invariant circle vanishes.

The disappearance of the invariant circle due to a global bifurcation is starkly apparent in a bifurcation diagram (Figs. 7 and 8). In both cases, we see the edges of the invariant circle, represented by the minimum and maximum x values of the invariant circle, run into the saddle point at ${E_2}$ or ${E_6}$, depending on the case, and the two saddle points ${E_1}$ and ${E_3}$ with coordinate $x = 0$.

The blue sky disappearance of an invariant circle due to a heteroclinic bifurcation is a global bifurcation not previously observed in a discrete-time pioneer–climax model, though Kim and Marlin (1999) noted this bifurcation for the first case in continuous time. This research is, however, the first time that anyone has observed the disappearance of an invariant circle due to a homoclinic bifurcation in a pioneer–climax system, whether in discrete or continuous time.

We have already noted that in both cases with bifurcations, ${E_5}$ was destabilized and an attracting quasiperiodic solution arose as a species’ competitive effect on the climax species grew stronger relative to its effect on the pioneer species. Now, we see that if this competitive effect on the climax species becomes stronger still, the attracting orbit vanishes entirely, leaving no stable coexistence state at all. This result is, perhaps, a ready conclusion to draw from observing the growth of the invariant circle. As ${v_2}$ decreases and competitive effect on the climax species becomes stronger, populations are regularly driven very close to either $x = 0$ or $y = 0$, such that a small stochastic effect on the populations may easily drive one species to extinction.

Discussion

As climate change affects ecosystems, causing effects from range shifts to increased habitat disturbances, successional communities are likely to become more common (Chen et al. 2011; Lenoir et al. 2008; Seidl et al. 2017). In this paper, we introduced a simple discrete-time pioneer–climax model where the species’ interaction dynamics are built upon principles of ecological succession. We characterized the long-term dynamics of the model in the first in-depth analytical stability analysis of a discrete-time pioneer–climax system, and the stability results were supported by numerical simulations. Our analysis is organized around nine geometric cases that are distinguished by different parameter regimes, corresponding to limits on the relative competition coefficients ${v_1}$ and ${v_2}$. The model results in a diverse set of long-term behaviors among the nine different cases, some of which are novel to a discrete-time pioneer–climax model. These behaviors were robust to other parameter sets.

Of these nine geometric configurations, six cases have an asymptotically stable climax-exclusion state in which the pioneer survives (Figs. 2a, c and 3b, c, e, f). For a long time, researchers assumed the asymptotic result of an undisturbed pioneer–climax system was pioneer exclusion (Harper 1977; Selgrade and Namkoong 1990). Many pioneer–climax models have since shown that pioneer exclusion is not the only possible long-term behavior, but this analysis was done primarily through proofs of the stabilities of coexistence states (e.g., Kim and Marlin (1999); Selgrade and Roberds (1996); Sumner (1996)). The existence of asymptotically stable climax-exclusion states was largely ignored in these analyses.

It is interesting, however, that in a majority of the nine parameter regimes the system does indeed have an asymptotically stable climax-exclusion state. This result tells us that even in an undisturbed environment, under some (perhaps limited) circumstances the pioneer species could not just coexist with the climax species, but it could even outcompete the climax species. The only cases in which a stable pioneer-only state was not present was for moderate values of ${v_1}$, i.e., when the effect of the pioneer on the climax species is relatively close to the effect of the pioneer species on itself.

Five of the six cases with an asymptotically stable climax exclusion state also have other asymptotically stable states (Figs. 2a and 3b, c, e, f). In addition to these five, the first case with bifurcations has multiple stable states for parameter values prior to the heteroclinic bifurcation (Fig. 3a). We therefore have six cases in total with multiple stable states that the system may tend to, depending on initial conditions and parameter values.

While the existence of multiple stable states is an ongoing debate in ecological theory, increasing evidence suggests that multiple stable states may indeed exist in nature (Beisner et al. 2003; Folke et al. 2004; Petraitis 2013; Scheffer and Carpenter 2003; Schröder et al. 2005). Being able to determine conditions under which an ecosystem may tend toward one stable state or another may be of particular interest in conservation and ecosystem management. Using this model, we provide a framework for understanding such conditions.

In all, the model presents four different cases, Fig. 3a, d, e, f, with stable coexistence states. In two of these cases, the interior fixed point is asymptotically stable for all relevant parameter values. In each of the other two cases, however, two bifurcations, one local and one global, change the nature of the coexistence state. For some parameter values, the coexistence state is an asymptotically stable fixed point. The Neimark–Sacker bifurcation that occurs in both cases destabilizes the interior fixed point and gives rise to an attracting quasi-periodic invariant circle, and a global bifurcation later causes the attractor to vanish entirely.

While many studies of pioneer–climax models have shown the existence and stability of the attractor that arises after a Hopf (or Neimark–Sacker) bifurcation (Selgrade and Namkoong 1990; Selgrade and Roberds 1996; Sumner 1996; Sumner 1998), only one considered the persistence of the attractor. One continuous-time system noted the abrupt disappearance of the attractor (Kim and Marlin 1999), but in only one geometric case. This disappearance corresponds to our first case with bifurcations and the associated heteroclinic bifurcation. In the second case with bifurcations, two coexistence states occur. We discovered that the formation of a homoclinic connection between the invariant manifolds of one coexistence equilibrium (${E_6}$) can result in the blue-sky disappearance of the attracting invariant circle that arose after the other coexistence equilibrium ${E_5}$ was destabilized. Thus, there are two cases for which the behavior of the system can change precipitously, causing a loss of any stable coexistence state and a sudden transition of the system to a new exclusion state. Knowledge of the conditions under which an attracting coexistence state vanishes may be important for ecosystem management.

While the structure of our model is simple, with rational fitness functions, we believe that this model captures the major dynamics of a pioneer–climax system and provides valuable insight into an inherently successional system. However, two species do not exist in a vacuum, and many more layers of interaction dynamics could be taken into account, as in Joshi and Blackmore's (2010) higher-dimensional model. To confirm the realism of this two-dimensional model, in the future it would be beneficial to compare our model predictions with data from pioneer–climax systems. Such data may be challenging to gather, given the long timescales involved in ecological succession. Additionally, the number of parameters in this model may make parameter estimation difficult.

A limitation of our model stems from a lack of survival from one generation to the next, despite many pioneer and climax plants living longer than a single generation. One simple way of accounting for survivorship is to add small positive linear terms, $s_1 x_t$ or $s_2 y_t$, to our model, where $s_1$ and $s_2$ represent the probability of surviving to the next generation. This addition would shift our per capita replacement rates (Fig. 1) upwards by a constant amount. Based on preliminary analytical and numerical exploration, the additional terms will not significantly change the behavior of our model, though some restrictions on the relative strength of survivorship may be necessary to ensure stability.

Spatiotemporal dynamics of successional communities are increasingly relevant in both conservation and invasion ecology. Through spatiotemporal modeling of a pioneer-climax system, we can approach problems like determining the critical patch size for population persistence or quantifying the rate of invasion into a new habitat by calculation of a species’ spreading speed. Additionally, moving habitat models can inform us whether or not populations can keep pace with climate-induced range shifts. With such germane and interesting problems, future work should focus on incorporating spatial components into our model.

After the many initial studies into the basic dynamics of pioneer–climax models and the potential for reversing bifurcations and restoring stability in those models (e.g., Kim and Marlin 1999; Selgrade and Namkoong 1990; Selgrade and Roberds 1997; Sumner 1996), research on pioneer–climax systems has indeed made the shift to spatiotemporal models, but it has largely focused on reaction-diffusion models. Much of this research has focused on proving the existence of traveling waves (Brown et al. 2005; Buonomo and Rionero 2009; Cao and Weng 2017; Huang and Weng 2013; Yu et al. 2013; Yuan and Zou 2010). Other studies have looked at Turing and Hopf bifurcations in reaction-diffusion systems (Buchanan 2005; Liu and Wei 2011). Two isolated studies have examined spreading speed within a reaction-diffusion context (Cao and Weng 2017; Weng and Zou 2014).

These studies provide a valuable starting point in the analysis of spatiotemporal dynamics of pioneer–climax systems. However, the problems mentioned above on critical patch size and moving habitat models are wholly unexamined, and only two studies have examined spreading speed (Cao and Weng 2017; Weng and Zou 2014). Furthermore, reaction-diffusion models are limited in their ecological realism, as dispersal is often non-diffusive and growth, particularly in plants, is frequently seasonal (Lewis et al. 2016; Lutscher 2019).

Using the discrete-time model presented here, however, we can develop a more ecologically realistic spatiotemporal integrodifference equation model to examine the same problems. These models have the benefit of separating growth and dispersal, and they take long-distance dispersal into account in a way reaction-diffusion models do not. With flexible dispersal kernels and separation of growth and dispersal phases, a common feature in plant populations, integrodifference equations provide a more ecologically realistic method of examining some of these spatial issues in invasions and conservation (Kot and Schaffer 1986; Lewis et al. 2016; Lutscher 2019).

Before moving to these more complex questions, it is critical that we understand the dynamics of the underlying foundational models, like the one analyzed here. Our work in analyzing this discrete-time system not only advances our knowledge of stationary pioneer–climax dynamics, but also paves the way for situating our model within a spatial context. With this framework, we have the ability to examine how these species will persist, invade new habitat, and survive in an ever-changing environment.

Appendices

Appendix

Preliminary organization and derivatives for the Jacobian

Throughout all appendices, we will refer to the fixed points as ${E_1}$–${E_6}$, and their individual components as $x_1$–${x_6}$ and $y_1$–${y_6}$. Equations 20–25 provide the explicit forms of the equilibria, which are used only in rare cases. We now drop all time-dependent notation, as we are analyzing fixed points $\left( x, y\right)$ that, by definition, do not change over time. In all of our calculations, we consider only nonnegative x and y, and therefore nonnegative m and n. All parameters are positive.

Table 2 Parameter restrictions for each of the nine geometric cases. For the six parameter regimes with at least one interior fixed point, these conditions are the same as those for existence of the interior fixed point(s) in the first quadrant. A dash in the “Case” column indicates that no specific label has been assigned to that case. A dash in the “Interior fixed point” column indicates no interior fixed point exists for that parameter regime. There is no Condition 3 for the three parameter regimes with no interior equilibria, as that condition was derived by restricting an interior fixed point to be within the first quadrant

We first recall that there are nine different geometric cases (illustrated in Figs. 2 and 3) for this system, and up to six equilibria that exist in each case. The parameter restrictions bounding each of the nine cases are summarized in Table 2. In three cases, only boundary fixed points exist. In the other six cases, one or both interior fixed points occur in addition to the boundary equilibria. Stability analyses for the boundary fixed points are relatively straightforward, and each boundary fixed point can be analyzed individually and its stability classified for each case based on the parameter restrictions from Table 2. The interior fixed points present more of a challenge, and are examined on a case-by-case basis. The inequalities presented in Table 2 are also heavily used in classifying the stabilities for the interior fixed points. The stabilities of the fixed points for all nine cases are summarized in Table 3.

Table 3 Stability of each equilibrium point in the nine different cases. A dash indicates the fixed point does not exist in the first quadrant for that parameter regime

To begin our stability analyses, we first present some general results that we will use extensively. After this initial setup is complete, we proceed to the boundary fixed points and classify the stability of each boundary fixed point for all nine cases by examining the magnitudes of the eigenvalues of the Jacobian (Eq. 26). We finally move to the interior fixed-points and lay out a few general results that are specific to the interior fixed points and the eigenvalues of the Jacobian (Eq. 26) evaluated at the interior fixed points. Our stability analyses for the interior fixed points proceed on a case-by-case basis, and for each case, we use the Jury conditions (Eqs. 27, 28 and 29) and Descartes’ Rule of Signs (Descartes 1886) to bound the eigenvalues of the Jacobian evaluated at the given interior fixed point.

Derivatives for the Jacobian

Recall that a fixed point is asymptotically stable if both eigenvalues of the Jacobian (Eq. 26) evaluated at the fixed point are less than 1 in magnitude. If stability of a fixed point is classified through use of the Jury conditions, the trace $\tau$ and determinant $\Delta$ of the Jacobian are important quantities for consideration. When evaluated at a fixed point, our Jacobian matrix (Eq. 26) has elements containing the derivatives of the fitness functions, including the particular quantities $1 + x_i \ f'\left( m_i\right)$ and $1 + y_i \ g'\left( n_i\right)$, where $x_i$ and $y_i$ are the coordinates of the fixed point and $m_i$ and $n_i$ are also evaluated at that fixed point. Before beginning our stability analyses, it is useful to present some general results pertaining to the derivatives of the fitness functions and these key quantities.

First, we explicitly calculate the forms of the derivatives $f'\left( m\right)$ and $g'\left( n\right)$, where

$$\begin{aligned} f'\left( m\right) = \frac{-a R_1}{\left( 1 + a m \right) ^2} < 0 \end{aligned}$$

(31)

for all m. We will write $g'\left( n\right)$ in a way that is particularly useful by first considering the expression

$$\begin{aligned} 1 - g\left( n\right) =&{} \ \frac{1 + {b_2} n^2 - {R_2} \left( 1 + {b_1} n\right) }{1 + {b_2} n^2} \\ =&{} \ \frac{{b_2} \left( n - {n_1}^*\right) \left( n - {n_2}^*\right) }{1 + {b_2} n^2}, \end{aligned}$$

(32)

where we factored the numerator to obtain the second equality. Differentiating both sides with respect to n, we find

$$\begin{aligned} g'\left( n\right) =&{} \ \frac{-{b_2} \left[ \left( n - {n_1}^*\right) + \left( n - {n_2}^*\right) \right] }{1 + {b_2} n^2} \\ &{} \ + \ \frac{2{b_2}^2 n \left( n - {n_1}^*\right) \left( n - {n_2}^*\right) }{\left( 1 + {b_2} n^2\right) ^2}, \end{aligned}$$

(33)

which may be either positive or negative. In particular, note that for $n = {n_1}^*$, $g'\left( n\right) > 0$, and when $n = {n_2}^*$, $g'\left( n\right) < 0$.

Next, note that

$$\begin{aligned} \frac{\partial }{\partial x} \left[ f(m) x \right] = f\left( m\right) + x \ f'\left( m\right) \frac{\partial m}{\partial x} = f\left( m\right) + x \ f'\left( m\right) , \end{aligned}$$

(34)

which is the upper-left entry in our Jacobian matrix (Eq. 26) for this system. Explicitly calculating this quantity, we find

$$\begin{aligned} \frac{\partial }{\partial x} \left[ f(m) x \right] = \frac{R_1 \left( 1 + a {v_2} y\right) }{\left[ 1 + a\left( x + {v_2} y\right) \right] ^2} > 0. \end{aligned}$$

(35)

Thus, for all nonnegative x and y,

$$\begin{aligned} f\left( m\right) + x \ f'\left( m\right) > 0. \end{aligned}$$

(36)

At fixed points ${E_2}$, ${E_5}$, and ${E_6}$, $f\left( m\right) = 1$, and so Eq. 36 simplifies even further to

$$\begin{aligned} 1 + x_i \ f'\left( m_i\right) > 0 \end{aligned}$$

(37)

for $i = 2, 5, 6$, $m_i = x_i + {v_2} y_i$.

We now consider the lower-right entry of the Jacobian matrix. As above,

$$\begin{aligned} \frac{\partial }{\partial y} \left[ g(n) y \right] = g\left( n\right) + y g'\left( n\right) \frac{\partial n}{\partial y} = g\left( n\right) + y g'\left( n\right) . \end{aligned}$$

(38)

Unlike Eq. 36, we do not prove this quantity is positive in general. However, it is positive at each of the four equilibria ${E_3}$–${E_6}$ for which $y \ne 0$ and $g\left( n\right) = 1$. We now prove this statement for each of the four fixed points individually.

At ${E_3}$, recall that $y = n = {n_1}^*$. Substituting this relation into Eq. 38 and using Eq. 33 to simplify the derivative, we obtain

$$\begin{aligned} g\left( n\right) + y g'\left( n\right) = 1 + \frac{{b_2} {n_1}^* \left( {n_2}^* - {n_1}^* \right) }{1 + {b_2} \left( {n_1}^* \right) ^2} > 0, \end{aligned}$$

(39)

as ${n_2}^* > {n_1}^*$.

Similarly, at ${E_4}$, $y = n = {n_2}^*$. Again substituting this relation into Eq. 38, we find

$$\begin{aligned} g\left( n\right) + y g'\left( n\right) = 1 - \frac{{b_2} {n_2}^* \left( {n_2}^* - {n_1}^* \right) }{1 + {b_2} \left( {n_2}^* \right) ^2} = \frac{1 + {b_2} {n_1}^* {n_2}^*}{1 + {b_2} \left( {n_2}^* \right) ^2} > 0, \end{aligned}$$

(40)

after finding a common denominator and simplifying.

For ${E_5}$, we have $y = {y_5}$ and $n = {v_1}{x_5} + {y_5} = {n_1}^*$. Using the same process as above,

$$\begin{aligned} g\left( n\right) + y g'\left( n\right) = 1 + \frac{{b_2} {y_5} \left( {n_2}^* - {n_1}^* \right) }{1 + {b_2} \left( {n_1}^* \right) ^2} > 0. \end{aligned}$$

(41)

Finally, at ${E_6}$, $y = {y_6}$ and $n = {v_1} {x_6} + {y_6} = {n_2}^*$. Instead of using the fact that $g\left( {n_2}^*\right) = 1$, we substitute the explicit form of $g\left( {n_2}^*\right)$, which has the same denominator as that of $g'\left( {n_2}^*\right)$. We find that

$$\begin{aligned} g\left( n\right) + y g'\left( n\right) = \frac{{R_2}(1 + {b_1} {n_2}^*) - {b_2} {y_6} \left( {n_2}^* - {n_1}^* \right) }{1 + {b_2} \left( {n_2}^* \right) ^2}. \end{aligned}$$

(42)

The denominator is positive, so if the numerator is positive then Eq. 42 is positive as desired. Considering the numerator of Eq. 42, we begin by substituting for ${n_2}^* - {n_1}^*$, yielding

$$\begin{aligned} &{}{R_2}\left( 1 + {b_1} {n_2}^*\right) - {b_2} {y_6} \left( {n_2}^* - {n_1}^* \right) \\ &{} = {R_2}\left( 1 + {b_1} {n_2}^*\right) - {y_6} \sqrt{{R_2}^2 {b_1}^2 - 4 {b_2} \left( 1 - {R_2} \right) }. \end{aligned}$$

(43)

For convenience, denote

$$\begin{aligned} c = \sqrt{{R_2}^2 {b_1}^2 - 4 {b_2} \left( 1 - {R_2} \right) } \end{aligned}$$

(44)

and substitute ${n_2}^* = {v_1} {x_6} + {y_6}$ into Eq. 43 to obtain

$$\begin{aligned} &{}{R_2}\left[ 1 + {b_1} \left( {v_1} {x_6} + {y_6} \right) \right] - c {y_6} \\ &{} = {R_2}\left( 1 + {b_1} {v_1} {x_6}\right) + {y_6} \left( {R_2} {b_1} - c\right) . \end{aligned}$$

(45)

Finally, using

$$\begin{aligned} &{} {n_1}^* = \frac{{R_2} {b_1} - c}{2 {b_2}} \\ &{} \implies 2{b_2} {n_1}^* = {R_2}{b_1} - c, \end{aligned}$$

(46)

we can rewrite Eq. 45 as

$$\begin{aligned} {R_2}\left( 1 + {b_1} {v_1} {x_6}\right) + {y_6} \left( 2 {b_2} {n_1}^*\right) , \end{aligned}$$

(47)

which is clearly positive. Thus, the numerator of Eq. 42 is positive and therefore Eq. 42 is positive.

We have just shown that $g\left( n\right) + y g'\left( n\right) > 0$ when evaluated at each of the fixed points ${E_3}$–${E_6}$, i.e.,

$$\begin{aligned} 1 + y_i \ g'\left( n_i\right) > 0 \end{aligned}$$

(48)

for $i = 3, 4, 5, 6$, $n_i = {v_1} x_i + y_i$.

Stability of boundary fixed points

With these general quantities and inequalities laid out, we now establish the stability of the four boundary fixed points for all nine parameter regimes. Refer to Table 1 for the eigenvalues for each fixed point ${E_1}$–${E_4}$.

First, note that all eigenvalues are positive and real for the boundary fixed points. It is clear that the eigenvalues for ${E_1}$ are positive and real, as the eigenvalues are the populations’ net reproductive rates, which must be positive. Furthermore, as the fitness functions $f\left( m\right)$ and $g\left( n\right)$ are always real and greater than zero, it is evident that ${\lambda _1} > 0$ and is real for ${E_2}$–${E_4}$. Next, consider ${\lambda _2}$ as given in Table 1. As demonstrated in Appendix A (Preliminary organization and derivatives for the Jacobian) by Eqs. 37, 39, and 40, we see that ${\lambda _2} > 0$ for ${E_2}$–${E_4}$. Thus, both eigenvalues are positive for all four boundary fixed points. To prove stability of ${E_1}$–${E_4}$, it is therefore sufficient to show whether each eigenvalue ${\lambda _1}$ and ${\lambda _2}$ is greater or less than $+1$.

With the knowledge that $f'\left( m\right) < 0$ always, $g'\left( {n_1}^*\right) > 0$, $g'\left( {n_2}^*\right) < 0$, and the results shown in Appendix A (Preliminary organization and derivatives for the Jacobian), we may continue with proving the stability of the boundary fixed-points for the nine cases shown in Figs. 2 and 3.

Stability of trivial extinction state ${E_1}$

As mentioned in Subsection 4.1 (Boundary equilibria), the equilibrium ${E_1}$ at the origin is always a saddle point, as the associated eigenvalues of the Jacobian are $R_1 > 1$ and $0< {R_2} < 1$.

Stability of climax-exclusion state ${E_2}$

By construction of the climax fitness function $g\left( n\right)$, there are two roots n, ${n_1}^*$ and ${n_2}^*$, where $g\left( n\right) = 1$. For values of n below and above those roots, $g\left( n\right) < 1$, and for values of n between the roots, $g\left( n\right) > 1$ (see Fig. 1b).

Consider $n = {v_1} x$. Then, $g\left( n\right) = g\left( {v_1} x \right)$. The two roots where $g\left( {v_1} x \right) = 1$ are given by

$$\begin{aligned} {v_1} {x_\alpha} = {n_1}^* \implies {x_\alpha} = {n_1}^*/{v_1},\end{aligned}$$

(49)

$$\begin{aligned} {v_1} {x_\beta} = {n_2}^* \implies {x_\beta} = {n_2}^*/{v_1}. \end{aligned}$$

(50)

Thus, the value of ${\lambda _1} = g\left( {v_1} x_2 \right)$ depends on where $x_2$ falls in relation to ${x_\alpha}$ and ${x_\beta}$ from Eqs. 49 and 50. Let ${n_2} = {v_1} x_2$. We see that

$$\begin{aligned} x_2< {x_\alpha} \implies {n_2}< {n_1}^*&\implies {\lambda _1}= g\left( {n_2} \right) < 1, \end{aligned}$$

(51)

$$\begin{aligned} x_2> {x_\beta} \implies {n_2} > {n_2}^*&\implies {\lambda _1} = g\left( {n_2} \right) < 1, \end{aligned}$$

(52)

$$\begin{aligned} {x_\alpha}< x_2< {x_\beta} \implies {n_1}^*< {n_2} < {n_2}^*&\implies {\lambda _1} = g\left( {n_2} \right) > 1. \end{aligned}$$

(53)

Observe that ${x_\alpha}$ and ${x_\beta}$ are the x-intercepts of the two zero-growth isoclines where $g\left( n\right) = 1$, while $x_2$ is the x-intercept of the zero-growth isocline $f\left( m\right) = 1$. Thus, we can read off the magnitude of ${\lambda _1}$ straight from the configuration of the zero-growth isoclines for a particular case.

Recall that ${\lambda _2} = f\left( x_2\right) + x_2 \ f'\left( x_2\right) = 1 + x_2 f'\left( x_2\right)$. By Eq. 37, ${\lambda _2} > 0$. As $f'\left( m\right) < 0$ for all m, ${\lambda _2} < 1$.

Combining our analysis for both eigenvalues, we see that stability depends on ${\lambda _1}$, as $0< {\lambda _2} < 1$. ${E_2}$ will be a stable node if the nontrivial pioneer zero-growth isocline intercepts the x-axis to the left or right of both nontrivial climax zero-growth isoclines ($x_2 < {x_\alpha}$ or $x_2 > {x_\beta}$). It will be a saddle point if the nontrivial pioneer zero-growth isocline intercepts the x-axis between the two nontrivial climax zero-growth isoclines (${x_\alpha}< x_2 < {x_\beta}$). See Figs. 2 and 3 for details. The stability of ${E_2}$ for each of the nine cases is summarized in Table 3.

Stability of pioneer-exclusion states ${E_3}$ and ${E_4}$

Now, consider the construction of the pioneer fitness function. There is one root m where $f\left( m\right) = 1$. For values of m below this root, $f\left( m\right) > 1$, and for values of m above this root, $f\left( m\right) < 1$ (see Fig. 1a).

Take $m = {v_2} y$. Then $f\left( m\right) = f\left( {v_2} y\right)$. The root where $f\left( {v_2} y \right) = 1$ is

$$\begin{aligned} {v_2} {y_\alpha} = \frac{R_1 - 1}{a} = {m^*} \implies {y_\alpha} = {m^*}/{v_2}. \end{aligned}$$

(54)

As with ${E_2}$, for ${E_3}$ the value of ${\lambda _1} = f\left( {v_2} {y_3}\right)$ depends on the size of ${y_3}$ relative to ${y_\alpha}$. Similarly, for ${E_4}$ the value of ${\lambda _1} = f\left( {v_2} {y_4}\right)$ depends on the size of ${y_4}$ relative to ${y_\alpha}$. Let ${m_3} = {v_2} {y_3}$ and ${m_4} = {v_2} {y_4}$.

For ${E_3}$,

$$\begin{aligned} {y_3}< {y_\alpha} \implies {m_3} < {m^*} \implies {\lambda _1} = f\left( {m_3}\right) > 1,\end{aligned}$$

(55)

$$\begin{aligned} {y_3}> {y_\alpha} \implies {m_3} > {m^*} \implies {\lambda _1} = f\left( {m_3}\right) < 1. \end{aligned}$$

(56)

For ${E_4}$,

$$\begin{aligned} {y_4}< {y_\alpha} \implies {m_4} < {m^*} \implies {\lambda _1} = f\left( {m_4}\right) > 1,\end{aligned}$$

(57)

$$\begin{aligned} {y_4}> {y_\alpha} \implies {m_4} > {m^*} \implies {\lambda _1} = f\left( {m_4}\right) < 1. \end{aligned}$$

(58)

As in our analysis for ${E_2}$, note that ${y_\alpha}$ is the y-intercept of the zero-growth isocline $f\left( m\right) = 1$, while ${y_3}$ and ${y_4}$ are the y-intercepts of the two zero-growth isoclines where $g\left( n\right) = 1$. We can again obtain the magnitude of ${\lambda _1}$ for both ${E_3}$ and ${E_4}$ directly from the configuration of the zero-growth isoclines for a particular case.

Turning to ${\lambda _2} = 1 + y_i \ g'\left( y_i\right)$, $i = 3$ or 4, we have shown that ${\lambda _2} > 0$ for ${E_3}$ (Eq. 39) and ${E_4}$ (Eq. 40). At ${E_3}$, $y = {n_1}^*$ and $g'\left( {n_1}^*\right) > 0$. Thus, ${\lambda _2} > 1$ at ${E_3}$. Conversely, at ${E_4}$, $y = {n_2}^*$ and $g'\left( {n_2}^*\right) < 0$, so ${\lambda _2} < 1$ at ${E_4}$.

The natures of the equilibrium points ${E_3}$ and ${E_4}$ depend on their respective ${\lambda _1}$. For ${E_3}$, ${\lambda _2} > 1$, so ${E_3}$ is an unstable node if the lower nontrivial climax zero-growth isocline intercepts the y-axis below the nontrivial pioneer zero-growth isocline (${y_3} < {y_\alpha}$), and is a saddle point if the lower nontrivial climax zero-growth isocline intercepts the y-axis above the nontrivial pioneer zero-growth isocline (${y_3} > {y_\alpha}$).

For ${E_4}$, $0< {\lambda _2} < 1$, so ${E_4}$ is a saddle point if the upper nontrivial climax zero-growth isocline intercepts the y-axis below the nontrivial pioneer zero-growth isocline (${y_4} < {y_\alpha}$), and is a stable node if the upper nontrivial climax zero-growth isocline intercepts the y-axis above the nontrivial pioneer zero-growth isocline (${y_4} > {y_\alpha}$). See Figs. 2 and 3 for illustration. The stability of ${E_3}$ and ${E_4}$ for each of the nine cases is summarized in Table 3.

Stability of interior fixed points

We now turn to the stability of the two interior fixed points ${E_5}$ and ${E_6}$, given by Eqs. 24 and 25. Rather than examining each fixed point separately and handling all cases for a particular fixed point at once, as in Appendix B (Stability of boundary fixed points), it is simpler to organize this section by individual case. The stability analyses proceed by checking the three Jury conditions (Eqs. 27, 28 and 29) for each of the six cases in which one or both interior fixed points occur. We also frequently use Descartes’ Rule of Signs (Descartes 1886) to bound the eigenvalues of the Jacobian. For convenience, we will refer to these six cases by the letter of their corresponding subfigure in Fig. 3, e.g., case a corresponds to Fig. 3a. Before beginning on a case-by-case basis, we present several results that simplify the overall analysis.

Conditions for existence of interior fixed points

For each of the six cases a–f, three conditions on ${v_1}$ and ${v_2}$ must be satisfied to guarantee existence of an interior fixed point in the first quadrant. The conditions are easily derived by ensuring that both numerator and denominator are of the same sign for the pairs $\left( {x_5}, {y_5}\right)$ and $\left( {x_6}, {y_6}\right)$ – one condition comes from setting the numerator of the x-coordinate to be positive (negative), one condition from setting the numerator of the y-coordinate to be positive (negative), and one condition from setting the shared denominator $1 - {v_1} {v_2}$ to be positive (negative). These conditions are presented in Table 2 and will be heavily used in proving stability of the interior fixed points. From now on, when we say that an interior fixed point exists, we are specifically referring to an interior fixed point that exists inside the first quadrant, i.e., $x > 0$ and $y > 0$.

Descartes’ rule of signs

In addition to the conditions from Table 2, we make liberal use of Descartes’ Rule of Signs (Descartes 1886) throughout our analyses, and provide a brief refresher on this useful technique here. Consider the characteristic polynomial for the Jacobian given by

$$\begin{aligned} \upchi \left( \lambda \right) = \lambda ^2 - \tau \lambda + \Delta , \end{aligned}$$

(59)

whose roots are the eigenvalues of our Jacobian (Eq. 26) and where the coefficients $\tau$ and $\Delta$ are the trace and determinant of the Jacobian. According to Descartes’ Rule of Signs, the number of positive real roots of Eq. 59— and therefore positive real eigenvalues of the Jacobian — is equal to the number of sign changes between consecutive coefficients of the polynomial, or less than it by some even number (Descartes 1886). Thus, if there are two sign changes, we have either two or zero positive roots.

In addition, we can find the number of negative real roots of the characteristic polynomial by counting the number of sign changes in the coefficients of

$$\begin{aligned} \upchi \left( -\lambda \right) = \lambda ^2 + \tau \lambda + \Delta , \end{aligned}$$

(60)

Again, the number of negative roots is equivalent to the number of sign changes or less by an even number, so if there are two sign changes in Eq. 60, there are two or zero negative roots.

By introducing a change of variable, we can use Descartes’ Rule of Signs once more to check how many roots are greater than $+1$, simplifying the stability analysis for several cases. To analyze eigenvalues around $+1$, we introduce the shift

$$\begin{aligned} \mu = \lambda - 1 \implies \lambda = 1 + \mu . \end{aligned}$$

(61)

Substituting this into Eq. 59, we obtain

$$\begin{aligned} \upchi \left( \mu \right) = \mu ^2 + \left( 2 - \tau \right) \mu + \left( 1 - \tau + \Delta \right) , \end{aligned}$$

(62)

and note that the constant term is the same as the left-hand side of the first Jury condition (Eq. 27). Thus, any roots $\mu > 0$ correspond to $\lambda > 1$.

General results for $\tau$ and $\Delta$

Finally, we present a few useful quantities and relationships pertaining to the interior fixed points. These are summarized in Table 4. In particular, we examine the trace $\tau$ and determinant $\Delta$ of our Jacobian (Eq. 26) and use the results of our examination to simplify the first and second Jury conditions as well as to bound the eigenvalues of the Jacobian for both interior fixed points.

Table 4 Quantities and relationships frequently used in Descartes’ Rule of Signs or in checking the Jury conditions

Recall that for both interior fixed points ${E_5}$ and ${E_6}$, $f\left( m\right) = 1$, $g\left( n\right) = 1$, and $m = {m^*}$ after substituting in the explicit forms of the equilibria. For ${E_5}$, $n = {n_1}^*$ and for ${E_6}$, $n = {n_2}^*$ after substitution of the equilibria.

We can write the trace of the Jacobian, $\tau$, as

$$\begin{aligned} \tau = 2 + x \ f'\left( m\right) + y \ g'\left( n\right) . \end{aligned}$$

(63)

Note that $\tau$ is positive for the Jacobian evaluated at both ${E_5}$ and ${E_6}$, as we have shown that both $1 + x f'\left( m\right) > 0$ and $1 + y g'\left( n\right) > 0$ for ${E_5}$ and ${E_6}$ (Eqs. 37, 41, 42).

The determinant of the Jacobian, $\Delta$, can be written as

$$\begin{aligned} \Delta = \left[ 1 + x \ f'\left( m\right) \right] \left[ 1 + y \ g'\left( n\right) \right] - {v_1} \ {v_2} \ x\ y \ f'\left( m\right) g'\left( n\right) \end{aligned}$$

(64)

for $\left( x, y\right) = \left( {x_5}, {y_5}\right)$ or $\left( {x_6}, {y_6}\right)$.

It is clear that $\Delta > 0$ when the Jacobian is evaluated at ${E_5}$. By the results of Eqs. 37 and 41, the first two terms of Eq. 64 are positive, and therefore their product is positive. At ${E_5}$, $g'\left( n\right) > 0$, and thus the entire third term of Eq. 64, $\left[ - {v_1} \ {v_2} \ x\ y \ f'\left( m\right) g'\left( n\right) \right]$, is positive. Therefore, $\Delta > 0$ at ${E_5}$.

It is not as straightforward to show that $\Delta > 0$ when the Jacobian is evaluated at ${E_6}$, but this statement is indeed true for all four cases c–f in which ${E_6}$ exists. We rewrite the determinant from Eq. 64, plugging in the explicit fixed point coordinates given in Eq. 25 and simplified derivatives from Eqs. 31 and 33. After expanding and finding a common denominator, then rearranging and grouping like terms, we obtain the very long expression

$$\begin{aligned} \Delta =& \ \frac{\left( {v_1}{v_2} - 1\right) + {b_2} {n_1}^*\left( {m^*} {v_1} - {n_2}^*\right) + {b_2} {v_1} {n_2}^* \left( {v_2} {n_2}^* - {m^*}\right) }{\left( 1 + a {m^*}\right) \left[ 1 + {b_2}\left( {n_2}^*\right) ^2\right] \left( {v_1} {v_2} - 1\right) } \\ &{}+ \frac{a {v_2} \left( {m^*} {v_1} - {n_2}^*\right) + a {b_2} {n_1}^* {n_2}^* {v_2} \left( {v_1} {m^*} - {n_2}^*\right) }{\left( 1 + a {m^*}\right) \left[ 1 + {b_2}\left( {n_2}^*\right) ^2\right] \left( {v_1} {v_2} - 1\right) }. \end{aligned}$$

(65)

In cases c and e, we see from condition three in Table 2 that the denominator of Eq. 65 is positive. Thus, we require the numerator be positive in these two cases for $\Delta > 0$. First, note that all quantities multiplying each parenthetical term are positive (this statement is true for all four cases c–f). Each parenthetical quantity in the numerator is in fact a condition from Table 2, and so we know the signs of these quantities. From the first condition for cases c and e, ${m^*} {v_1} - {n_2}^* > 0$, from the second, ${v_2} {n_2}^* - {m^*} > 0$, and we have already used the third condition. Thus, the numerator is also positive, and $\Delta > 0$ for the Jacobian evaluated at ${E_6}$ in cases c and e.

For cases d and f, the third condition from Table 2 shows that the denominator of Eq. 65 is negative, so we would like to show that the numerator is negative for $\Delta > 0$. As with cases c and e, each parenthetical quantity corresponds to a condition from Table 2, and we now know ${m^*} {v_1} - {n_2}^* < 0$ and ${v_2} {n_2}^* - {m^*} < 0$, in addition to the third condition we have already used. Thus, the numerator is a sum of negative quantities, and is therefore negative. Once again, $\Delta > 0$ for the Jacobian evaluated at ${E_6}$ in cases d and f.

Given that $\tau > 0$ and $\Delta > 0$ for the Jacobian evaluated at both ${E_5}$ and ${E_6}$, the second Jury condition $1 + \tau + \Delta > 0$ will clearly always be satisfied for both interior fixed points for all cases a–f (assuming the interior fixed point exists).

Additionally, observe that $\tau > 0$ and $\Delta > 0$ means there are two sign changes between coefficients in Eq. 59 and zero sign changes in Eq. 60. Therefore, by Descartes’ Rule of Signs, there are zero negative eigenvalues and either two or zero positive eigenvalues of the Jacobian. Thus, for both interior fixed points in all cases a through f (assuming the interior fixed point exists), all real eigenvalues are positive. Furthermore, we either have two real, positive eigenvalues, or a pair of complex-conjugate eigenvalues for both interior fixed points in all cases.

Finally, we simplify the expression for the first Jury condition using our knowledge of $\tau$ and $\Delta$. After substituting Eqs. 63 and 64 into the left-hand side of the first Jury condition (Eq. 27), we obtain

$$\begin{aligned} 1 - \tau + \Delta = x \ y \ f'\left( m\right) g'\left( n\right) \left( 1 - {v_1} {v_2} \right) \end{aligned}$$

(66)

for $\left( x, y\right) = \left( {x_5}, {y_5}\right)$ or $\left( {x_6}, {y_6}\right)$. As already stated, we consider only positive equilibria (see Subsection 4.2, Interior equilibria, for a discussion of the border cases where an interior fixed point hits an axis as a parameter is varied). We know that $f'\left( m\right) < 0$ always. Thus, ascertaining whether the first Jury condition is satisfied comes down to the sign of $g'\left( n\right)$, which is positive at ${E_5}$ and negative at ${E_6}$, and the sign of $1 - {v_1} {v_2}$, which is given by condition three in Table 2. We now proceed to a case-by-case stability analysis, recalling that a summary of the general results presented in this section are found in Table 4. This information will be particularly useful when examining cases e and f.

Case A

In case a, we have one interior fixed point ${E_5}$. The first Jury condition

$$\begin{aligned} 1 - \tau + \Delta = {x_5} \ {y_5} \ f'\left( {m_5}\right) g'\left( {n_5}\right) \left( 1 - {v_1} {v_2} \right) > 0 \end{aligned}$$

(67)

by condition three in Table 2. Thus, the first Jury condition is always satisfied. We already know that the second Jury condition is always satisfied.

The third Jury condition is $\Delta < 1$. After expanding the expression for $\Delta$ from Eq. 64 we see that the third Jury condition $\Delta < 1$ is equivalent to

$$\begin{aligned} x \ f'\left( m\right) + y \ g'\left( n\right) + x\ y \ f'\left( m\right) g'\left( n\right) \left( 1 - {v_1} \ {v_2} \right) < 0. \end{aligned}$$

(68)

This expression does not simplify further in a useful manner. Instead, we look to express the left-hand side of Eq. 68 in terms of our relative competition coefficients ${v_1}$ and ${v_2}$, and find the curve in ${v_1}$–${v_2}$ space where $\Delta = 1$ (or Eq. 68 is equal to zero). In other words, we want to find the boundary where we violate the third Jury condition, such that crossing that boundary causes a change of stability of ${E_5}$ in case a.

To find this boundary in ${v_1}$–${v_2}$ space, we will substitute the explicit forms of ${x_5}$ and ${y_5}$ into the left-hand side of Eq. 68. Recall from Table 4 that ${m_5} = {m^*}$ and ${n_5} = {n_1}^*$. As the derivatives $f'\left( {m^*}\right)$ and $g'\left( {n_1}^*\right)$ do not have ${v_1}$ or ${v_2}$ dependence, we can leave them in their general form. After this substitution and minor simplification, the boundary where we violate the third Jury condition is given by

$$\begin{aligned} 0 = & \ \frac{\left( {m^*} - {v_2} {n_1}^*\right) f'\left( {m^*}\right) }{1 - {v_1} {v_2}} + \frac{\left( {n_1}^* - {v_1} {m^*} \right) g'\left( {n_1}^*\right) }{1 - {v_1} {v_2}} \\ &{}+ \frac{\left( {m^*} - {v_2} {n_1}^*\right) \left( {n_1}^* - {v_1} {m^*}\right) f'\left( {m^*}\right) g'\left( {n_1}^*\right) }{1 - {v_1} {v_2}}. \end{aligned}$$

(69)

From condition three in Table 2, we know $1 - {v_1} {v_2} \ne 0$ while in this case, and so we safely multiply both sides by $1 - {v_1} {v_2}$ to obtain

$$\begin{aligned} 0 \ = &\;\left( {m^*} - {v_2} {n_1}^*\right) f'\left( {m^*}\right) + \left( {n_1}^* - {v_1} {m^*} \right) g'\left( {n_1}^*\right) \\ &{}+ \left( {m^*} - {v_2} {n_1}^*\right) \left( {n_1}^* - {v_1} {m^*}\right) f'\left( {m^*}\right) g'\left( {n_1}^*\right) . \end{aligned}$$

(70)

Rearranging and solving for ${v_2}$ as a function of ${v_1}$ gives us the curve

$$\begin{aligned} {v_2} = h({v_1}) =&{} \frac{{v_1} \left[ {m^*} g'\left( {n_1}^*\right) + \left( {m^*}\right) ^2 f'\left( {m^*}\right) g'\left( {n_1}^*\right) \right] }{{n_1}^* \ f'\left( {m^*}\right) \left[ {v_1} {m^*} g'\left( {n_1}^*\right) - {n_1}^* g'\left( {n_1}^*\right) - 1\right] } \\ &{}- \frac{{m^*} f'\left( {m^*}\right) + {n_1}^* g'\left( {n_1}^*\right) + {m^*} {n_1}^* \ f'\left( {m^*}\right) g'\left( {n_1}^*\right) }{{n_1}^* \ f'\left( {m^*}\right) \left[ {v_1} {m^*} g'\left( {n_1}^*\right) - {n_1}^* g'\left( {n_1}^*\right) - 1\right] }. \end{aligned}$$

(71)

While Eq. 71 does not simplify further, it should be noted that all quantities are constants except for ${v_1}$ and ${v_2}$, assuming that all other parameters are held constant. Thus, ${v_2} = h({v_1})$ from Eq. 71 gives us the boundary of the region in ${v_1}$–${v_2}$ space where ${E_5}$ satisfies the third Jury condition (see Fig. 4a).

We will briefly note that this curve does yield both a vertical and horizontal asymptote. The vertical asymptote occurs at

$$\begin{aligned} {v_1} = {v_1}^* = \frac{{n_1}^* g'\left( {n_1}^*\right) + 1}{{m^*} g'\left( {n_1}^*\right) } \end{aligned}$$

(72)

and the horizontal asymptote at

$$\begin{aligned} {v_2} = {v_2}^* = \frac{{m^*} f'\left( {m^*}\right) + 1}{{n_1}^* \ f'\left( {m^*}\right) }. \end{aligned}$$

(73)

Equation 71 therefore has two branches. The left branch is valid for ${v_1} < {v_1}^*$ and ${v_2} > {v_2}^*$, and the right branch is valid for ${v_1} > {v_1}^*$ and ${v_2} < {v_2}^*$, with ${v_1}^*$ and ${v_2}^*$ given by Eqs. 72 and 73. However, note that

$$\begin{aligned} {v_2}^* = \frac{{m^*} f'\left( {m^*}\right) + 1}{{n_1}^* \ f'\left( {m^*}\right) } = \frac{{m^*}}{{n_1}^*} + \frac{1}{{n_1}^* f'\left( {m^*}\right) } < \frac{{m^*}}{{n_1}^*}. \end{aligned}$$

(74)

Recall by condition two from Table 2 that case a requires ${v_2} > \frac{{m^*}}{{n_1}^*}$. Therefore, the right branch of Eq. 71, for which ${v_2}< {v_2}^* < \frac{{m^*}}{{n_1}^*}$ is not a valid solution for this parameter regime, and so we only consider the left branch of Eq. 71 when determining the region of stability for ${E_5}$. On one side of the left branch of ${v_2} = h({v_1})$, ${E_5}$ is stable, and on the other side of the curve ${E_5}$ is unstable.

Case B

In case b, we have the same interior fixed point ${E_5}$ as in case a. For the first Jury condition, we see that the third condition from Table 2 is reversed, and

$$\begin{aligned} 1 - \tau + \Delta = {x_5} \ {y_5} \ f'\left( {m_5}\right) g'\left( {n_5}\right) \left( 1 - {v_1} {v_2} \right) < 0. \end{aligned}$$

(75)

The first Jury condition is never satisfied. By failing the first Jury condition, we know there must be at least one eigenvalue ${\lambda _1} > 1$. We can confirm this fact by examining the shifted characteristic polynomial in Eq. 62. By our above analysis, we know the constant term $1 - \tau + \Delta < 0$, so no matter the sign of $2 - \tau$, the coefficient of the linear term, there will be one sign change between consecutive coefficients in Eq. 62. Therefore we have one positive real root ${\mu _1}$, and thus only one eigenvalue ${\lambda _1} > 1$.

We have already demonstrated that in general, if the eigenvalues associated with an interior fixed point are real, they will both be positive. As we have one real $\lambda$, both must be real, as complex eigenvalues come in pairs. We therefore have ${\lambda _1} > 1$, $0< {\lambda _2} < 1$, and thus ${E_5}$ is a saddle point in case b.

Case C

In case c, we have one interior fixed point ${E_6}$. As in case b, the first Jury condition is never satisfied as

$$\begin{aligned} 1 - \tau + \Delta = {x_6} \ {y_6} \ f'\left( {m_6}\right) g'\left( {n_6}\right) \left( 1 - {v_1} {v_2} \right) < 0 \end{aligned}$$

(76)

using the third condition from Table 2. Precisely as before, we will use Descartes’ Rule of Signs to bound our eigenvalues. We have at least one eigenvalue ${\lambda _1} > 1$, by the results of the first Jury condition. To determine whether both eigenvalues are larger than $+1$, we examine the shifted characteristic polynomial (Eq. 62). As before, we know the constant term $1 - \tau + \Delta < 0$, so no matter the sign of $2 - \tau$, there is one sign change between the polynomial coefficients. Thus, there is one positive real root ${\mu _1}$, and therefore only one root ${\lambda _1} > 1$. Given one real $\lambda$, the other eigenvalue must be real as well, and we know all real eigenvalues are positive for the interior fixed points. Thus, ${\lambda _1} > 1$, $0< {\lambda _2} < 1$, and ${E_6}$ is a saddle point for case c.

Case D

In case d, we have the same interior fixed point ${E_6}$ as in case c. Using the third condition from Table 2, we now see

$$\begin{aligned} 1 - \tau + \Delta = {x_6} \ {y_6} \ f'\left( {m_6}\right) g'\left( {n_6}\right) \left( 1 - {v_1} {v_2} \right) > 0, \end{aligned}$$

(77)

and the first Jury condition is always satisfied in this case. The second Jury condition is also always satisfied.

Rather than checking the third Jury condition $\Delta < 1$, we instead show that both eigenvalues are real and bounded between 0 and $+1$. This statement is evidently sufficient to classify the fixed point, and proving realness eliminates the need to check the third Jury condition, which tests complex eigenvalues.

Consider the generalized eigenvalues of the Jacobian matrix evaluated at ${E_6}$, which are the roots of the characteristic polynomial Eq. 59:

$$\begin{aligned} \lambda _{1,2} = \frac{\tau \pm \sqrt{\tau ^2 - 4\Delta }}{2}. \end{aligned}$$

(78)

If $\tau ^2 - 4\Delta > 0$, then both eigenvalues are real. Using the trace from Eq. 63 and the expanded form of the determinant in Eq. 64, after combining like terms we find

$$\begin{aligned} \tau ^2 - 4\Delta =&\; \left[ {x_6} \ f'\left( {m_6}\right) \right] ^2 + \left[ {y_6} \ g'\left( {n_6}\right) \right] ^2 - 2 \ {x_6} \ {y_6} \ f'\left( {m_6}\right) \ g'\left( {n_6}\right) \\ &{}+ 4 {v_1} {v_2} {x_6} {y_6} \ f'\left( {m_6}\right) g'\left( {n_6}\right) \\ =&\; \left[ {x_6} \ f'\left( {m_6}\right) - {y_6} \ g'\left( {n_6}\right) \right] ^2 \\ &{}+ 4 {v_1} {v_2} {x_6} {y_6} \ f'\left( {m_6}\right) g'\left( {n_6}\right) > 0, \end{aligned}$$

(79)

as both terms are clearly positive under the conditions of this case. In fact, this statement holds for ${E_6}$ for any parameter regime, as we are not using conditions pertaining to a particular case, only conditions that are true at ${E_6}$. Thus, both eigenvalues are real for all cases where ${E_6}$ exists.

Having proved the eigenvalues are real, we again use Descartes’ Rule of Signs to bound the eigenvalues. By our earlier result from Appendix C.3 (General results for 𝜏 and Δ), since our eigenvalues are real we must have two positive eigenvalues, $\lambda _{1,2} > 0$.

We will use the shifted characteristic polynomial (Eq. 62) to determine whether $\lambda _{1,2} > 1$. By the first Jury condition, we know the constant term $1 - \tau + \Delta > 0$. The coefficient of the linear term

$$\begin{aligned} 2 - \tau = -{x_6} f'\left( {m_6}\right) - {y_6} g'\left( {n_6}\right) > 0, \end{aligned}$$

(80)

as both derivatives are negative at ${E_6}$. All coefficients of Eq. 62 are therefore positive, so there are zero sign changes and zero roots $\mu > 0$, meaning there are no eigenvalues $\lambda > 1$ (by Eq. 61). Thus, $0< \lambda _{1,2} < 1$, and ${E_6}$ is a stable node in case d.

Case E

In case e, we have both interior fixed points ${E_5}$ and ${E_6}$. We will first consider ${E_5}$. It is clear that

$$\begin{aligned} 1 - \tau + \Delta = {x_5} \ {y_5} \ f'\left( {m_5}\right) g'\left( {n_5}\right) \left( 1 - {v_1} {v_2} \right) > 0 \end{aligned}$$

(81)

by the conditions from Tables 2 and 4, so the first Jury condition is always satisfied for ${E_5}$ in case e. The second Jury condition is also always satisfied for ${E_5}$ in case e.

The third Jury condition is $\Delta < 1$. The expression for the third Jury condition is precisely the same as Eq. 68 from Appendix C.4 (Case A), and indeed the entirety of the analysis is identical to analysis of the third Jury condition for ${E_5}$ from case a. We do not repeat it here, and refer the reader back to Appendix C.4 (Case A) for details. Instead, we will simply recall that we can express the boundary of the region of stability for ${E_5}$ in terms of ${v_1}$ and ${v_2}$. This boundary yields a rational function, given by Eq. 71, with two branches. As the bounds on ${v_2}$ are the same between case a and case e (see Table 2), only the left branch is valid for determining the boundary of the region of stabiliy in case e. To one side of this branch, ${E_5}$ is stable, and on the other side ${E_5}$ is unstable. See Fig. 4b for illustration of the region of stability for ${E_5}$ in case e.

We now turn to the second interior fixed point, ${E_6}$. We see that

$$\begin{aligned} 1 - \tau + \Delta = {x_6} \ {y_6} \ f'\left( {m_6}\right) g'\left( {n_6}\right) \left( 1 - {v_1} {v_2} \right) < 0, \end{aligned}$$

(82)

by the conditions from Tables 2 and 4, so the first Jury condition is never satisfied for ${E_6}$ in case e.

As should now be familiar, we will use Descartes’ Rule of Signs to bound the eigenvalues. By the results of analyzing the first Jury condition, at least one eigenvalue is larger than $+1$. Examining the shifted characteristic polynomial (Eq. 62), we see there is only one sign change between the polynomial’s coefficients no matter the sign of $2 - \tau$, as the constant term $1 - \tau + \Delta < 0$. There is one positive real root ${\mu _1}$, and thus only one eigenvalue ${\lambda _1} > 1$. The other eigenvalue ${\lambda _2}$ must also be real, and as all real eigenvalues are positive for interior fixed points, we know ${\lambda _1} > 1$, $0< {\lambda _2} < 1$. Therefore, ${E_6}$ is a saddle point for case e.

Case F

In case f, we again have both interior fixed points ${E_5}$ and ${E_6}$. We first consider ${E_5}$.

Using the conditions and relationships in Tables 2 and 4,

$$\begin{aligned} 1 - \tau + \Delta = {x_5} \ {y_5} \ f'\left( {m_5}\right) g'\left( {n_5}\right) \left( 1 - {v_1} {v_2} \right) < 0. \end{aligned}$$

(83)

The first Jury condition is never satisfied for ${E_5}$ in case f. Considering the shifted polynomial Eq. 62, there is one sign change between coefficients of the polynomial and therefore one root ${\mu _1} > 0$. Thus, we have only one eigenvalue ${\lambda _1} > 1$. As before, the second eigenvalue ${\lambda _2}$ must be real given that ${\lambda _1}$ is real, and we know all real eigenvalues of the interior fixed points are positive. We find that ${\lambda _1} > 1$, $0< {\lambda _2} < 1$, and ${E_5}$ is a saddle point for case f.

Finally, we consider ${E_6}$ for this final case. The first Jury condition is obviously always satisfied as

$$\begin{aligned} 1 - \tau + \Delta = {x_6} \ {y_6} \ f'\left( {m_6}\right) g'\left( {n_6}\right) \left( 1 - {v_1} {v_2} \right) > 0 \end{aligned}$$

(84)

for ${E_6}$ in case f, by the conditions from Tables 2 and 4. The second Jury condition is always satisfied.

As with case d, we will bound the eigenvalues between 0 and $+1$ rather than checking the third Jury condition. Recall that in the stability analysis for case d, we demonstrated that both eigenvalues associated with ${E_6}$ are real for all cases where ${E_6}$ exists. Thus, we already know our eigenvalues to be real for ${E_6}$ in case f.

Employing the oft-used result a final time, we know that since the eigenvalues for ${E_6}$ are real, we have $\lambda _{1,2} > 0$. Identically to the analysis for case d, we consider the shifted characteristic polynomial Eq. 62. As the first Jury condition is always satisfied for ${E_6}$ in this case, the constant term $1 - \tau + \Delta > 0$. It is evident that the coefficient of the linear term

$$\begin{aligned} 2 - \tau = -{x_6} f'\left( {m_6}\right) - {y_6} g'\left( {n_6}\right) > 0 \end{aligned}$$

(85)

as well. There are zero sign changes between coefficients and therefore zero roots $\mu > 0$. Thus, there are no eigenvalues $\lambda > 1$ (by Eq. 61). With this result, we see that $0< \lambda _{1,2} < 1$, and ${E_6}$ is a stable node in case f.

The results of the above stability analyses for all six cases with one or more interior fixed point in the first quadrant are summarized in Table 3. This appendix completes our analytic stability analysis; analyses of the global bifurcations present in cases a and e were performed numerically.