Mathematics in Applied Sciences and Engineering https://ojs.lib.uwo.ca/mase
Volume 2, Number 2, June 2021, pp.103-122 https://doi.org/10.5206/mase/13569

SEASONAL DYNAMICS OF A GENERALIST AND A SPECIALIST PREDATOR

ON A SINGLE PREY

NOAH BOLOHAN, VICTOR LEBLANC, AND FRITHJOF LUTSCHER

Abstract. 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 mathe-

matical 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 com-

munity 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.

1. 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

Received by the editors 16 January 2021; accepted 14 April 2021; published online 16 April 2021.

2000 Mathematics Subject Classification. Primary 92D40, 34C23; Secondary 34A36.
Key words and phrases. predator-prey model, seasonal variation, temporal averaging, qualitative analysis, bifurcation,

invasion analysis.
V. LeBlanc was supported in part by a Discovery Grant from the Natural Resources and Engineering Council of

Canada (RGPIN-2016-04318).

F. Lutscher was supported in part by a Discovery Grant from the Natural Resources and Engineering Council of

Canada (RGPIN-2016-04795).

103



104 N. BOLOHAN, V. LEBLANC, AND F. LUTSCHER

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



SEASONAL DYNAMICS 105

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 :

τ ∈ [n,n + T), n ∈ Z,
dN

dτ
= rN

(
1 −

N

K

)
−

cNP

b + N
−

aN2Q

B2 + N2
,

dP

dτ
= f

cNP

b + N
−mP,

dQ

dτ
= h

aN2Q

B2 + N2
+

sQ

1 + vQ
−uQ,

Winter :

τ ∈ [n + T,n + 1), n ∈ Z,
dN

dτ
= −

cNP

b + N
−

αNQ

β + N
,

dP

dτ
= f

cNP

b + N
−mP,

dQ

dτ
= h

αNQ

β + N
−uQ.

(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 B2. 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 param-

eters (namely 15), and the vector field is not continuous at the switches of the seasons. Our first step



106 N. BOLOHAN, V. LEBLANC, AND F. LUTSCHER

Parameter Description Units

T proportion of the year that is summer dimensionless

r hare growth rate 1/year

K hare carrying capacity hare

c lynx maximum kill rate hare/lynx/year

b lynx half-saturation constant hare

f lynx conversion efficiency lynx/hare

m lynx death rate 1/year

a owl maximum kill rate (summer) hare/owl/year

B owl half-saturation constant (summer) hare

α owl maximum kill rate (winter) hare/owl/year

β owl half-saturation constant (winter) hare

h owl conversion efficiency owl/hare

u owl death rate 1/year

s owl growth rate on alternative prey 1/year

v owl self-interference on alternative prey 1/owl

Table 1. Parameters in the seasonal model (1.1).

at reducing model complexity is a standard nondimensionalization of (1.1) by means of the following

change of variables:

x =
N

K
,

y =
cP

rK
,

z =
aQ

rK
,

t = rτ,

b =
b

K
,

B =
B2

K2
,

a =
α

a
,

d =
β

K
,

f =
fc

rK
,

m =
m

r
,

h =
ha

r
,

g =
hα

r
,

s =
sa

r2K
,

v =
a

vKr
,

u =
u

R
.

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,

ẋ = x(1 −x) −
xy

b + x
−

x2z

B + x2
,

ẏ =
fxy

b + x
−my,

ż =
hx2z

B + x2
+

sz

v + z
−uz,

Winter

t ∈ [r(n + T),r(n + 1)), n ∈ Z,

ẋ =
xy

b + x
−

axz

d + x
,

ẏ =
fxy

b + x
−my,

ż =
gxz

d + x
−uz.

(1.2)

While the reduction of parameters helps, the model is clearly still too complex for a complete qual-

itative 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.



SEASONAL DYNAMICS 107

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 equations

ẋ = rT

(
x(1 −x) −

xy

b + x
−

x2z

B + x2

)
+ r(1 −T)

(
−

xy

b + x
−

axz

d + x

)
,

ẏ = r

(
fxy

b + x
−my

)
,

ż = rT

(
hx2z

B + x2
+

sz

v + z
−uz

)
+ r(1 −T)

(
gxz

d + x
−uz

)
.

(1.3)

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.

Lemma 1.1. System (1.3) is well-posed and preserves positivity. Furthermore, if u > max{ h
B+1

, g
d+1
}

then there is a compact, forward invariant, attractive set in the positive orthant. The x-axis, as well as

the planes {y = 0} and {z = 0} are invariant for the dynamics.

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

limt→∞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 limt→∞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 (x0,y0,z0) are bounded

above by the solution with initial condition (x0,y0, 0). Hence, limt→∞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. �

2. 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



108 N. BOLOHAN, V. LEBLANC, AND F. LUTSCHER

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

λ1 = rT > 0, λ2 = −rm < 0, and λ3 = r
(
Ts

v
−u
)
. (2.1)

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

T > T∗ =
uv

s
. (2.2)

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

λ1 = −rT, λ2 = r
(

f

b + 1
−m

)
and

λ3 = rT

(
h

B + 1
+
s

v

)
+
r(1 −T)g
d + 1

−ru.

Hence, this state is locally asymptotically stable if and only if

f < m(b + 1) (2.3)

and

T < T∗∗ =
u− g

d+1
h

B+1
+ s

v
− g

d+1

. (2.4)

These two conditions can again be interpreted by invasion analysis. If (2.3) is reversed, then the lynx

can invade the hare-only state, and if (2.4) is reversed, then the owl can invade that state. We will

come back to these conditions in Sections 2.2 and 2.3, respectively, where we will study steady states

where two of the three species are present. Here, we continue with the case that all three species have

nonzero densities.

Lemma 2.1. System (1.3) has a unique steady state where all densities are nonzero. It is given by

x∗ =
mb

f −m
,

y∗ = (b + x∗)

(
T

(
1 −x∗ −

x∗z∗

B + (x∗)2

)
−

(1 −T)az∗

d + x∗

)
,

z∗ =
−Ts

T h(x∗)2

B+(x∗)2
+

(1−T )gx∗
d+x∗

−u
−v.

(2.5)



SEASONAL DYNAMICS 109

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


J11 − rx
∗

b+x∗
− rT (x

∗)2

B+(x∗)2
− r(1−T )ax

∗

d+x∗

rfby∗

(b+x∗)2
r
(

fx∗

b+x∗
−m

)
0

2rT Bhx∗z∗

(B+(x∗)2)2
+

r(1−T )gdz∗
(d+x∗)2

0 J33


 , (2.6)

where

J11 = rT

(
1 − 2x∗ −

2Bx∗z∗

(B + (x∗)2)2

)
−
r(1 −T)adz∗

(d + x∗)2
−

rby∗

(b + x∗)2
,

J33 = rT

(
h(x∗)2

B + (x∗)2
+

sv

(v + z∗)2

)
+
r(1 −T)gx∗

d + x∗
−ru.

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.

2.2. Owl Invasion. In the absence of the owl, system (1.3) becomes the classical Rosenzweig–MacArthur

model [23]

ẋ = rTx(1 −x) −r
xy

b + x
,

ẏ = r

(
fxy

b + x
−my

)
,

z = 0,

(2.7)

whose dynamics are well understood. If

0 <
mb

f −m
< 1, (2.8)

there exists a unique positive steady state for system (2.7), given by

(x∗,y∗, 0) =

(
mb

f −m
,T(1 −x∗)(b + x∗), 0

)
. (2.9)

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.



110 N. BOLOHAN, V. LEBLANC, AND F. LUTSCHER

The steady state (2.9) is stable when

f(1 − b) −m(1 + b) < 0. (2.10)

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

ż =

[
rT

(
h(x∗)2

B + (x∗)2
+
s

v

)
+
r(1 −T)gx∗

d + x∗
−ru

]
z. (2.11)

For successful owl invasion, we require that z = 0 be unstable in (2.11). Hence, the owl invasion

condition is

T

(
h(x∗)2

B + (x∗)2
+
s

v

)
+

(1 −T)gx∗

d + x∗
−u > 0. (2.12)

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

ż =

[
rT

(
h(x∗(t))2

B + (x∗(t))2
+
s

v

)
+
r(1 −T)gx∗(t)
d + x∗(t)

−ru
]
z ≡ A(t)z. (2.13)

We evaluate the stability of z = 0 in (2.13) by means of the Floquet multiplier [6]

µ =
1

T

∫ T
0

A(t)dt. (2.14)

The owl can invade the periodic orbit of the hare-lynx system if

µ =
1

T

∫ T
0

[
rT

(
h(x∗(t))2

B + (x∗(t))2
+
s

v

)
+
r(1 −T)gx∗(t)
d + x∗(t)

−ru
]

dt > 0. (2.15)

2.3. Lynx Invasion. In the absence of the lynx, system (1.3) becomes

ẋ = rT

(
x(1 −x) −

x2z

B + x2

)
+ r(1 −T)

(
−
axz

d + x

)
,

y = 0,

ż = rT

(
hx2z

B + x2
+

sz

v + z
−uz

)
+ r(1 −T)

(
gxz

d + x
−uz

)
,

(2.16)



SEASONAL DYNAMICS 111

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 ex-

plicit calculations are tedious and uninformative [29]. Instead, we use invasion analysis on this two-

dimensional 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

(0,z∗(T)) =
(

0,
s

u
(T −T∗)

)
, T∗ =

uv

s
(2.17)

is stable if

T −
a

d
z∗(T)(1 −T) < 0 (2.18)

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

T± =
du

2as


−(1 − a(s + uv)

du

)
±

√(
1 −

a(s + uv)

du

)2
−

4a2sv

d2u


 . (2.19)

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 (2.18) are not real.

This is the case if the discriminant in (2.19) is negative, i.e., if(
1 −

a(s + uv)

du

)2
<

4a2sv

d2u
. (2.20)

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.



112 N. BOLOHAN, V. LEBLANC, AND F. LUTSCHER

Figure 1. Steady states of model (2.16) as a function of T. Blue and green solid

curves represent the hare and owl densities at the coexistence state(s), respectively.

The green dashed line represents the owl-only steady state. For T < T∗∗, the only

nontrivial steady state is the hare-only state (1, 0). Fixed parameters are B = 0.0625,

a = 0.2, d = 0.08, h = 0.07, u = 0.5, g = 0.07, v = 0.3333. We set s = 0.75 in the left

plot and s = 0.74 in the right.

Case 1: Lynx invasion at a steady state. When we linearize the system at the steady state (x∗, 0,z∗),

the linearized lynx equation decouples and reads

ẏ =

(
fx∗

b + x∗
−m

)
y. (2.21)

Hence, the lynx invasion condition at a steady state is

fx∗

b + x∗
−m > 0, (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 reads

ẏ =

(
fx∗(t)

b + x∗(t)
−m

)
y. (2.23)

We use Floquet theory to obtain the lynx invasion conditions at a periodic orbit as

1

T

∫ T
0

[
fx∗(t)

b + x∗(t)
−m

]
dt > 0. (2.24)

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∗∗,



SEASONAL DYNAMICS 113

Figure 2. Regions in the T-s parameter plane according to threshold values (2.25)

and (2.26). The blue dashed curve represents T∗∗. Below it, the owl cannot invade

the system under any circumstances. Even if introduced at high density, it will asymp-

totically 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 represents s̄, 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

s∗ =
uv

T
, s∗∗ =

v

T

(
u−

g(1 −T)
d + 1

−
Th

B + 1

)
, s± =

uv

T
+

du

a(1 −T)
, (2.25)

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

s̄ =
v

T

(
u−

gx∗(1 −T)
d + x∗

−
Th(x∗)2

B + (x∗)2

)
. (2.26)

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.

3. 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.



114 N. BOLOHAN, V. LEBLANC, AND F. LUTSCHER

Figure 3. Invasion conditions (left), coexistence region (middle) and stability of the

coexistence state (right) in the case that the hare-lynx submodel has a globally stable

coexistence state. Left: The lynx invades the hare-owl state (where it exists) in regions

A and C. The owl invades the hare-lynx state in regions A and B. Middle: The three-

species coexistence state has all positive components in the black region. Only the

owl term is non-positive in the magenta region, and only the lynx term is non-positive

in the cyan region. Right: All three eigenvalues have negative real part in the black

region, while only two have negative real part in the blue region. Fixed parameters are

r = 1, B = 0.0625, a = 0.2, d = 0.08, h = 0.07, u = 0.5, g = 0.07, v = 1
3
, b = 0.3,

f = 3.2 and m = 2.

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 three-

species 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



SEASONAL DYNAMICS 115

Figure 4. In the left plot, mutual invasion occurs between the red and grey curves,

and the steady state is positive in the black region. In the right plot, mutual invasion

occurs between the cyan and magenta curves, and the steady state is stable in the black

region. Fixed parameters are r = 1, B = 0.0625, a = 0.2, d = 0.08, h = 0.07, u = 0.5,

g = 0.07, v = 1
3
, b = 0.3, f = 3.2 and m = 1.6.

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



116 N. BOLOHAN, V. LEBLANC, AND F. LUTSCHER

0 0.2 0.4 0.6 0.8 1
0

0.05

0.1

0.15

0.2

0.25

0.3

hare

o
w

l

−0.5 0 0.5
0.1

0.11

0.12

0.13

0.14

0.15

0.16

0.17

lynx

o
w

l

Figure 5. Left plot: Multistability in the hare-owl phase plane. The hare nullcline

(blue) intersects the owl nullcline (red) three times in the positive quadrant. The

smaller and the larger intersection are locally stable, the middle intersection is a saddle.

Right plot: The owl and lynx components of the coexistence state as a function of

lynx death rate m. The intersections with the line where the lynx component is zero

correspond to the three positive steady states of the hare-owl subsystem. The lynx

component is zero when the lynx invasion condition is an equality, i.e., for m ≈ 0.27
(top), m ≈ 0.093 (middle) and m ≈ 0.011 (bottom). Parameters are B = 0.0025, a =
0.2, d = 0.08, h = 0.07, u = 0.5, g = 0.07, v = 1/3, b = 2, f = 1, s = 1/3, T = 0.65.

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



SEASONAL DYNAMICS 117

0 100 200 300 400 500
0

0.2

0.4

0.6

0.8

time

d
e

n
s
it
ie

s

0 100 200 300 400 500
0

0.2

0.4

0.6

0.8

time

d
e

n
s
it
ie

s

0 100 200 300 400 500

0

0.2

0.4

0.6

0.8

time

d
e

n
s
it
ie

s

0 100 200 300 400 500

0.2

0.25

0.3

0.35

0.4

0.45

0.5

time

d
e

n
s
it
ie

s

Figure 6. Time course of the system as the lynx invades a hare-owl state for m = 0.25

(top left), m = 0.18 (top right), m = 0.173 (bottom). The lynx initial conditions for

all plots except bottom right are 0.02. At the bottom right, initial conditions are near

the coexistence state. Parameters are as in 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].

4. 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 com-

munities. 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.



118 N. BOLOHAN, V. LEBLANC, AND F. LUTSCHER

0 0.2 0.4 0.6 0.8
0

0.02

0.04

0.06

0.08

0.1

hare

o
w

l

0 100 200 300 400 500
0

0.2

0.4

0.6

0.8

time

d
e

n
s
it
ie

s

Figure 7. Left plot: Multistability in the hare-owl plane. Right plot: successful

lynx invasion leads to limit cycle coexistence of the three species. Parameters are

B = 0.0025, a = 0.2, d = 0.1, h = 0.05, u = 0.4, g = 0.1, v = 1/3, b = 2, f =

1, m = 0.05, s = 1/3 and T = 0.414.

Figure 8. Comparing solutions of the averaged model (1.3) to those of the seasonal

model (1.2). Averaged solutions are plotted in black, and seasonal dynamics are plotted

in blue, red and green, representing hare, lynx and owl, respectively. Condition (2.10)

is satisfied in the first row, yielding a steady state in the averaged model. The condition

is reversed in the second row, yielding a periodic solution. Fixed parameters are r = 1,

B = 0.0625, a = 0.2, d = 0.08, h = 0.07, u = 0.5, g = 0.07, v = 1
3
, b = 0.3, f = 3.2,

s = 0.4 and T = 0.5. We set m = 2 in the first row and m = 1.7 in the second.

Even though we represented a continuously varying environment by simply distinguishing two sea-

sons, 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



SEASONAL DYNAMICS 119

Figure 9. Comparing solutions of the averaged model (1.3) to those of the seasonal

model (1.2). Averaged solutions are plotted in black, and seasonal dynamics are plotted

in blue, red and green, representing hare, lynx and owl, respectively. B = 0.0025, a =

0.2, d = 0.08, h = 0.07, u = 0.5, g = 0.07, v = 1/3, b = 2, f = 1, s = 1/3, T = 0.65.

m = 0.1753 in the first row, m = 0.17513 in the second and m = 0.1749 in the third.

Initial conditions for all rows are (x0,y0,z0) = (0.75, 0.01, 0.168).

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



120 N. BOLOHAN, V. LEBLANC, AND F. LUTSCHER

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



SEASONAL DYNAMICS 121

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.

References

[1] M. Allen, O. Boucher, D. Chanbers, J. Christensen, P. Chiais, P. Clark, M. Collins, J. Comiso, V. de Menezes, R. Feely,

T. Fichefet, A. Fiore, G. Flato, J. Fuglestvedt, G. Hegerl, P. Hezel, G. Johnson, G. Kases, V. Kattsov, J. Kennedy,

A. Tank, C. OL. Quéeré, G. Myhre, T. Osborn, A. Payne, J. Perlwitz, S. Power, M. Prather, S. Rintoul, J. Rogelj,

M. Rusticucci, M. Schulz, J. Secláček, P. Stott, R. Sutton, P. Thorne, and D. Wuebbles, Climate change 2013: The

physical science bases. Contribution of working group I to the fifth assessment report of the intergovernmental panel

on climate change, Cambridge University Press, 2013.

[2] R.A. Armstrong and R. McGehee, Competitive exclusion, Am Nat 115 (1980), no. 2, 151–170.

[3] M.S. Boyce, A.R.E. Sinclair, and G.C. White, Seasonal compensation of predation and harvesting, Oikos 87 (1999),

419–426.

[4] P.L. Chesson, The stabilizing effect of a random environment, J Math Biol 15 (1982), 1–36.

[5] P.L. Chesson and S. Ellner, Invasibility and stochastic boundedness in monotonic competition models, J Math Biol

27 (1989), 117–138.

[6] C. Chicone, Ordinary differential equations with applications, 2nd edition, Texts in Applied Mathematics, Springer,

2006.

[7] P. De Mottoni and A. Schiaffino, Competition systems with periodic coefficients: a geometric approach, J Math Biol

11 (1981), no. 3, 319–335.

[8] A. Erbach, F. Lutscher, and G. Seo, Bistability and limit cycles in generalist predator–prey dynamics, Ecol Com-

plexity 14 (2013), 48–55.

[9] J.M. Fryxell and P. Lundberg, Diet choice and predator–prey dynamics, Evol Ecol 8 (1994), 407–421.

[10] J. Guckenheimer and P. Holmes, Nonlinear oscillations, dynamical systems, and bifurcations of vector fields, Applied

Mathematical Sciences, Springer, 1983.

[11] A. Hastings and T. Powell, Chaos in a three-species food chain, Ecology 72 (1991), no. 3, 896–903.

[12] C.S. Holling, The components of predation as revealed by a study of small-mammal predation of the European pine

sawfly, Canad Entomol 91 (1959), 293–320.

[13] , Some characteristics of simple types of predation and parasitism, Canadian Entomol 91 (1959), 385–398.

[14] R.D. Holt, Predation, apparent competition, and the structure of prey communities, Theor Pop Biol 12 (1977), no. 2,

197–229.

[15] S.-B. Hsu and X.-Q. Zhao, A Lotka–Volterra competition model with seasonal succession, J Math Biol 64 (2012),

109–30.

[16] M. Kot, Elements of mathematical ecology, Cambridge University Press, 2001.

[17] C.J. Krebs, S. Boutin, and R. Boonstra, Ecosystem dynamics of the boreal forest: The Kluane project, Oxford

University Press, 2001.

[18] S.A. Levin, Community equilibria and stability, and an extension of the competitive exclusion principle., Am Nat

104 (1970), 413–423.

[19] C. Magal, C. Cosner, S. Ruan, and J. Casas, Control of invasive hosts by generalist parasitoids, Math Med Biol 25

(2008), 1–20.

[20] T. Namba, Competitive co-existence in a seasonally fluctuating environment, J Theor Biol 111 (1984), no. 2, 369–386.

[21] K Norrdahl and E Korpimäki, Seasonal changes in the numerical responses of predators to cyclic vole populations,

Ecography 25 (2002), 428–438.



122 N. BOLOHAN, V. LEBLANC, AND F. LUTSCHER

[22] L.A. Real, The kinetics of functional response, Am Nat 111 (1977), 289–300.

[23] M.L. Rosenzweig and R.H. MacArthur, Graphical representation and stability conditions of predator–prey interac-

tions, Am Nat 97 (1963), 209–223.

[24] Y. Samia and F. Lutscher, Coexistence and spread of competitors in heterogeneous landscapes, Bull Math Biol 72

(2010), 2089–2112.

[25] H.L. Smith, The interaction of steady state and Hopf bifurcations in a two-predator–one-prey competition model,

SIAM J Appl Math 42 (1982), no. 1, 27–43.

[26] A. Szilágyi and G. Meszéna, Coexistence in a fluctuating environment by the effect of relative nonlinearity: A

minimal model, J Theor Biol 267 (2010), 502–512.

[27] M. Turelli, A reexamination of stability in randomly varying versus deterministic environments with comments on

the stochastic theory of limiting similarity, Theor Pop Biol 13 (1978), 244–267.

[28] , Niche overlap and invasion of competitors in random environments I: Models without demographic stochas-

ticity, Theor Pop Biol 20 (1981), 1–56.

[29] R. Tyson and F. Lutscher, Seasonally varying predation behaviour and climate shifts are predicted to affect predator–

prey cycles, Am Nat 188 (2016), 539–553.

Department of Mathematics and Statistics, University of Ottawa, Ottawa, ON, K1N 6N5, Canada

E-mail address: nbolo094@uottawa.ca

Department of Mathematics and Statistics, University of Ottawa, Ottawa, ON, K1N 6N5, Canada

E-mail address: vleblanc@uottawa.ca

Corresponding author, Department of Mathematics and Statistics, and Department of Biology, Univer-

sity of Ottawa, Ottawa, ON, K1N 6N5, Canada

E-mail address: flutsche@uottawa.ca