Seasonal dynamics of a generalist and a specialist predator on a single prey

In ecological communities, the behaviour of individuals and the interaction between species may change between seasons, yet this seasonal variation is often not represented explicitly in mathematical models. As global change is predicted to alter season length and other climatic aspects, such seasonal variation needs to be included in models in order to make reasonable predictions for community dynamics. The resulting mathematical descriptions are nonautonomous models with a large number of parameters, and are therefore challenging to analyze. We present a model for two predators and one prey, whereby one predator switches hunting behaviour to seasonally include alternative prey when available. We use a combination of temporal averaging and invasion analysis to derive simplified models and determine the behaviour of the system, in particular to gain insight into conditions under which the two predators can coexist in a changing climate. We compare our results with numerical simulations of the temporally varying model.


Introduction
In predator-prey communities, population dynamics strongly depend on the way in which a prey responds to a predator, i.e., how many prey are killed and how frequently [9,23]. It is therefore important to understand the behaviour of the predator in a system in order to understand the dynamics of the system. Species behaviour often varies in response to seasonal phenomena [3,21]. For example, resource availability may alter how a predator hunts for food. In this paper, we consider two distinct behaviours of predators. When a predator specializes in hunting a single species and other food sources are negligible, we classify them as specialists. There are also instances where a predator may prefer a certain species if available, but as this species becomes hard to find, the predator will switch to alternative food sources. We classify these predators as generalists. Changes in season length due to climate change and global warming is expected to have significant impact on prey availability, and consequently on the ecological dynamics of predator-prey systems, especially in the presence of generalist predators.
A change in predation behaviour of the great horned owl (Bulbo virginialis) in the boreal forest, an area which is susceptible to climate change [1], is documented empirically [17]. The authors observe that gut content of the great horned owl indicates a change from a specialist in the winter to a generalist in the summer. In the winter, the snowshoe hare (Lepus americanus) represents a high percentage of the diet, and therefore the owl is a specialist predator of the hare in this season. In the summer, the hare is still prominent in the diet of the owl, but there is a significant amount of other species present. This indicates that the behaviour of the owl has switched to that of a generalist predator. This scenario of a behaviour switch has previously been modelled by [29]. The authors divided the year into two seasons and modelled the owl as a generalist in the summer and a specialist in the winter. They found that relatively small changes in summer season length can have a profound impact on the system. In particular, the predator can drive prey to extinction, there can be coexisting stable steady states, and there can be large-amplitude limit cycles coexisting with a stable steady state.
In this paper, we extend the two-species model of [29] to include the most important predator of the hare, the Canadian lynx (Lynx canadensis). The lynx is a specialist predator of the hare whose behaviour does not vary throughout the year [17]. Therefore, in the seasonal model the lynx will behave identically in both seasons.
We model the three-species system with a set of ordinary differential equations (ODEs), where we express the rate of change in a species density as a function of the density itself. To accurately model the seasonal dependence of predation, we create one set of ODEs for the summer season and a separate set for the winter season. The differences between these sets of equations are mainly due to the predators and their behaviour: the owl is a generalist in the summer and a specialist in the winter. We model the different behaviours by using two different functional responses [12]. The functional response represents the rate of decrease of the prey population due to predation. The commonly used Holling type II (for specialist predators) and III (for generalist predators) were derived mathematically and compared to empirical data ( [12,13] and [22], respectively). The seasonally dependent model we obtain is difficult to study analytically. The model is essentially non-autonomous, periodically and discontinuously forced. This renders any type of analytic study extremely challenging. To compensate, we take the annual average of the seasonal equations and study the resulting autonomous model. In particular, we are interested in stable coexistence of the species. In the framework of ODEs, this corresponds to a stable equilibrium or limit cycle of the system in which all components are strictly positive.
Although it is possible to derive formulas to determine regions in parameter space in which we have stable coexistence, these are too complicated to understand the biological implications. Therefore, we use a different technique that can predict when there is stable coexistence, but does so in a way that sheds light on the biological mechanisms at play, namely invasion analysis. Rather than studying coexistence in the three-species model directly, we ask whether one of the predators can invade (i.e., increase when rare) if the other two species coexist stably, for example at a steady state or a periodic orbit. With linear stability analysis of ODEs, we find invasion conditions, i.e., conditions under which each predator can simultaneously invade the other two species. This scenario is called mutual invasion. It is a longstanding tenet in population biology that "mutual invasion implies coexistence" [27], which has found numerous applications; see, e.g., [5,24,28]. However, explicit scenarios have been found in which mutual invasion does not imply coexistence [4,26]. We will test this hypothesis in our model. We note that the species can coexist stably at a steady state or at an oscillatory state, and the invasion conditions may not specify which coexistence we obtain.
In most of our analysis, we highlight and distinguish the parameter T , which represents in our model the proportion of the year corresponding to summer, and we discuss the effects of this parameter on such ecologically relevant issues as species coexistence and extinction. In this way, it is possible to make predictions on the potential effects of global warming and climate change on these issues. Most climate change scenarios predict that the length of the summer season will increase, in particular in northern climates where our study system is located.
The paper is organized as follows. In section 2, we derive our model and discuss the significance of the parameters, as well as the well-posedness of the model. In section 3, we use the tools and techniques of the qualitative theory of ODEs, including invasion analysis, to derive information regarding equilibria, periodic orbits, stability and bifurcations in our model. In section 4 we complement our theoretical analysis from section 3 with numerical analysis of the model. In particular we will present results on the positivity of the coexistence steady state, its stability, and relate these issues to predictions made from our invasion analysis. We will include a numerical comparison between the averaged and the seasonal model dynamics in our discussion in section 5.
We model the population dynamics of the snowshoe hare (Lepus americanus) and two of its predators, the Canadian lynx (Lynx canadensis) and the great horned owl (Bubo virginialis). We extend the model in [29] that did not include the lynx, which is the most important predator of the showshoe hare. Because the hunting behaviour of the owl changes seasonally, we divide the year into two distinct seasons (summer and winter) and we define appropriate systems of differential equations for each season. We denote time by τ and the population density of the hare by N (·), that of the lynx by P (·) and that of the owl by Q(·). The hare is assumed to grow logistically throughout the summer, whereas it does not reproduce in the winter. The lynx is a specialist predator year round. We use the standard Holling type II functional response [12,13] to model its impact on the prey and its resulting growth rate. In the absence of the hare, the lynx declines exponentially. The owl is a specialist predator in the winter and a generalist predator in the summer. We model its impact and growth in the winter also by a Holling type II functional response, but in the summer, we use a type III functional response, which better captures the behaviour of generalist predators [22,29]. During the winter, the owl declines exponentially in the absence of the hare. During the summer, the owl has an additional growth term that reflects the presence of alternative prey. Hare mortality in the winter occurs exclusively through predation.
We denote subsequent years by n, n + 1, . . . and the fraction of the year that is summer by T . Then τ ∈ [n, n + T ) corresponds to the summer in year n while τ ∈ [n + T, n + 1) corresponds to winter. Our seasonal model in dimensional form is given by Summer : Winter : τ ∈ [n + T, n + 1), n ∈ Z, (1.1) Population densities are continuous from the end of one season to the beginning of the next. Parameter r is the intrinsic hare growth rate and K the carrying capacity. Parameters c and α are the kill rates in the type II functional responses; b and β are the corresponding half-saturation constants. The kill rate of the type III functional response is a; the corresponding half-saturation constant is B 2 . The conversion efficiency of prey into predator biomass is f for the lynx and h for the owl. Parameter s is the maximal growth rate of the owl on their alternative prey, and v is the corresponding density-dependence. Finally, m and u are per capita death rates for the lynx and owl, respectively. We summarize these parameters, their interpretations and their units in table 1.
It is clear that this model, while still only a crude approximation to reality, is too complex for a complete qualitative analysis. The model is temporally periodic, it contains a large number of parameters (namely 15), and the vector field is not continuous at the switches of the seasons. Our first step at reducing model complexity is a standard nondimensionalization of (1.1) by means of the following change of variables: We drop the overbars to simplify notation and write the rescaled system for hare (x), lynx (y) and owl (z) with the remaining 12 parameters as Summer t ∈ [rn, r(n + T )), n ∈ Z, Winter t ∈ [r(n + T ), r(n + 1)), n ∈ Z, (1.2) While the reduction of parameters helps, the model is clearly still too complex for a complete qualitative analysis. The summer model alone (i.e., setting T = 1) contains a two-dimensional invariant subsystem of hare and owl dynamics (i.e., y = 0) that was independently analyzed in [8]; see [19] for a related model. It exhibits periodic orbits, bistability, and nonlocal bifurcations.
Nonetheless, our goal is to gain insights into the dynamics of the three populations in a seasonally varying environment. We continue to simplify model (1.2) by taking a weighted average of summer and winter equations, with respective weights rT (length of summer) and r(1 − T ) (length of winter). Averaging theory is, of course, well developed and applies to periodic vector fields with a small parameter [10]. Formally, we do not identify a small parameter in our model. However, previous studies have shown that this seasonal averaging can often predict non-averaged dynamics quite well [15,29], although these results certainly depend on chosen parameter values. Tyson and Lutscher found that for biologically reasonable parameter values in the hare-owl model (with the lynx absent), their averaged model either slightly over-predicted averaged densities, or it was virtually indistinguishable from averaged densities [29]. We speculate that one reason for the unreasonable effectiveness of temporal averaging here is that the birth-death dynamics of the two predators are typically somewhat slower than the yearly time scale, although not necessarily much smaller.
After averaging, we obtain the equationṡ This averaged model is the main focus of our analysis here. We will use several strategies to determine as much as possible of the qualitative behaviour of this model. We then numerically compare how well the averaged model predicts the dynamics of the seasonal model (1.2). Before we start the qualitative analysis of the averaged model, we establish some basic properties of the averaged system. Proof. The vector field is smooth in the positive orthant, hence, solutions exist at least locally and are unique. When x = 0 thenẋ = 0, and the same is true for the other two variables. Hence the nonnegative orthant is invariant. Sinceẏ = 0 whenever y = 0 andż = 0 whenever z = 0, each of the coordinate planes, where one of the state variables is zero, is invariant. When x > 1, thenẋ < 0. Hence, for any δ > 0 and any nonnegative initial condition, the x-coordinate of the solution satisfies lim t→∞ x(t) ≤ 1 + δ. In the invariant plane {z = 0}, the system is the well-known Rosenzweig-MacArthur model. Hence, there is some positive constant Y (depending on parameters but not on initial conditions), so that each nonnegative solution in this plane satisfies lim t→∞ y(t) ≤ Y . Since the first equation in (1.3) is decreasing in z and the second equation in (1.3) is increasing in x, the first two components of the solution of (1.3) with any nonnegative initial condition (x 0 , y 0 , z 0 ) are bounded above by the solution with initial condition (x 0 , y 0 , 0). Hence, lim t→∞ y(t) ≤ Y holds not only in the {z = 0} plane but also in the orthant. Finally, we turn to the equation for z. Since x is bounded by unity and the middle term in the first parentesis is bounded independently of z, the conditions in the Lemma are sufficient for there to exist some Z > 0 such that for z > Z, we haveż ≤ 0. Hence, the region where 0 ≤ x ≤ 1, 0 ≤ y ≤ Y and 0 ≤ z ≤ Z is compact and forward invariant.

Qualitative analysis of the averaged model
At first glance, the qualitative analysis of system (1.3) should simply follow established techniques. In fact, many ecological models exist for three-species dynamics; see, e.g., [11,14]. One difficulty arises from the still large number of parameters that reflects having two seasons. Another difficulty arises from the dynamics in one of the invariant planes alone being quite complex [8]. We begin our analysis with some relatively simple observations and then use reduction to the invariant planes combined with invasion analysis (see introduction) to obtain insights into the qualitative behaviour of the model.
2.1. The coexistence state and its stability. Clearly, the trivial state (x * , y * , z * ) = (0, 0, 0) is a steady state of system (1.3). Linearization at this state yields a diagonal Jacobi matrix with eigenvalues Hence, the zero state is unstable. More specifically, the first eigenvalue is related to the x coordinate.
Since it is positive, the hare can grow when all species are at low densities. The second eigenvalue is related to the y coordinate. Since it is negative, the lynx cannot grow when all species are at low densities. From the third eigenvalue, which corresponds to the z coordinate, we see that the owl can grow when all species are at low densities, provided that In terms of invasion analysis, we say that the hare will successfully invade when the two other species are absent, the lynx will never successfully invade, and the owl can invade provided (2.2) holds. A second, relatively simple steady state is the hare-only state (x * , y * , z * ) = (1, 0, 0). The associated Jacobi matrix is upper triangular with eigenvalues Hence, this state is locally asymptotically stable if and only if and These two conditions can again be interpreted by invasion analysis.
The computations that lead to this result are fairly straight forward. Equally straightforward is the condition that x * > 0 if and only if f > m. Formal conditions for positivity of the other two components can be written down but are much more cumbersome and do not lead new insights. Later, we will evaluate the conditions numerically. We call this state the coexistence state.
The linearization of (1.3) at the coexistence state is given by the Jacobi matrix where If parameter values were known, we could evaluate whether the coexistence state has all positive components and whether it is stable. Unfortunately, there are too many parameters and for some of those, reasonable ranges are too large or unknown. We also know that the invariant two-dimensional subsystem of hare and owl only can have complicated dynamics, including bistability [8,29], so that linear stability analysis will not give us the complete picture. For those reasons, we will combine invasion analysis with the study of the two invariant two-species submodels to obtain insights about the full three-species model. More precisely, we will consider each of the two-species subsystems at its attractor and use invasion analysis to find conditions under which the missing species can successfully invade, i.e., grow when rare. If invasion conditions for all species are satisfied simultaneously, we say the system exhibits mutual invasion. It is an old tenet of theoretical ecology that mutual invasion implies coexistence [27], but it is known that this approach may fail [2]. We will see whether and to what extent this principle holds in our system. Moreover, we investigate how the dynamics of the two-species subsystems influence the qualitative behaviour of the full system. Since the hare-lynx system is well understood, we begin by analyzing the owl-invasion conditions in that system.

Owl Invasion.
In the absence of the owl, system (1.3) becomes the classical Rosenzweig-MacArthur model [23] whose dynamics are well understood. If there exists a unique positive steady state for system (2.7), given by We note that the y-component is positive precisely when condition (2.3) is satisfied. Therefore, as the hare-only state loses stability via a transcritical bifurcation, the steady state (2.9) enters the positive quadrant of the (x, y)-plane.
The steady state (2.9) is stable when When the inequality is reversed, the positive steady state is unstable and there exists a unique positive and globally stable periodic orbit [16] that arises through a supercritical Hopf bifurcation where the inequality turns into an equality. Hence, there are only two scenarios at which we have to perform our invasion analysis for the owl: at the state in (2.9) and at a limit cycle.
Case 1: Owl invasion at steady state. When condition (2.10) is satisfied, the steady state (2.9) is stable in the hare-lynx plane of (2.7). After we linearize, the equation for the owl decouples from the other two. It is given byż For successful owl invasion, we require that z = 0 be unstable in (2.11). Hence, the owl invasion condition is Case 2: Owl invasion at a periodic orbit. When (2.10) is reversed, there is a globally stable periodic orbit of the form x * (t), y * (t), 0 with some period T . Invasion analysis at a periodic orbit asks the exact same question as invasion analysis at a steady state, namely, whether a species (in this case the owl) can grow when at low density in the existing community (in this case hare and lynx). Since the existing community is oscillatory, the linearized equation for the species at low density is a periodically forced equation. Whether the invading species can grow depends on the conditions that it encounters during one period. If the invading species can grow when rare, one expects a periodic orbit with all species present to emerge.
When we linearize our model along the hare-lynx orbit, the owl equation again decouples. It is given byż We evaluate the stability of z = 0 in (2.13) by means of the Floquet multiplier [6] 14) The owl can invade the periodic orbit of the hare-lynx system if In the absence of the lynx, system (1.3) becomeṡ which, up to a scaling of parameters, is the hare-owl model of Tyson and Lutscher [29]. The nonzero steady states of this equation are given by the roots of a complicated cubic polynomial, so that explicit calculations are tedious and uninformative [29]. Instead, we use invasion analysis on this twodimensional submodel to find conditions for the existence of steady states. The trivial state is unstable as in the three-dimensional model. The hare-only state is stable if T < T * * ; see (2.4), but unstable to owl invasion if T > T * * . The owl-only state and unstable to hare invasion if the inequality is reversed. This invasion condition is a quadratic equation in T , which can be solved explicitly to find the range of values of T , where the hare can invade the owl-only state. The result of invasion analysis is summarized by the following theorem.
Theorem 2.2. Define T * as in (2.2) and T * * as in (2.4). Then T * * < T * . The roots of (2.18) are given by If T ± are real, we have T * * < T * < T − < T + . If 0 < T < T * * , the hare-only state is globally asymptotically stable in the hare-owl system (2.16). If T > T * * , the hare-only state is unstable and the owl can invade it. If T > T * , both species can invade the trivial state independently. For T * * < T < T − , each species can invade the other at the single-species steady state. For T − < T < T + , the hare cannot invade the owl-only state. For T + < T < 1, we have again mutual invasion.
This theorem is a significant extension of the results in Tyson and Lutscher [29]. We illustrate its statements in Figure 1, where we plot the numerically calculated steady-state values for hare and owl. There are two cases. In the right plot, there is a stable coexistence state for all values T > T * * ; in the left plot, the hare goes extinct for some T inside the interval (T − , T + ). Tyson and Lutscher [29] did not find the scenario in the right plot, nor did they realize that there is a second branch of coexistence for large enough T in the left plot.
The plots show regions of bistability: the owl-only state is stable for T − < T < T + , yet there exists a positive stable coexistence state at least for some values of T in that interval. In the right plot, stable coexistence is possible for all T > T * * . In the left plot, there is a pair of saddle-node bifurcations where the coexistence state disappears and the reappears as T increases. The critical value between these two scenarios cannot be calculated explicitly, but a combination of algebraic geometry and bifurcation analysis can be used to obtain implicit equations that determine the critical case where the two saddle-node bifurcations touch. We can, however, give one simple condition to ensure that there is no bistability possible in the hare-owl system, namely when the solutions ( We now come to the question of lynx invasion into the hare-owl system. In addition to steady states and bistability, Tyson and Lutscher found limit cycles in the hare-owl model, and even bistability between a limit cycle and a steady state [29]. Therefore, we consider two cases of invasion conditions of the lynx in the hare-owl model.  Hence, the lynx invasion condition at a steady state is (2.22) which is the same as (2.3), where x * = 1.
Case 2: Lynx invasion at a periodic orbit. When we linearize the system at a periodic orbit (x * (t) , 0, z * (t)) with period T , the lynx equation decouples and readṡ We use Floquet theory to obtain the lynx invasion conditions at a periodic orbit as In the next section, we use these results and insights about the two-species subsystems to study the three-species model in the positive orthant. We choose the two parameters T and s for illustration purposes. As mentioned previously, the impact of T is of great interest when interpreting the results in terms of global change since it represents the fraction of the year that is summer. We chose s as our second parameter for a number of reasons. The hare-lynx system is independent of T and s, so that the dynamics in that invariant subspace do not change throughout 2-parameter plane. Also, s is the parameter that we have the least information on, while for many other parameter values, estimated ranges can be found in the literature [29]. It also turns out that we can express the curves of T * , T * * , Even if introduced at high density, it will asymptotically approach extinction. Above the blue dashed curve, the owl can invade the hare-only state, but not the hare-lynx state nor the trivial state. Above the black curve, which representss, the owl can invade the hare-lynx state, but not trivial state. Above the blue solid curve, which represents T * , the owl can invade all possible steady states in the hare-lynx plane. Above the red curve, which represents T ± , the hare cannot invade the owl-only state, so that there is bistability in the hare-owl plane. Condition (2.20) corresponds to a horizontal line that touches the minimum of the red curve. and T ± explicitly for s = s(T ), which makes plotting these quantities easy. The explicit expressions are respectively. Furthermore, with the hare density at the hare-lynx steady state (2.9) denoted as x * = mb f −m < 1, the owl invasion condition (2.12) at the hare-lynx state for s as a function of T is From the expressions, it is obvious that s * * <s < s * < s ± . We illustrate these quantities in Figure  2. While the absolute locations of these curves change with parameter values, their relative locations do not. In the next section, we take this figure as the blueprint and we superimpose the region of three-species coexistence and its stability properties for difference scenarios.

Numerical results and comparison with the seasonal model
In this section, we evaluate the invasion conditions, the steady-state expression and its stability conditions numerically and display them in plots similar to that in Figure 2. We are particularly interested in whether mutual invasion implies coexistence and how the behaviour in the two-species subsystems in the invariant planes affects the behaviour of the three-species system in the positive orthant. We also compare the behaviour of the averaged model with simulations of the seasonal model. We consider different scenarios, depending on the dynamics of the two-species subsystems. We continue to use T and s as our two bifurcation parameters. 3.1. Globally stable coexistence steady state in the hare-lynx system. We begin with the simplest case where (2.8) and (2.10) are satisfied so that the hare-lynx subsystem has a globally stable coexistence state as in (2.9). The qualitative dynamics of the hare-lynx subsystem are independent of parameters T and s. We choose s small enough so that the hare-owl subsystem does not show bistability according to Theorem 2.2. In addition to the owl invasion condition at the hare-lynx state (2.25) we calculate the the lynx invasion condition at the hare-owl state; see Figure 3. For the chosen parameters, the owl cannot invade when T and s are low (region C), and the lynx cannot invade when T and s are high (region B), but there is mutual invasion in some intermediate range (region A). The green curve in this plot is exactly the same as the black curve in Figure 2. The coexistence state (2.5) has all positive components in the black region in the middle panel of Figure 3. This region exactly equals region A in the left plot. Finally, the eigenvalues at the coexistence state have all negative real parts in the black region in the right panel of Figure 3. Again, this region exactly equals region A in the left plot. (We observed that the regions are identical by superimposing the plots.) Hence, we see that mutual invasion implies (stable) coexistence here. We conjecture that the coexistence state is globally stable if it is positive, but we cannot prove this.

3.2.
Stable periodic orbit in the hare-lynx system. In the next simplest case, (2.8) is satisfied but (2.10) is not so that the hare-lynx subsystem has a globally stable limit cycle, independent of T and s. As above, we choose s small enough so that the hare-owl subsystem does not show bistability according to Theorem 2.2. We calculate the respective invasion conditions, evaluate whether the threespecies coexistence state has all positive components, and calculate the eigenvalues at that state; see Figure 4. We do not display the invasion conditions separately since they look qualitatively exactly like in Figure 3, left plot. However, we plot the two curves into the steady state (left) and stability (right) plots here. The region of mutual invasibility (i.e., the lynx invades the hare-owl steady state and the owl invades the hare-lynx periodic orbit) corresponds exactly to the region where a positive coexistence state exists (black region in the left plot). The positive coexistence state is unstable and has two complex conjugate eigenvalues with positive real part in the green region in the right plot. We expect a positive periodic orbit in this region. This orbit disappears in a Hopf bifurcation at the boundary between the green and black regions. The positive coexistence state is stable in the black region and disappears in a transcritical bifurcation at the boundary with the blue region where one positive real eigenvalue emerges.

3.3.
More complex dynamics in the hare-owl plane. All the cases considered thus far had very simple steady-state dynamics in the hare-owl plane, yet, it is known that more complex scenarios are possible. Tyson and Lutscher found that the hare-owl system alone can have multistability (three positive steady states of which two are locally stable) or a homoclinic loop bifurcation to a large limit cycle; see Figure 5B in [29]. Since that system is two dimensional, chaotic behaviour cannot occur. Steady states and periodic orbits in the hare-owl plane can both undergo transcritical bifurcations if the lynx invasion conditions (2.22) or (2.24) are satisfied. These conditions can be evaluated numerically. Rather than attempting a complete analysis of all potential scenarios (which is beyond the scope of this work because of the sheer number of parameters), we consider several interesting cases.
We begin with the case where the hare-owl subsystem shows multiple steady states ( Figure 5, left plot). The hare nullcline (blue) and the owl nullcline (red) intersect three times in the positive quadrant. The black curves show solutions that approach the left and the right steady state and reveal the saddle structure of the middle state. The lynx can invade the largest hare-owl state if its death rate is small enough, here m < 0.27. It can invade the largest two if m < 0.093 and all three if m < 0.011. However, we know that there can be at most one coexistence state with all three coordinates positive. If the lynx can invade all three states, what happens to the other two? The lynx death rate does not affect the hare-owl subsystem. It turns out that the lynx component of the coexistence state changes sign in this range of m ( Figure 5, right plot).
We illustrate what happens to lynx invasion along the upper positive branch of the coexistence state in Figure 5, right plot. This corresponds to the case where 0.093 < m < 0.27, i.e., the lynx can invade all three positive hare-owl states. For m = 0.25, the lynx invades from low density and the system settles at the stable coexistence point (Figure 6, top left). For m = 0.18, the system reaches the stable coexistence point by decaying oscillations. For m = 0.173, the lynx invades from low density, decreases the hare density so low that the lynx itself cannot survive. The hare-owl system settles at its lower stable state (bottom left). However, for the same value of m, all three species can persist in the system if initial conditions are chosen close to the coexistence state. The coexistence state can be stable or unstable, in which case a limit cycle exists (bottom right). Decreasing m further increases the amplitude of the limit cycle, so that eventually the hare and lynx drop below the threshold where the lynx can persist, and the hare-owl system settles at its lower stable state as in the bottom left plot.
Hence, the region of multistability in the hare-owl subsystem creates a region of multistability in the three-species system. There are two cases: two locally stable states may coexist (one with and one without the lynx), or one locally stable state (without the lynx) may coexist with a locally stable limit cycle (of all three species).
The second case that we consider also includes multistability in the hare-owl subsystem, but this time between a periodic orbit and a steady state (Figure 7, left plot). The nullclines intersect only once in the positive quadrant, and the steady state is locally stable. In addition, there is a locally stable limit cycle that arose through a homoclinic loop bifurcation; for details, see [8,29]. The hare density at the steady state is approximately 0.0417, so that the invasion condition for the lynx at this state requires m < 0.02. The lynx can invade the limit cycle according to condition (2.24). Numerical evaluation gives a Floquet multiplier of almost 6. Hence, we observe a limit cycle of coexistence of all three species in the absence of a coexistence steady state with positive densities for all species. In addition, we have multistability since the semitrivial state where the lynx is absent and the hare and owl are at their steady state, is also locally stable. As in the previous case, initial conditions can lead to the system moving from one basin of attraction to another and thereby lead to a surprising outcome. If the lynx is introduced in low density near the hare-owl steady state, it will decline to extinction because its death rate (m = 0.4) is much higher than the threshold for invasion (m = 0.02). However, if the lynx is  Figure 5 introduced at a very high density, it will decimate the hare initially and thereby push the entire system into the basin of attraction of the large limit cycle. At that limit cycle, the lynx can invade, and it will therefore stay in the system. For the given parameter values, the hare-lynx system (in the absence of the owl) will approach a globally stable positive steady state. The owl can invade this steady state from low density. Therefore, the owl will never be excluded from the system in this scenario. This latter example shows that three species may coexist indefinitely even in the absence of a positive steady state. In other words, we would not have found this coexistence periodic orbit from a local stability analysis of the (explicitly available) steady-state expression. This finding illustrates the power of invasion analysis. We note that periodic coexistence of a two-predator-one-prey system in the absence of a coexistence steady state has been proven before in a different system [25].

Discussion
Mathematical models for population dynamics in a seasonally varying environment are becoming more prominent and more important for studying the effects of a changing climate on biological communities. A particularly interesting aspect arises when the interaction between species changes with seasons. Our study is inspired by the snowshoe hare, Canadian lynx and great horned owl system in western Canada, where the owl acts as a generalist predator in the summer and as a specialist in the winter, while the lynx is a specialist predator year round [17]. Our model is an extension of the two-species model in [29], where the lynx was not considered, even though it is the dominant predator of the showshoe hare.  Even though we represented a continuously varying environment by simply distinguishing two seasons, our model is still a periodically forced system of equations. We use temporal averaging to derive an autonomous model that still is complicated since the model has in some sense "more" parameters than "normal" models due to the fact that two seasons are represented. We could show that our model has a unique coexistence steady state (with all three components nonzero), but the precise conditions for positivity and stability require numerical evaluation. Instead, we used the structure of the system with two invariant planes (the hare-owl and the hare-lynx subsystem) in combination with invasion analysis to obtain a number of more general qualitative results on our system. In particular, we considered the principle that "mutual invasion implies coexistence" [27] to find conditions of coexistence of the two predators.
We found that the length of the summer season (T ) can have a profound effect on system dynamics with various bifurcations occurring as T changes. General conclusions, however, are difficult to draw since the dynamics overall can be hugely varying. A longer summer season tends to benefit the generalist predator because it has additional food sources then. In particular, there is a threshold value of T , above which the generalist does not require the focal prey for persistence. It is unclear how likely this scenario is in reality, both in terms of T ever getting that large and winter mortality being modeled by a linear process. On the other hand, for the parameter values chosen, the specialist predator has an advantage at short season lengths. Intermediate values of T seem to promote coexistence between the two species. It was not our goal to analyze system behavior in all of parameter space since this would be impossible. Instead, the techniques that we presented translate to other seasonal population models. If and when estimates of parameter are available, our method allows a relatively inexpensive analysis of the averaged model, which can then be followed up with simulations of the seasonal model in parameter ranges that the averaged model identified as interesting.
It is an old tenet in population biology that two predator species cannot persist at a stable steady state on a single limiting resource [18], but they can coexist in a stable limit cycle [2]. Since one of our predators has an alternate food source during the summer, this tenet does not apply to our model, and we found indeed that coexistence of the two species is possible at a stable steady state as well as a stable limit cycle. We also found that mutual invasion implied coexistence when the dynamics in each of the two invariant subsystems were simple enough. When the dynamics in the hare-owl system showed bistability, some ecologically surprising outcomes were possible. In particular, an invasion attempt by the lynx could be successful in the short run but cause the system to move to the basin of attraction of another state, at which the lynx could not invade, so that it would eventually go extinct. Vice versa, an invasion attempt with a sufficient number of lynx could push the system from a non-invasible state to the basin of attraction of an invasible state and the lynx could persist in the long run. Such an initial invasion at large numbers may not be expected to occur naturally but could be used a a management tool (e.g., through a captive breeding and release program) to resettle a species that has become locally extirpated.
While our model included some variability between seasons, it did not include all that is possible. For example, it is conceivable that the specialist predation success varies with season when snow cover gives refuges to the snowshoe hare. We also did not include any winter mortality of the prey. Such extensions can easily be included into the model, and the process of averaging can still be applied. The analysis of the resulting model is complicated not only by the increased number of parameters but also by the fact that the dynamics on the invariant hare-lynx plane then depend on season length T as well.
The most obvious question that remains from our analysis is to what extent the averaged model captures and represents the dynamics of the seasonal model. From an abstract point of view, if there is a small parameter that represents the separation of time scales, averaging theory can be applied. If no obvious small parameter exists, there are no general theorems about the applicability of any kind of averaging. However, it is a typical observation in averaging theory that the averaged model provides a good approximation even beyond the case where it can be mathematically proven to be valid. While we cannot analytically identify a small parameter in our model, biological considerations indicate different time scales. While the hare population dynamics are relatively fast with several reproductive events within a single summer, lynx and owls reproduce only once a year so that their dynamics could be considered slow compared to seasonal changes. When we compared numerical simulations of the averaged model with the seasonal model, we found a good agreement for the simpler dynamics from Sections 3.1 and 3.2; see Figure 8. The solution of the averaged model (in black) is very close to the solution of the seasonal model (in colour). The corresponding analysis of the averaged model is in Figures 3 and 4.
With some more complex dynamical behaviour, the correspondence between the averaged and the seasonal model is not as perfect; see Figure 9 and corresponding Figure 5. As lynx mortality m decreases, the system passes through a Hopf bifurcation, but as limit cycles grow in amplitude, solutions get trapped into the semitrivial state where the lynx is extinct. Here, the solution to the averaged model captures the solution to the seasonal model at least for some transient period but tends to collapse to the semitrivial state earlier. The discrepancy might not be satisfying for mathematicians, yet, decent short-term predictions and transients are still highly valuable in applications in ecosystem management.
These findings are in line with previous results [15,29]. It is clear, however, that this method is limited and cannot work always. For example, it is known that periodic two-species competition systems possess properties that no autonomous competition system has [7,20]. More specifically to our system, if the dynamics were fast within each season and seasons changed slowly (i.e., the opposite of what happens), the system would approach a stable state within a season. If that state happens to exclude one of the predators, there would be no chance for this predator to return to the system in the subsequent season. Hence, the averaging technique could not possibly predict the correct result. Nonetheless, our study as well as others show that temporal averaging can give important insights and act to guide more detailed simulations of the full seasonal model.