Science, Irving K. Barber Faculty of (Okanagan)
Computer Science, Mathematics, Physics and Statistics, Department of (Okanagan)
DSpace
UBCO
Hewa Thudallage, Sanuri
Master of Science - MSc
University of British Columbia
Bylot Island, in the Canadian Arctic, hosts a predator-prey system that includes snowy owls, arctic foxes, and lemmings. The latter are prey for both the snowy owls and foxes. Snowy owls reside temporarily on the island, their presence or absence depending on the availability of the prey only. Hence, unlike many other predator-prey systems, the composition of the food web (one or two predators) depends on the internal system dynamics. We begin with an ordinary differential equation model, that we reduce to a simple onespecies discrete map. This simplification is possible because we assume that the fox and owl population on Bylot Island remains constant because of a fixed number of territories and nests. The simplified model allows us to prove several key aspects of solution behaviour. We investigate the asymptotic behaviour of solutions. We also fully characterise the cyclic and pseudo-cyclic dynamics in the lemming-arctic fox-snowy owl (LFO) system. Next, since arctic regions exhibit drastic seasonal changes, we extend the model to account for some behavioural changes. Specifically, sufficient lemmings are also a necessity for arctic foxes to breed. As in the previous case, we now obtain a model with two thresholds. The second threshold represents the lemming density above which arctic foxes breed. We observe and characterise more complex behaviour in pseudo-cyclic dynamics relative to our first model. We also discover bistable behaviour: long-term dynamics vary between two- or three-cycle solutions depending on the initial conditions. The availability of lemmings determines the presence or absence of snowy owls and the breeding behaviour of arctic foxes in our model. We implement these seasonal changes as thresholds. Using these threshold values and lemming growth rates, we explain the necessary conditions to observe cyclic dynamics in this predator prey system. Moreover, we demonstrate that our simple mathematical model can exhibit complex behaviour that nonetheless contains defined structure. Specifically, we present upper bounds on the number of consecutive increases or decreases in a positive solution of the LFO model based on parameter values. This qualitative behaviour differs from both overcompensation cycles and consumer-resource cycles.
Predator-Prey Interactions in ArcticCanada: A Discrete-Time Model ofLemmings, Foxes and Snowy Owlson Bylot Island.bySanuri Hewa ThudallageA THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF SCIENCEinThe College of Graduate Studies(Mathematics)THE UNIVERSITY OF BRITISH COLUMBIA(Okanagan)September 2024© Sanuri Hewa Thudallage, 2024The following individuals certify that they have read, and recommend to the College ofGraduate Studies for acceptance, a thesis entitled:Predator-Prey Interactions in Arctic Canada: A Discrete-Time Model ofLemmings, Foxes and Snowy Owls on Bylot Island.submitted by Sanuri Hewa Thudallage in partial fulfilment of the requirements of thedegree of Master of ScienceDr. Rebecca C. Tyson, Irving K. Barber Faculty of ScienceSupervisorDr. Frithjof Lutscher, Faculty of Science, University of OttawaCo-SupervisorDr. Eric Foxall, Irving K. Barber Faculty of ScienceSupervisory Committee MemberDr. Wayne Broughton, Irving K. Barber Faculty of ScienceSupervisory Committee MemberDr. Sarah P. Otto, Department of ZoologyUniversity ExamineriiAbstractBylot Island, in the Canadian Arctic, hosts a predator-prey system that includes snowyowls, arctic foxes, and lemmings. The latter are prey for both the snowy owls and foxes.Snowy owls reside temporarily on the island, their presence or absence depending on theavailability of the prey only. Hence, unlike many other predator-prey systems, the compo-sition of the food web (one or two predators) depends on the internal system dynamics.We begin with an ordinary differential equation model, that we reduce to a simple one-species discrete map. This simplification is possible because we assume that the fox and owlpopulation on Bylot Island remains constant because of a fixed number of territories andnests. The simplified model allows us to prove several key aspects of solution behaviour.We investigate the asymptotic behaviour of solutions. We also fully characterise the cyclicand pseudo-cyclic dynamics in the lemming-arctic fox-snowy owl (LFO) system.Next, since arctic regions exhibit drastic seasonal changes, we extend the model toaccount for some behavioural changes. Specifically, sufficient lemmings are also a necessityfor arctic foxes to breed. As in the previous case, we now obtain a model with two thresholds.The second threshold represents the lemming density above which arctic foxes breed. Weobserve and characterise more complex behaviour in pseudo-cyclic dynamics relative to ourfirst model. We also discover bistable behaviour: long-term dynamics vary between two- orthree-cycle solutions depending on the initial conditions.The availability of lemmings determines the presence or absence of snowy owls and thebreeding behaviour of arctic foxes in our model. We implement these seasonal changesas thresholds. Using these threshold values and lemming growth rates, we explain thenecessary conditions to observe cyclic dynamics in this predator prey system. Moreover,we demonstrate that our simple mathematical model can exhibit complex behaviour thatnonetheless contains defined structure. Specifically, we present upper bounds on the numberof consecutive increases or decreases in a positive solution of the LFO model based onparameter values. This qualitative behaviour differs from both overcompensation cyclesand consumer-resource cycles.iiiLay SummaryWe study the lemming-arctic fox-snowy owl (LFO) system on Bylot Island in northernCanada. This predator-prey system involves lemming as the common prey for both thearctic fox and snowy owl. This food web interests us because the presence of the migratorypredator snowy owl on the island depends on the availability of lemmings. This meansthat the structure of the food web is sensitive to the internal state of the prey population.Our aim is to provide a mathematical understanding of this food web and the differentfluctuation patterns in the lemming population.First, we develop a mathematical model that represents the species interactions. Thismodel includes all the 3-species but assumes a constant fox population due to a fixednumber of territories and a constant owl population due to a fixed number of owl nests.These predators predate lemmings at a constant rate throughout the year. As a result ofthese assumptions, we can reduce the model to a single species equation. This allows us todiscuss the behaviours of this food web by focusing on a single species, the lemming. Inour model, there are certain threshold population values that change the structure of thefood web. For example, if the lemming population has grown and crossed a threshold value,then snowy owls begin nesting on the island, i.e., the lemming population is large enoughto support another predator. With this model, we find that the only possible behaviour ofthe population in this system is cyclic behaviour. We also find the needed mathematicalconditions to observe particular patterns in the population cycles in terms of some factorssuch as population growth rates.Next, we extend this model to include fox reproduction in the food web. Assuming thatthe lemming population density required for arctic foxes to reproduce is lower than thatneeded for the snowy owls to start nesting on the island, there are now two thresholds inthe model. The results of this study demonstrate that a simple single species model canexhibit population cycles and certain repetitive patterns. We also find the maximum and theminimum population values in the long run, as a function of the growth rates and thresholdvalues in the system. Furthermore, our work shows that a population in this system canexhibit two types of cycles depending on the population density at the beginning. Forexample, we observe lemming population increase in a year is followed by a decrease in thefollowing year before the next increase. Or else, two increases in the lemming populationduring two consecutive years are followed by a decrease in the following year before the nextincrease.ivPrefacei This thesis is my own original work, carried out under the supervision of my supervisor,Dr. Rebecca Tyson, and my co-supervisor, Dr. Frithjof Lutscher.ii The author of this thesis owns certain copyright or related rights in it (the Copyright)and she has given The University of British Columbia Okanagan certain rights to usesuch Copyright, including for administrative purposes.iii The author has not used generative artificial intelligence for this work.iv The content of this thesis has not been published elsewhere.vContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vContents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiChapter 1: Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Cyclic Populations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Mathematical Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.3 Bylot Island Predator-Prey System . . . . . . . . . . . . . . . . . . . . . . . 31.4 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4Chapter 2: Derivation of the model equations . . . . . . . . . . . . . . . . . 62.1 Scenario 1: Foxes have a constant effect . . . . . . . . . . . . . . . . . . . . 62.2 Scenario 2: Foxes have a seasonally variable effect . . . . . . . . . . . . . . 82.3 Notation and terminology for discrete dynamical systems . . . . . . . . . . 10Chapter 3: Analysis of the one-threshold model . . . . . . . . . . . . . . . . 113.1 Basic model properties and fixed points . . . . . . . . . . . . . . . . . . . . 113.2 Two-cycles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133.3 Three-cycles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143.4 Cycles of length k . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 173.5 Pseudo-cycles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20Chapter 4: Analysis of the two-threshold model . . . . . . . . . . . . . . . . 224.1 Basic model properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22viCONTENTS4.2 Pseudo-cycles with two thresholds . . . . . . . . . . . . . . . . . . . . . . . 254.3 Illustration of the dynamics with two thresholds . . . . . . . . . . . . . . . 274.3.1 Shape of the cycles . . . . . . . . . . . . . . . . . . . . . . . . . . . . 274.3.2 Coexistence of m and n cycles . . . . . . . . . . . . . . . . . . . . . 28Chapter 5: Discussion and conclusions . . . . . . . . . . . . . . . . . . . . . . 35Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37viiList of FiguresFigure 2.1 Illustration of the model with one threshold. Solid line: low lemmingdensity at the beginning of the summer. Dashed line: high lemmingdensity at the beginning of the summer. . . . . . . . . . . . . . . . . 7Figure 2.2 Illustration of the model with two thresholds. Solid line: very lowlemming density at the beginning of the summer. Dotted line: lowlemming density at the beginning of the summer. Dashed line: highlemming density at the beginning of the summer. . . . . . . . . . . 8Figure 3.1 Illustration of the proof of Lemma 3.6. All solutions of (3.2) convergeto a fixed point. Red dashed lines: threshold. Grey dashed line:identity. Black line: updating function f(yn). Blue line: cobweb.Parameter values: R = 1, r = 0.8 (left plot) and R = 1.5, r = 1(right plot). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13Figure 3.2 Illustration of the second-iterate map in the proof of Lemma 3.7.All solutions of (3.5) converge to a fixed point. Red dashed lines:thresholds. Grey dashed line: identity. Black line: updating functionf2(yn). Blue line: cobweb. Parameter values are R =32 and r =23 . . 15Figure 3.3 Illustration of the proof of Lemma 3.8. All solutions of (3.6) convergeto a fixed point when rR2 = 1 (left plot) and when r2R = 1 (rightplot). Red dashed lines: thresholds. Grey dashed line: identity.Black line: updating function f3(yn). Blue line: cobweb. Parametervalues are r = 0.6 (left plot) and R = 1.5 (right plot). . . . . . . . . 16Figure 3.4 Density vs time plot. Two increases followed by a decrease whenrR2 = 1 (left plot). An increase followed by two decreases whenr2R = 1 (right plot). Dashed lines: thresholds. Blue line: updatingfunction f3(yn). Parameter values are r = 0.8 (left plot) and R = 3(right plot). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17Figure 3.5 Illustration of possible and impossible 5-cycles. (a) Rr4 = 1. (b)R2r3 = 1. (c) Impossible 5-cycle. The cases R4r = 1 and R3r2 = 1are in some sense “inverted” from cases (a) and (b): increases becomedecreases and vice versa. . . . . . . . . . . . . . . . . . . . . . . . . 19Figure 3.6 Illustration of how 2-cycles shift depending on initial conditions.Dashed lines: thresholds. Parameter values are r = 0.25 and R = 4. 19viiiLIST OF FIGURESFigure 4.1 Illustration of Lemma 4.1. All solution sequences to Equation (4.1)remain in (0, L′) when the condition h2 ≤ 1 and 1 ≤ h1 < L′ issatisfied. Red dashed lines: thresholds. Grey dashed line: identity.Black line: updating function f(yn). Blue line: cobweb. Parametervalues are L′ = 1.7, h1 = 1.3, h2 = 0.8 < 1 and h3 = 0.4. . . . . . . . 23Figure 4.2 Illustration of Lemma 4.2. All solution sequences to Equation (4.1)are bounded below by 1 when the condition h2 ≥ 1 and 1 ≥ h3 ≥ 1L′is satisfied. Red dashed lines: thresholds. Grey dashed line: identity.Black line: updating function f(yn). Blue line: cobweb. Parametervalues are L′ = 3, h1 = 4, h2 = 1.5 and h3 = 13 . . . . . . . . . . . . . 24Figure 4.3 Sample cycles illustrating the maximum number of consecutive in-creases (Lemma 4.6). Note that for all of the pseudo-cycles shown,there are sometimes multiple consecutive increase steps, but nevermore than one decrease step. According to (4.4), there cannot bemore than 3 (a), 5 (b), and 17 (c) consecutive increases. This re-quirement is fulfilled by all of the pseudo-cycles shown. Parametervalues are h2 = 2 (a), h2 = 1.4 (b), and h2 = 1.1 (c). Fixed param-eters are h1 = 2.3, h3 = 0.2 and L′ = 3. Black dashed lines indicatethe thresholds. The pink dashed and solid lines give the upper andlower bound of the invariant interval. . . . . . . . . . . . . . . . . . 31Figure 4.4 Sample cycles illustrating the maximum number of consecutive de-creases (Lemma 4.7). Note that for all of the pseudo-cycles shown,there are sometimes multiple consecutive decrease steps, but nevermore than one increase step. According to (4.5), for the cases shownthere cannot be more than 3 (a), 2 (b), and 2 (c) consecutive de-creases. This requirement is fulfilled by all of the pseudo-cyclesshown. Parameter values are h3 = 0.4 (a), h3 = 0.3 (b), and h3 = 0.2(c). Fixed parameters are h1 = 15, h2 = 2 and L′ = 2. Black dashedlines indicate the thresholds. The pink dashed and solid lines givethe upper and lower bound of the invariant interval. . . . . . . . . . 32Figure 4.5 Illustration of coexisting cycles. (a) Solutions form a 2-cycle froman initial value at 0.4. (b) Solutions form a 3-cycle from an initialvalue at 0.7. (c) Asymptotic outcome depending on initial condition.Parameters are h1 = 4, h2 = 2, h3 = 0.25 and L′ = 3. Black dashedlines indicate the thresholds. Pink dashed and solid lines correspondto the upper and lower bound of the invariant interval. . . . . . . . 33Figure 4.6 Illustration of the calculations to determine the initial values that giverise to 2-cycles for the following parameter values: h1 = 4, h2 = 2,h3 = 0.25 and L′ = 3. Blue dots: Initial values that give rise to2-cycles. Red dots: Initial values that give rise to 3-cycles. . . . . . 34ixAcknowledgementsMy deepest acknowledgment goes to my supervisors, who are the main pillars of thisgraduate research journey. Dr. Frithjof Lutscher, I am grateful for your guidance and forproviding me with your inspiring knowledge in the field. I am also thankful to Dr. RebeccaTyson for giving me this opportunity to join her research group and, especially for enhancingmy writing skills throughout this project.I extend my sincere thanks to my research group members and committee members,Dr. Eric Foxall, for your continuous encouragement. Additionally, I thank Dr. WayneBroughton for his recommendation on this thesis. This achievement could not be completedwithout both of you.I appreciate my husband, Dinura, for providing much needed emotional support through-out this journey. Lastly, a big thanks to my colleagues for cheering me up always.xTo future myself!For every experience you gathered and for the strength to not to give up.xiChapter 1Introduction1.1 Cyclic PopulationsExploring cyclic behaviour in populations is a study that interests ecologists (Myers,2018). There exist many populations that exhibit cycles over multi-annual periods (Smithet al., 2007; Murdoch et al., 2002). For example, the vole populations in North Amer-ica show both annual and multi-annual cyclicity (Hanski et al., 1991). In northern Europe,there are small rodent populations that exhibit multi-annual cycles (Hanski and Korpima¨ki,1995). In arctic Canada, the snowshoe hare and Canada lynx display multi-annual cyclic-ity (Lindstro¨m et al., 2001), as well as arctic foxes (Elton, 1924), coyotes, and great hornedowls (O’Donoghue et al., 1997). The periodicity of these cyclic behaviours varies from smallvalues like two and three years (Lindstro¨m et al., 2001) to large values like 32 years (El-ton, 1924). Generally, the rodent species exhibit 3–4 year cycles (Ho¨rnfeldt, 1978). Forinstance, voles have a 4-year cycle (Fauteux et al., 2015), and lemmings have a 3–4 yearcycle (Gauthier et al., 2011). Comparatively, the snowshoe hare and many of its predatorshave a lengthy cycle period of 10 years (Adams, 1959). Cycle amplitude may vary fromplace to place (Angerbjo¨rn et al., 2001; Adams, 1959) and species to species (Ims et al.,2008; Stenseth, 1999). Many of the high-amplitude cycles belong to rodents living underextreme weather conditions in the high arctic (Hansson and Henttonen, 1988; Gruyer et al.,2009).For many ecological reasons, it is important to understand these cyclic dynamics. Cyclic-ity in a key resource or prey population can define the structural behaviour of the entireecosystem (Krebs et al., 2001; Ims and Fuglei, 2005). Indeed, collapse of the cycles is a seri-ous issue that affects the food web dynamics in some populations (Ims and Fuglei, 2005; Imset al., 2008). A collapse of the population cycle blocks the periodic delivery of resources tothe ecosystem, potentially leading to extinctions. This collapse is clearly visible in northernconsumer-resource systems (Ims et al., 2008).There are several internal and external forces identified by ecologists that drive a certainsystem to exhibit cyclic behaviour (Barraquand et al., 2017; Myers, 2018; Elton, 1924;Gauthier et al., 2011; Hanski et al., 2001). Factors causing snowshoe hare cycles mayinclude food availability, parasites, diseases, predation, social behaviour, or combinationsof these (Adams, 1959; Keith and Rongstad, 1984; Krebs et al., 2017). Similar mechanismsaffect northern lemming and vole populations (Klemola et al., 1997), though disease maybe the more important driver of the vole cycle (Hanski et al., 2001). The Moran effect isalso considered to be a cause for synchronized population variations (Abbott, 2007). This11.2. MATHEMATICAL MODELSeffect describes how correlated environmental factors such as temperature or precipitationpatterns can synchronize population variations across large regions. Climate change andfragmentation of habitats can disturb these cyclic behaviours.(Ims et al., 2008; Strohm andTyson, 2009; Gauthier, 2004). When all these factors coexist, it is challenging to decidewhich factor influences the cyclicity.1.2 Mathematical ModelsCyclic populations of interacting species have also interested mathematicians over manydecades (Turchin, 2013), and mathematical studies have helped to understand the causes ofcyclic behaviour (Dostalkova et al., 2002). Here we briefly review 2- and 3-species predator-prey models.We begin by discussing some classical mathematical models with two interacting speciesin an ecosystem. The Lotka-Volterra model is a foundational tool in this area (Kot, 2001).This model has a linear functional response between the prey density and prey consump-tion by the predator (Murray, 1993; Kot, 2001), that is, the predator-prey interaction rateis proportional to the product of the population sizes in the system. Solutions of theLotka-Volterra model are periodic, and the model is structurally unstable to small pertur-bations (Murray, 2002; Kot, 2001). An extension of the Lotka-Volterra model is the classicalRosenzweig-MacArthur model (Rosenzweig and MacArthur, 1963; Kot, 2001). This modelincludes a prey growth rate that is regulated by its own population, that is, there is den-sity dependence in prey growth (Turchin, 2013; Murray, 1993). In addition, instead of alinear functional response, which implies that the predator’s consumption rate increaseslinearly with prey density, this model allows a type II functional response. This meansthe consumption rate increases with prey density but at a decreasing rate. This combi-nation of functional responses controls the population density, and the system dynamicsexhibit stable limit cycles under certain parameter values (Murray, 1993). There are manyother well-known mathematical models derived from the Rosenzweig-MacArthur model byvarying the functional response. Some of these include the hyperbolic response (DeAngelismodel, Bazykin model), sigmoid response (Yodzis model, Ludwig model), predator interfer-ence responses, ratio-dependent responses etc. (Turchin, 2013). Thus, there exist various2-species mathematical models depending on the features of the predator-prey interactions.Adding a third species to a 2-species model introduces more complexity to the system,and so there are many different types of three-species predator-prey models (Lin and Kahn,1978; Freedman and Waltman, 1984; Rosenzweig, 1973; Hastings and Powell, 1991; Alebra-heem and Hasan, 2013; Holt and Polis, 1997). One such is the resource-consumer-predatormodel (Rosenzweig, 1973; Lin and Kahn, 1978; Hastings and Powell, 1991; Pattanayak et al.,2021). According to Rosenzweig (1973), this type of food chain has similar interaction out-comes as that of 2-species models. Other authors have found oscillations involving two orall three species (Lin and Kahn, 1978), bistability (Pattanayak et al., 2021), and chaoticbehaviour (Hastings and Powell, 1991). Another common model is that of a three-speciessystem with a predator feeding on two different prey (Freedman and Waltman, 1984; Ter-21.3. BYLOT ISLAND PREDATOR-PREY SYSTEMamoto et al., 1979). In the case where both prey species coexist, the system is balancedby the switching process of the predator between the two available prey depending on theprey abundance. This switching allows all three species to coexist under certain conditions:If the predator consumes the prey with higher intrinsic growth when it becomes abundant,then the prey with lower growth rate can persist (Freedman and Waltman, 1984). Re-searchers have also modeled systems with two predator populations feeding on one prey(Khaliq et al., 2022; Alebraheem and Hasan, 2013; Freedman and Waltman, 1984). Dueto the presence of two predators in the same system, direct or indirect competition amongthese two species may occur (Alebraheem and Hasan, 2013). For instance, higher searchand capture efficiency in one predator can drive the other predator species extinct throughindirect competition. In some cases, one predator may consume the other predator in adirect interaction process called intraguild predation (Vance-Chalcraft et al., 2007; Holtand Polis, 1997; Ji and Wang, 2022). Larger food web models also exist (Hutchison et al.,2020; Farhan et al., 2023; Krebs et al., 2017) and have been useful to explore the connectionbetween biodiversity and higher order trophic levels in an ecosystem (Montoya et al., 2003;Berlow et al., 2004; Tyson et al., 2010).The classic 2- and 3-species model dynamics may vary continuously or discretely if thereare factors that change the species interactions in the ecosystem. Growth rates, survivalrates, and movement rates are all parameters that can vary with season and impact thesystem dynamics (Stenseth, 1999; Hang et al., 2020; Gakkhar et al., 2009; Ims and Fuglei,2005; Gilg et al., 2003). Migration of certain species is another factor that can changefood web dynamics over time (Jin et al., 2019; Bauer and Hoye, 2014; Chen and Zhang,2013; Wang, 2019; Sun et al., 2008). Furthermore, seasonal switching of a predator betweengeneralist and specialist behaviour can lead to new dynamics such as prey extinction, aswell as two types of bistability behaviours where dynamics vary between two stable statesor between a limit cycle and a steady state (Tyson and Lutscher, 2016). Other factors thataffect species interactions or food web structure include, for example, fear of the predatorin prey species (Wang and Zou, 2017; Debnath et al., 2023). Finally, external factors suchas climate, disease, landscape configuration, and human disturbance can also cause changesin an ecosystem (Harvell et al., 2002; Gao et al., 2013; Kumar et al., 2021).1.3 Bylot Island Predator-Prey SystemHere we model a predator prey system on Bylot Island, located next to Baffin Island inthe Nunavut territory of Arctic Canada. A part of Bylot Island is protected as a migratorybird sanctuary, and there exists a study center to keep track of the bird or mammal speciesin these ecosystems (Gauthier, 2004; Bolduc et al., 2023).The predator-prey system that interests us is the lemming, arctic fox, and snowy owl(LFO) system. The only prey in this food web is the lemmimg; both the arctic fox andsnowy owl are lemming predators (Hutchison et al., 2020; Kharitonov et al., 2008; Elton,1924; Lindstro¨m et al., 2001; Gilg et al., 2009; Ims and Fuglei, 2005). The lemming is akey prey species across the arctic (Gilg et al., 2009; Ims and Fuglei, 2005; Gauthier, 2004),31.4. OVERVIEWand so understanding lemming predator-prey systems gives us important information onthe balance of arctic species (Post et al., 2009; Latimer and Zuckerberg, 2019; Norrdahl andKorpima¨ki, 1995). In particular, the high amplitude multi-annual cycles of the lemming areinfluential. Our goal is to use mathematical analysis to study the cyclic dynamics of thelemmings within the LFO system. Below, we describe the biological details needed to buildour model.The abundant lemming plays a key role in the arctic food web on Bylot Island (Gauthier,2004). It is a common prey for arctic predators because of its rapid population growth in theabsence of predators (Hanski et al., 2001; Elton, 1924). During the winter, lemmings buildtunnels under the snow, which protect them from predators. Therefore, lemmings have anintense population growth in the winter season (Fauteux et al., 2015). Next, at the beginningof the summer season, predators begin breeding, and also migratory predators start visitingthe island (Gauthier et al., 2011; Gruyer et al., 2009; Duchesne et al., 2011; Fauteux et al.,2015; Reid et al., 1995). Consequently, lemmings experience high predation rates duringthe summer. In addition to these annual population changes, lemming populations exhibitmulti-annual periodic fluctuations across the arctic (Elton, 1924; Angerbjo¨rn et al., 2001;Stenseth, 1999; Gruyer et al., 2009).Arctic foxes live on Bylot Island throughout the year (Gauthier et al., 2011; Elton,1924). They are territorial, and the number of territories on the island is approximatelyconstant. The foxes will breed only if there exist sufficient lemmings to support reproduction(Samelius and Alisauskas, 2017; Gauthier, 2004; Tannerfeldt et al., 1998). Foxes can leavethe island in the winter when the sea becomes ice (Gravel et al., 2023) which regulates thefox population density on the island.The other predator in our study is the snowy owl, which visits the island in earlysummer in search of a land to begin nesting. This migratory bird also feeds on lemmings(Ims and Fuglei, 2005; Elton, 1924; Gauthier, 2004; Hutchison et al., 2020; Therrien et al.,2014). Because snowy owls prefer nesting in areas with high prey density, their presenceor absence on the island is influenced by the dynamics in the lemming population. Thatmeans, that only if the lemming population is high enough snowy owls will then settle andbreed on the island. This predator exhibits territorial behaviour (Boxall and Lein, 1982).This behaviour avoids population overcrowding. The end of the snowy owl breeding season(and the onset of winter) marks the end of their stay on Bylot Island.1.4 OverviewIn this thesis, we investigate the dynamic changes in the LFO food web with lemmingdensity and how they affect the lemming population cycles. Previous work tends to focuson models with structurally fixed predator-prey interactions, i.e., fixed trophic structure(Ims and Fuglei, 2005; Legagneux et al., 2012; Hastings and Powell, 1991; Turchin, 2013),or models where the parameters (growth, mortality, or predation rates) change seasonally(Tyson and Lutscher, 2016; Gilg et al., 2003; Ims and Fuglei, 2005; Stenseth, 1999; Hanget al., 2020). The LFO system is new because the structure of the predator-prey food41.4. OVERVIEWweb varies depending on the presence of one predator (i.e., the owl) and the trigger forthis change is one or two threshold values in prey abundance, rather than some externalforce such as seasonality. The numerical study in Hutchison et al. (2020) presents the LFOsystem including another prey, the snow goose. They note the need for theory in this area,hence the motivation for our work. In Chapter 3 we present the mathematical analysis ofthe LFO system when the predation pressure by foxes is constant and owls breed if thelemming density is high enough (one-threshold case). In Chapter 4 we introduce a lowerthreshold lemming density that triggers fox breeding, i.e., lower than that which triggersowl breeding (two-threshold case). Finally, in Chapter 5, we discuss our results.5Chapter 2Derivation of the model equationsIn this section, we present the assumptions and derivation of our model equations. Wedivide this section into two parts: in the first part, we consider the impact of the arctic foxon lemmings to be constant throughout the year. Only snowy owls respond to the densityof lemmings. In the second part, we include the response of the arctic fox to the density oflemmings as well.2.1 Scenario 1: Foxes have a constant effectWe make the following assumptions about the growth rate of lemmings and the inter-actions of foxes and owls with lemmings.1. The lemming population shows exponential growth when there are no predators.2. The number of arctic foxes is constant due to a fixed number of available territories.3. Foxes consume lemmings at a constant per capita rate.4. Snowy owls nest on the island only if the density of lemmings exceeds a certainthreshold at the beginning of the summer season. If that is the case, the number ofowls is constant due to a fixed number of nesting sites.5. Snowy owls, if present, consume lemmings at a constant per capita rate.6. Predation by foxes alone is insufficient to prevent growth of the lemming population.Predation by foxes and owls combined does prevent it.To formulate our model based on these assumptions, we denote by t ∈ [0, 1] the timewithin a year and let n = 0, 1, 2, ... denote years. Time t = 0 corresponds to the beginningof the summer season (see Figure 2.1). We denote the density of lemmings in year n attime t by Ln(t). We write F for the number of foxes and aF for their (linear) attack rate onlemmings. Similarly, we use S for the number of snowy owls and aS for their attack rate.Finally, we write L∗ for the threshold density when snowy owls start nesting on the island.We have two equations. When the density of lemmings is too low for owls, we writeL˙n = bLn(t)− aFFLn(t), t ∈ [0, 1], Ln(0) < L∗. (2.1)62.1. SCENARIO 1: FOXES HAVE A CONSTANT EFFECTFigure 2.1: Illustration of the model with one threshold. Solid line: low lemming density atthe beginning of the summer. Dashed line: high lemming density at the beginning of thesummer.When the density of lemmings is sufficient, we haveL˙n = bLn(t)− aFFLn(t)− aSSLn(t), t ∈ [0, 1], Ln(0) ≥ L∗. (2.2)Equations (2.1) and (2.2) are linear and therefore can be solved explicitly from the beginningof one year to the next. The general solution of (2.2) isLn(t) = Ln(0) exp{(b− aFF − aSS)t}, (2.3)and similarly for (2.1). The resulting year-to-year map is given by setting t = 1 so thatLn+1 ={RLn, Ln(0) < L∗,rLn, Ln(0) ≥ L∗,(2.4)where R = eb−aFF and r = eb−aFF−aSS .By our model assumptions and in particular assumption 6, we will generally assumethat0 < r < 1 < R. (2.5)72.2. SCENARIO 2: FOXES HAVE A SEASONALLY VARIABLE EFFECT2.2 Scenario 2: Foxes have a seasonally variable effectFox predation of lemmings varies depending on season (see Figure 2.2). During thewinter, lemmings in their snow galleries are partially protected from fox attacks. Duringthe summer, foxes can choose to reproduce (if the lemming density is high enough) or not(Lai et al., 2016; Goltsman et al., 2005). We state our assumptions as follows.1. During the winter, arctic foxes are the only predators of lemmings. Their attack rateis low because snow partially protects lemmings.2. If the lemming density at the beginning of summer is too low, foxes will not breedso that their energetic needs (and attack rate) will be low. If the lemming density atthe beginning of summer is high enough, foxes will breed so that their energetic needs(and attack rate) will be higher.3. The snowy owls behave as in Scenario 1 above. The threshold for the owls to breedis higher than for the foxes.Figure 2.2: Illustration of the model with two thresholds. Solid line: very low lemmingdensity at the beginning of the summer. Dotted line: low lemming density at the beginningof the summer. Dashed line: high lemming density at the beginning of the summer.82.2. SCENARIO 2: FOXES HAVE A SEASONALLY VARIABLE EFFECTTo formulate our model based on these assumptions, we denote by q the fraction ofthe year that is winter. As above, time within a year is denoted by t ∈ [0, 1] with t = 0corresponding to the beginning of summer. Hence, t ∈ [0, 1 − q) corresponds to summerand t ∈ [1 − q, 1] corresponds to winter. We modify the notation above by introducingsuperscripts s and w for “summer” and “winter,” respectively. We indicate the thresholdsfor arctic foxes and snowy owls to breed by L∗F and L∗S , respectively. Then we have thefollowing equations for the lemming density in year n:dLndt=(rs − a˜sFF )Ln, 0 ≤ t < 1− q, Ln(0) < L∗F ,(rs − asFF )Ln, 0 ≤ t < 1− q, L∗F ≤ Ln(0) < L∗S ,(rs − asFF − asSS)Ln, 0 ≤ t < 1− q, Ln(0) ≥ L∗S ,(rw − awFF )Ln, 1− q ≤ t < 1,(2.6)where a˜sF denotes the attack rate of foxes in a year when they are not breeding. We havea˜sF < asF . We can solve these linear equations as before. For 0 ≤ t < 1 − q, we have threecases:1. If Ln(0) < L∗F :Ln(t) = Ln(0) exp{(rs − a˜sFF )t}. (2.7)2. If L∗F ≤ Ln(0) < L∗S :Ln(t) = Ln(0) exp{(rs − asFF )t}. (2.8)3. If Ln(0) ≥ L∗S :Ln(t) = Ln(0) exp{(rs − asFF − asSS)t}. (2.9)For 1− q ≤ t < 1, all cases are identical, namelyLn(t) = Ln(0) exp{(rw − awFF )(t− (1− q))}. (2.10)We reduce this system to the following year-to-year map:Ln+1 =h1Ln, Ln < L∗F ,h2Ln, L∗F ≤ Ln < L∗S ,h3Ln, Ln ≥ L∗S .(2.11)whereh1 = exp((rs − a˜sFF )(1− q)) exp((rw − awFF )q), (2.12a)h2 = exp((rs − asFF )(1− q)) exp((rw − awFF )q), (2.12b)h3 = exp((rs − asFF − asSS)(1− q)) exp((rw − awFF )q). (2.12c)We have reduced our two impulsive models to difference-equation models (2.4) and(2.11). Before we analyze the qualitative behaviour of these two models, we recall somebasic definitions and results from discrete dynamical systems.92.3. NOTATION AND TERMINOLOGY FOR DISCRETE DYNAMICAL SYSTEMS2.3 Notation and terminology for discrete dynamicalsystemsIn general, a discrete dynamical system is a recursion of the formxn+1 = f(xn), (2.13)usually equipped with an initial condition x0. We recall some basic definitions and resultsthat can be found in most elementary textbooks on discrete dynamical systems, for exampleKot (2001); Allen (2007).1. We call the function f the updating function.2. A state x∗ is called a fixed point if f(x∗) = x∗, so that the constant sequence xn = x∗for all n is a solution of the recursion in (2.13).3. A fixed point is called locally asymptotically stable (l.a.s.) if all nearby solutionsconverge to the fixed point. If at least one solution that starts arbitrarily close to afixed point but does not remain close to the fixed point, then the point is unstable.4. A fixed point is l.a.s. if |f ′(x∗)| < 1. If |f ′(x∗)| > 1, the fixed point is unstable.5. A sequence of p distinct points x¯1, . . . , x¯p is called a p-cycle if f(x¯j) = x¯j+1 forj = 1, . . . , p− 1 and f(x¯p) = x¯1.6. A p-cycle of (2.13) is a fixed point of the p-th iterate map fp = f ◦ · · · ◦ f .7. A p-cycle is stable if the corresponding fixed point of the p-th iterate map is stable.The updating functions of our two models (2.4) and (2.11) is piecewise linear but notcontinuous. Most techniques for analyzing the qualitative behaviour of discrete dynamicalsystems assume that the updating function is smooth (for example the above stabilitycriterion). Hence, while the definitions of fixed points and p-cycles still apply to our models,the usual techniques do not. We will find other techniques in our proofs.10Chapter 3Analysis of the one-thresholdmodelWe study the dynamical behaviour of the model in Equation (2.4). After applying thescaling Ln = Lˆyn, this system becomesyn+1 ={Ryn, yn <L∗Lˆ,ryn,L∗Lˆ≤ yn.(3.1)We choose Lˆ = L∗ and obtain the nondimensionalized systemyn+1 = f(yn) ={Ryn, yn < 1,ryn, 1 ≤ yn.(3.2)We will determine the asymptotic behaviour of solutions of this model, which depend onthe value of the parameters R and r. We begin with a general observation that we will userepeatedly throughout our analysis. The solution to the difference equationzn+1 = Rzn (3.3)is given by the geometric sequencezn = Rnz0. (3.4)In particular, for a positive initial value z0, we find that if R > 1 then zn → ∞, while if0 < R < 1 then zn → 0.3.1 Basic model properties and fixed pointsWe begin with some general considerations about model (3.2) without imposing theassumption that 0 < R < 1 < r, but only the weaker assumption that R, r > 0, which iswarranted by the model derivation. It is obvious that the zero state is a fixed point of themodel. Our first observation is as follows.Lemma 3.1. Solutions of (3.2) with nonnegative initial condition remain nonnegative.Proof. Since R, r > 0, we have f(y) ≥ 0 whenever y ≥ 0. Since yn+1 = f(yn), it is clearthat solutions remain nonnegative.113.1. BASIC MODEL PROPERTIES AND FIXED POINTSLemma 3.2. If 0 ≤ r,R < 1, then all nonnegative solutions converge to zero, i.e., zero isglobally asymptotically stable.Proof. Let r∗ = max{r,R} < 1. Then the solution of zn+1 = r∗zn converges to zero forall initial conditions; see (3.4). Furthermore, since yn+1 = f(yn) ≤ r∗yn, the solution yn isbounded above by the solution zn. Hence, it also has to converge to zero.Lemma 3.3. If r,R > 1, then all nonnegative solutions grow without boundsProof. Let r∗ = min{r,R} > 1. Then the solution of zn+1 = r∗zn grows without bound forall non-negative initial conditions; see (3.4). Furthermore, since yn+1 = f(yn) ≥ r∗yn, thesolution yn is bounded below by the non-negative solution zn. Therefore, all non-negativesolutions yn grow without bound.Lemma 3.4. If 0 ≤ R < 1 < r, then all solutions with 0 ≤ y0 < 1 converge to zero whileall solutions with y0 > 1 grow without boundProof. Let y0 < 1. Then y1 = Ry0 < y0 < 1 since R < 1. Therefore, solutions with initialconditions between 0 and 1 will remain between 0 and 1. Hence, they satisfy yn+1 = Rynwith R < 1 so that yn → 0 according to (3.4). The steps in the case for y0 > 1 aresimilar.Lemma 3.5. If 0 ≤ r < 1 < R, then zero is an unstable fixed point and all positive solutionsremain bounded.Proof. Let 0 ≤ yn < 1. Then yn+1 = Ryn > yn since R > 1 and therefore zero is unstable.For boundedness of the solution, consider 0 < yn < 1. Then yn+1 = Ryn < R. If yn ≥ 1,then yn+1 = ryn < yn since r < 1. Hence, a solution sequence will decrease while its valueexceeds unity, but it cannot increase above R once its value is below unity. Hence, for anyinitial condition, y0, the solution sequence is bounded by max{y0, R}.Lemma 3.6. System (3.2) has positive fixed points if and only if R = 1 or r = 1. In eithercase, there is a continuum of fixed points. If R = 1 and r < 1, all nonnegative solutionsconverge to a fixed point. Similarly, if R > 1 and r = 1, all nonnegative solutions convergeto a fixed point.Proof. Let y∗ be a positive fixed point for System (3.2). Then y∗ = Ry∗ with 0 < y∗ < 1or y∗ = ry∗ with y∗ > 1, i.e., R = 1 or r = 1. Now let R = 1 or r = 1. If R = 1, thenyn+1 = yn for all 0 < yn < 1. Therefore, yn = y∗ is a positive fixed point. If r = 1, thenyn+1 = yn for all yn > 1. Hence, yn = y∗ is a positive fixed point. Thus, the system hasfixed points if only if R = 1 or r = 1.If R = 1, then each y ∈ (0, 1) is a fixed point. If r = 1, then each y ∈ (1,∞) is a fixedpoint. Therefore, in either case, there is a continuum of fixed points.Next, let R = 1 and r < 1 (see Figure 3.1(a)). Let {yn} be a nonnegative solution. If0 < y0 < 1 then y0 is a fixed point (see above) and hence the sequence converges to a fixedpoint. If y0 > 1, then yn+1 = ryn as long as yn ≥ 1. Hence, by (3.4), the sequence will123.2. TWO-CYCLEScontinue to decrease until yk < 1 for some k ≥ 1. But then yk+1 = yk. Hence, the solutionconverges to a fixed point.Similarly, let R > 1 and r = 1 (see Figure 3.1(b)). Let {yn} be a nonnegative solution.If y0 > 1 then y0 is a fixed point (see above) and hence the sequence converges to a fixedpoint. If y0 < 1, then yn+1 = Ryn as long as yn < 1. Hence, by (3.4), the sequence willcontinue to increase until yk > 1 for some k ≥ 1. But then yk+1 = yk. Hence, the solutionconverges to a fixed point.(a) (b)Figure 3.1: Illustration of the proof of Lemma 3.6. All solutions of (3.2) converge to afixed point. Red dashed lines: threshold. Grey dashed line: identity. Black line: updatingfunction f(yn). Blue line: cobweb. Parameter values: R = 1, r = 0.8 (left plot) andR = 1.5, r = 1 (right plot).3.2 Two-cyclesGiven the basic properties of the model as well as the biological interpretation, we willhenceforth assume that 0 < r < 1 < R. Hence, all solutions remain bounded, the zero state133.3. THREE-CYCLESis unstable, and there is no positive fixed point. Here, we study two-cycles, which are fixedpoints of the second-iterate map of our system.Lemma 3.7. System (3.2) exhibits two-cycles if and only if rR = 1. If there are two-cycles,then there is a continuum of them, and all positive solutions approach a two-cycle.Proof. The second-iterate map of System (3.2) is given byyn+2 =R2yn if yn <1R ,rRyn if1R ≤ yn < 1r ,r2yn if yn ≥ 1r .(3.5)Studying two-cycles in System (3.2) is equivalent to studying fixed points in System (3.5).Let y∗ be a positive fixed point for System (3.5). Then y∗ = R2y∗ with 0 < y∗ < 1R ory∗ = rRy∗ with 1R < y∗ < 1r or y∗ = r2y∗ with 1r < y∗, i.e., R2 = 1 or rR = 1 or r2 = 1.Only rR = 1 is possible since 0 < r < 1 < R. Now let rR = 1. Then yn+1 = yn for all1R ≤ yn < 1r . Hence, yn = y∗ is a positive fixed point in System (3.5). Thus, System (3.5)has fixed points if only if rR = 1.Now assume that rR = 1. Then each y ∈ [ 1R , 1r ) is a fixed point of (3.5). Therefore, ifSystem (3.2) has two-cycles, then there is a continuum of them.Next, let 1R ≤ y0 < 1r . Then y0 is a fixed point and hence the corresponding solution{yn} converges to a fixed point. If y0 < 1R , then y2 = R2y0 > y0 since R > 1. As in (3.4),the sequence generated by (3.5) will increase geometrically until yk−2 < 1R and yk >1R forsome k ≥ 2. We claim that yk < 1r . Since yk−2 < 1R , we find yk = R2yk−2 < R2 1R = R = 1r .But that means that yk is a fixed point. Therefore, the solution of System (3.5) convergesto a fixed point (see Figure 3.2). The argument for y0 >1r is similar. Hence, all solutionsof System (3.5) converge to a fixed point.3.3 Three-cyclesThe considerations about two-cycles from the preceding section can be applied to studythree-cycles, which are fixed points of the third-iterate map of our system.Lemma 3.8. The system exhibits three-cycles if and only if rR2 = 1 or r2R = 1. If thereare three-cycles, then there is a continuum of them, and all solutions approach a three-cycle.Proof. The third-iterate map for System (3.2) is given byyn+3 =R3yn if yn <1R2,rR2yn if1R2≤ yn < 1rR ,Rr2yn if1Rr ≤ yn < 1r2 ,r3yn if yn ≥ 1r2 .(3.6)143.3. THREE-CYCLESFigure 3.2: Illustration of the second-iterate map in the proof of Lemma 3.7. All solutionsof (3.5) converge to a fixed point. Red dashed lines: thresholds. Grey dashed line: identity.Black line: updating function f2(yn). Blue line: cobweb. Parameter values are R =32 andr = 23 .Exploring three-cycles in System (3.2) is equivalent to studying fixed points in Sys-tem (3.6).Let y∗ be a positive fixed point for System (3.6). Then y∗ = R3y∗ with 0 < y∗ < 1R2or y∗ = rR2y∗ with 1R2≤ y∗ < 1rR or y∗ = Rr2y∗ with 1Rr ≤ y∗ < 1r2 or y∗ = r3y∗ with1r2≤ y∗, i.e., R3 = 1 or rR2 = 1 or Rr2 = 1 or r3 = 1. Only rR2 = 1 or Rr2 = 1 are possiblesince 0 < r < 1 < R. Now let rR2 = 1 or r2R = 1. Then yn+3 = yn when1R2≤ yn < 1rRor 1rR ≤ yn < 1r2 . Hence, yn = y∗ is a positive fixed point in System (3.6). Thus, theSystem (3.6) has fixed points if only if rR2 = 1 or r2R = 1.If rR2 = 1 or r2R = 1, then each y ∈ [ 1R2, 1rR) or y ∈ [ 1rR , 1r2 ) is a fixed point of (3.6).Therefore, if System (3.2) has three-cycles, then there is a continuum of them.Next, let rR2 = 1 and 1R2≤ y0 < 1rR . Then y0 is a fixed point and hence the correspond-ing solution y0 converges to a fixed point. If y0 <1R2, then y3 = R3y0 > y0 since R > 1.As in (3.4), the sequence will increase geometrically until yk−3 < 1R2 and yk ≥ 1R2 for somek ≥ 3. We claim that yk < 1rR . Since yk−3 < 1R2 , we find yk = R3yk−3 < R3 1R2 = R = 1rR .But that means that yk is a fixed point. Therefore, the solution of System (3.6) convergesto a fixed point (see Figure 3.3(a)). The argument for 1rR ≤ y0 < 1r2 and y0 ≥ 1r2 are similar.153.3. THREE-CYCLESOr else, let Rr2 = 1 when 1rR ≤ y0 < 1r2 . Then y0 is a fixed point and hence thecorresponding solution y0 converges to a fixed point. If y0 <1R2, then y3 = R3y0 > y0 sinceR > 1. As in (3.4), the sequence will increase geometrically until yk−3 < 1R2 and yk ≥ 1R2 forsome k ≥ 3. We claim that yk < 1r2 . Since yk−3 < 1R2 , we find yk = R3yk−3 < R3 1R2 = R =1r2. But that means yk is a fixed point. Therefore, the solution of System (3.6) converges toa fixed point (see Figure 3.3 (b)). The argument for 1R2≤ y0 < 1Rr and y0 ≥ 1r2 are similar.Hence, all solutions of System (3.6) converge to a fixed point.(a) (b)Figure 3.3: Illustration of the proof of Lemma 3.8. All solutions of (3.6) converge to afixed point when rR2 = 1 (left plot) and when r2R = 1 (right plot). Red dashed lines:thresholds. Grey dashed line: identity. Black line: updating function f3(yn). Blue line:cobweb. Parameter values are r = 0.6 (left plot) and R = 1.5 (right plot).While there is only one type of two cycle, we observe two different types of three-cycles.When rR2 = 1 we see two increases followed by a decrease (Figure 3.4(a)). When r2R = 1we see an increase followed by two decreases (Figure 3.4(b)). In the next section, we willsee that the above results can be generalized to cycles of arbitrary length, but that thereare constraints in terms of the qualitative shapes of cycles.163.4. CYCLES OF LENGTH k(a) (b)Figure 3.4: Density vs time plot. Two increases followed by a decrease when rR2 = 1 (leftplot). An increase followed by two decreases when r2R = 1 (right plot). Dashed lines:thresholds. Blue line: updating function f3(yn). Parameter values are r = 0.8 (left plot)and R = 3 (right plot).3.4 Cycles of length kLemma 3.9. The system exhibits a k-cycle if and only if rRk−1 = 1 or r2Rk−2 = 1 or . . .rk−1R = 1 for k ≥ 1. If there is a k-cycle, then there is a continuum of them.The proof of this statement follows exactly as the proofs or the preceding lemmas. Infact, the kth-iterate map for System (3.2) can be calculated explicitly asyn+k =Rkyn if yn <1Rk−1 ,rRk−1yn if 1Rk−1 ≤ yn < 1rRk−2 ,r2Rk−2yn if 1rRk−2 ≤ yn < 1r2Rk−3 ,...rk−1Ryn if 1rk−2R ≤ yn < 1rk−1 ,rkyn if yn ≥ 1rk−1 .(3.7)We conjecture that it is also true in general that all solutions converge to a k-cycle ifthe condition rnRk−n = 1 is satisfied, but the explicit calculations are more cumbersome.Instead, we study the qualitative patterns that can and cannot appear in a k-cycle in thelemma below.First, however, we briefly discuss the conditions for the existence of k-cycles in termsof the original parameters of our model. Let us assume that rnRk−n = 1. We recall thatR = eb−aFF and r = eb−aFF−aSS . Substituting these expressions into the conditions for173.4. CYCLES OF LENGTH kk-cycles, we obtain (eb−aFF−aSS)n (eb−aFF)k−n= 1. (3.8)We can solve for the intrinsic growth rate in terms of the attack rates on lemmings andfoxes and find the relationshipb = aFF +nkaSS. (3.9)We now turn to the qualitative study of the cycle patterns in system (3.7). The followinglemma imposes some surprising constraints on the type k-cycles that are possible.Lemma 3.10. If a k-cycle contains two or more consecutive increases from one generationto the next, then it cannot contain more than one consecutive decrease. Likewise, if a k-cycle contains two or more consecutive decreases from one generation to the next, then itcannot contain more than one consecutive increase.Proof. We prove this lemma by contradiction. Assume that there is a k-cycle, denoted by{y1, . . . yk}, with at least two consecutive increases and at least two consecutive decreases.Let ymin denote the minimum and ymax the maximum of the k-cycle. Since there are twoconsecutive increases, there must be a value yn < 1 such that yn+1 = Ryn < 1. In particular,Rymin < 1. Since one can only arrive at the minimum by a decrease and since the lowestvalue from which a decrease occurs is the threshold, we must have r × 1 = r ≤ ymin. ThenRr ≤ Rymin < 1.Since there are two consecutive decreases, there must be a value yn ≥ 1 such thatyn+1 = ryn ≥ 1. In particular, rymax ≥ 1. Since one can only arrive at the maximum byan increase and since the highest value from which an increase occurs is the threshold, wemust have R× 1 = R ≥ ymax. Then rR ≥ rymax ≥ 1.Thus, the presence of two consecutive increases and two consecutive decreases in a k-cycle gives a contradiction.We discuss and illustrate the implications of this lemma in the cases of k = 4 andk = 5. Just like in the case of 3-cycles (compare Figure 3.4), there can be 4-cycles with3 consecutive increases followed by a single decrease (if R3r = 1) or with 3 consecutivedecreases followed by one increase (if Rr3 = 1). However, there cannot be a 4-cycle withtwo consecutive increases followed by two decreases as might be expected if R2r2 = 1.Instead, since R2r2 = 1 implies Rr = 1, there will be 2-cycles.The case of 5-cycles is more interesting. As in the case of 3- and 4-cycles, we find 5-cycleswith four consecutive increases (decreases) followed by one decrease (increase) if R4r = 1(Rr4 = 1). In the case that R3r2 = 1, we could expect a 5-cycle with three consecutiveincreases followed by two consecutive decreases, but those are excluded by Lemma 3.10.Instead, we find a 5-cycle with two consecutive increases followed by alternating decreasesand increases. The case R2r3 = 1 is similar. We illustrate what is possible and what is notin Figure 3.5.183.4. CYCLES OF LENGTH k(a) (b) (c)Figure 3.5: Illustration of possible and impossible 5-cycles. (a) Rr4 = 1. (b) R2r3 = 1. (c)Impossible 5-cycle. The cases R4r = 1 and R3r2 = 1 are in some sense “inverted” fromcases (a) and (b): increases become decreases and vice versa.We close this section with a short comment about stability of cycles in our model.Cycles occur only under very specific conditions on parameters (Rnrk−n = 1), but if theyoccur, then there exists a continuum of them. Hence, we cannot expect the cycles to beasymptotically stable. In fact, depending on initial conditions, solutions approach cyclesthat are shifted up or down. This is illustrated in the case of 2-cycles with three differentinitial conditions in Figure 3.6. This behaviour is possible because our model is piecewiselinear. By contrast, 2-cycles that appear from a flip bifurcation in classical nonlinear discretemodels such as the logistic or the Ricker model are unique and locally or even globallyasymptotically stable. Solutions from all (nearby) initial conditions converge to the unique2-cycle.(a) (b) (c)Figure 3.6: Illustration of how 2-cycles shift depending on initial conditions. Dashed lines:thresholds. Parameter values are r = 0.25 and R = 4.193.5. PSEUDO-CYCLES3.5 Pseudo-cyclesFinally, we consider the dynamical behaviour of Equation (3.2) without the assumptionthat solutions form cycles. Hence, we shall always assume that 0 < r < 1 < R but we donot necessarily have integers n, k such that Rnrk−n = 1.From our observation (3.4), we know that solutions increase geometrically as long asthey are below the threshold and decrease geometrically above the threshold. Since thereare no fixed points under the assumption of this section, solutions will continue to oscillate:they increase until they first exceed the threshold and then decrease until they fall belowthe threshold again. Only if the condition Rnrk−n = 1 holds for some integers n, k can theseoscillations be true cycles in the sense of Section 2.3. In general, we call such oscillatorysolutions “pseudo-cycles”. The main goal of this section is to show that even for such pseudo-cycles, there are constraints regarding the number of consecutive increases and decreases inthe same spirit as Lemma 3.10.We begin with the simple observation that if Rr < 1, then it follows that there existssome k ∈ N such that Rr < R2r < R3r · · · < Rkr < 1 ≤ Rk+1r < . . . . Similarly, if Rr > 1,then it follows that there exists some l ∈ N such that Rr > Rr2 > Rr3 · · · > Rrl ≥ 1 >Rrl+1 > . . . .Lemma 3.11. Let {yn} be a positive solution of System (3.2). If y0 ≥ 1, then there existsa smallest n¯ ∈ N such that yn¯ ≥ 1 > yn¯+1. Similarly, if y0 < 1, then there exists a smallestn¯ ∈ N such that yn¯ < 1 ≤ yn¯+1. Furthermore, r ≤ yn¯+k < R for any k ≥ 1.Proof. Let y0 ≥ 1. Since r < 1, the sequence rky0 → 0 as k → ∞. Therefore, thereexists exactly one k such that rky0 ≥ 1 and rk+1y0 < 1. We set n¯ = k. Wheneveryn ≥ 1 and yn+1 < 1, then yn+1 = ryn ≥ r. Similarly, let y0 < 1. Since R > 1, thesequence Rk+1y0 →∞ as k →∞. Therefore there exists exactly one k such that Rky0 < 1and Rk+1y0 ≥ 1. In this case, we set n¯ = k. Whenever yn < 1 and yn+1 ≥ 1, thenyn+1 = Ryn < R.To study the asymptotic behaviour of solutions as n → ∞, we may discard transientsfor finite n. Hence, we may assume from hereon that r ≤ yn < R.Theorem 3.12. Suppose rR < 1 and Rr < R2r < R3r · · · < Rkr < 1 ≤ Rk+1r < . . . forsome k ∈ N. Then a positive solution cannot contain two or more consecutive decreasesand not more than k + 1 consecutive increases.Proof. Assume that there exist two consecutive decreases. Then there exists y0 ≥ 1 suchthat y1 = ry0 ≥ 1. We have ry0 < rR since y0 ∈ [r,R). Therefore, we get ry0 < rR < 1.This contradicts y1 = ry0 ≥ 1. Therefore, a positive solution cannot contain two or moreconsecutive decreases.Next, assume that there exist k + 2 consecutive increases. Then there exists y0 < 1 suchthat yk+1 = Rk+1y0 < 1. We have Rk+1y0 ≥ Rk+1r since y0 ∈ [r,R). Therefore, we getRk+1y0 ≥ Rk+1r ≥ 1. This contradicts yk+1 = Rk+1y0 < 1. Therefore, a positive solutioncannot contain more than k + 1 consecutive increases.203.5. PSEUDO-CYCLESTheorem 3.13. Suppose rR ≥ 1 and Rr > Rr2 > Rr3 · · · > Rrl ≥ 1 > Rrl+1 > . . . forsome l ∈ N. Then a positive solution cannot contain two or more consecutive increases andnot more than l + 1 consecutive decreases.Proof. Assume that there exist two consecutive increases. Then there exists y0 < 1 suchthat y1 = Ry0 < 1. We have Ry0 ≥ Rr since y0 ∈ [r,R). Therefore, we get Ry0 ≥ Rr ≥ 1.This contradicts y1 = Ry0 < 1. Therefore, a positive solution cannot contain two or moreconsecutive increases.Next, assume that there exist l + 2 consecutive decreases. Then there exists y0 ≥ 1 suchthat yl+1 = rl+1y0 ≥ 1. We have rl+1y0 < rl+1R since y0 ∈ [r,R). Therefore, we getrl+1y0 < rl+1R < 1. This contradicts yl+1 = rl+1y0 ≥ 1. Therefore, a positive solutioncannot contain more than l + 1 consecutive decreases.In summary, in this chapter, we studied the qualitative behaviour of solutions in theone-threshold model (3.2). We began with some basic properties about positivity andboundedness to ensure that our model satisfies the minimum requirements. Because of themodel assumption 0 < R < 1 < r, our model exhibits oscillatory solutions. In special cases,namely when Rnrj = 1 for positive integers n, j, our model has cycles of length n + j.Because the model is linear, these cycles are not unique, but there exists a continuum ofthem. In particular, they cannot be asymptotically stable for the dynamics. Most impor-tantly, there is a severe and somewhat surprising constraint on the number of consecutiveincreases or decreases in a cycle. When no true cycles exist, we still observe oscillatorysolutions, which we call pseudo-cycles. It turns out that, the constraints on the number ofconsecutive increases or decreases carries over to pseudo-cycles. While this might seem sur-prising at first, we could have speculated that such constraints exist based on the followingobservation: The condition Rnrj = 1 is equivalent to∣∣∣∣ lnRln r∣∣∣∣ = jn.Hence, if a given pair of parameters, R, r, does not allow for true cycles, because the fractionon the left-hand side is not rational, then there are cycles “arbitrarily close” because therational numbers are dense in the real numbers. However, to be arbitrarily close to theexpression on the left-hand side, the fraction j/n might be made up of very large integersj, n, so that the cycle would have very long period j + n. We leave open the question ofwhether there is chaos in our system.21Chapter 4Analysis of the two-thresholdmodelIn this section, we consider the two-threshold model for the lemming population, givenby Equation (2.6). We set the scaled lemming density yn =LnLˆ. Then we choose Lˆ = L∗f ,so that the model takes the form:yn+1 =h1yn, yn < 1,h2yn, 1 ≤ yn < L′,h3yn, yn ≥ L′.(4.1)where L′ = L∗oL∗f> 1.4.1 Basic model propertiesWith regard to the basic properties and biological dynamics, assume that 0 ≤ h3 < h2 <h1 and h3 < 1 < h1 in Equation (4.1). Under these assumptions and as in the precedingsection, zero is an unstable fixed point and all positive solutions remain bounded.Lemma 4.1. If h2 ≤ 1 and 1 < h1 < L′, then all positive solutions of System (4.1)eventually enter the invariant interval (0, L′).Proof. We start by showing that (0, L′) is invariant. Let 0 < ym < L′ for some m > 0. Ifym ∈ (0, 1), then ym+1 = h1ym > ym since h1 > 1. Hence, solutions cannot converge tozero. Furthermore, ym+1 = h1ym < L′ since h1 < L′. If ym ∈ [1, L′) then ym+1 ≤ ym sinceh2 ≤ 1. Hence, the interval is invariant.Next, we show that all solutions eventually enter the invariant interval. We only needto show this for initial conditions outside the interval. Suppose y0 ≥ L′. Then, because0 < h3 < 1, we have ym+1 = h3ym < ym as long as ym > L′. Since the sequence hm3 → 0,there is a smallest m such that hm3 y0 ≥ L′ > hm+13 y0. At that point, the solution hasentered the invariant interval. Figure 4.1 illustrates this behaviour.Lemma 4.2. If h2 ≥ 1 and 1 > h3 ≥ 1L′ , then all positive solutions eventually enter theinvariant interval [1, h2L′).Proof. We show that [1, h2L′) is an invariant interval. Let 1 ≤ ym < h2L′ for some m > 0.If ym ∈ [1, L′), then ym+1 = h2ym ≥ ym ≥ 1 since h2 ≥ 1. Also, ym+1 < h2L′ since ym < L′.224.1. BASIC MODEL PROPERTIESFigure 4.1: Illustration of Lemma 4.1. All solution sequences to Equation (4.1) remainin (0, L′) when the condition h2 ≤ 1 and 1 ≤ h1 < L′ is satisfied. Red dashed lines:thresholds. Grey dashed line: identity. Black line: updating function f(yn). Blue line:cobweb. Parameter values are L′ = 1.7, h1 = 1.3, h2 = 0.8 < 1 and h3 = 0.4.If ym ∈ [L′, h2L′), then ym+1 = h3ym < ym < h2L′ since h3 < 1 and ym+1 ≥ 1 sinceh3 ≥ 1/L′. Hence, the interval is invariant.Next, we show that all solutions eventually enter the invariant interval. We only needto show that for initial conditions outside the interval. Suppose 0 < y0 < 1. Then ym+1 =h1ym > ym as long as ym < 1 because h1 > 1. Since hm1 → ∞, there is a smallest msuch that hm1 y0 < 1 ≤ hm+11 y0. At that point, the solution has either entered the invariantinterval or reached above h2L′.Suppose y0 ≥ h2L′. Then, since h3 < 1, we have ym+1 = h3ym < ym as long asym > h2L′. Since hm3 → 0, there is a smallest m such that hm3 y0 ≥ h2L′ > hm+13 y0.Furthermore, when ym > h2L′, we have ym+1 = h3ym ≥ h3h2L′ ≥ h2 ≥ 1. Hence, thesolution enters the invariant interval. Figure 4.2 illustrates this behaviour.Lemma 4.1 states that when h2 ≤ 1 and h1 < L′, System (4.1) asymptotically reduces234.1. BASIC MODEL PROPERTIESFigure 4.2: Illustration of Lemma 4.2. All solution sequences to Equation (4.1) are boundedbelow by 1 when the condition h2 ≥ 1 and 1 ≥ h3 ≥ 1L′ is satisfied. Red dashed lines:thresholds. Grey dashed line: identity. Black line: updating function f(yn). Blue line:cobweb. Parameter values are L′ = 3, h1 = 4, h2 = 1.5 and h3 = 13 .to the following one-threshold model:yn+1 ={h1yn, yn < 1,h2yn, 1 ≤ yn < L′.(4.2)Similarly, Lemma 4.2 proves that when h2 ≥ 1 and 1 > h3 ≥ 1L′ , System (4.1) asymptoticallyreduces to the one-threshold model:yn+1 ={h2yn, 1 ≤ yn < L′h3yn, L′ ≤ yn < h2L′.(4.3)We studied the qualitative behaviour of models with a single threshold in the precedingchapter. Our interest here is in cycles that do not reduce to the one-threshold case that we244.2. PSEUDO-CYCLES WITH TWO THRESHOLDShave already studied. More specifically, we are interested in cycles that include both thefirst region (below the first threshold) and the third region (above the second threshold).For the remainder of this thesis, we will restrict our attention to cycles of this type.4.2 Pseudo-cycles with two thresholdsIn order for the population cycles to include the first and third population regions, werequire that at least one of the conditions in Lemma 4.1 or Lemma 4.2 above be violated.There are two cases, namely,1. h2 ≤ 1 and h1 ≥ L′,2. h2 ≥ 1 and h3 < 1L′ .We now study the properties of our system under these conditions. We begin by showingthat solutions eventually enter some invariant bounded interval. We consider the two casesseparately.Lemma 4.3. If h2 ≤ 1 and h1 ≥ L′, then all positive solutions of System (4.1) eventuallyenter the invariant interval [min{h2, h3L′}, h1).Proof. We start by showing that [min{h2, h3L′}, h1) is invariant. Let min{h2, h3L′} ≤ ym 0. If ym ∈ [min{h2, h3L′}, 1), then ym+1 = h1ym ≥ ym ≥ min{h2, h3L′}since h1 ≥ L′ > 1. Also, ym+1 = h1ym < h1 since ym < 1. If ym ∈ [1, L′) then ym+1 =h2ym ≤ ym < L′ ≤ h1 since h2 ≤ 1. Also, ym+1 = h2ym ≥ h2 ≥ min{h2, h3L′} since ym ≥ 1.If ym ∈ [L′, h1), then ym+1 = h3ym < ym < h1 since h3 < 1. Also, ym+1 = h2ym ≥ h2 ≥min{h2, h3L′} since ym ≥ h2. Hence, the interval is invariant.Next, we show that all solutions eventually enter the invariant interval. We only needto show this for initial conditions outside the interval. Suppose y0 < min{h2, h3L′}. Thenym+1 = h1ym ≥ ym as long as ym < 1, because h1 ≥ 1. Since the sequence hm1 →∞, thereis a smallest m such that hm1 y0 < 1 ≤ hm+11 y0. Furthermore, ym+1 = h1ym < h1 sinceym < 1. Hence, the solution will enter the invariant interval.Now suppose y0 ≥ h1. Then ym+1 = h3ym < ym as long as ym ≥ L′, because 0 < h3 < 1.Since the sequence hm3 → 0, there is a smallest m such that hm3 y0 ≥ L′ > hm+13 y0. At thatpoint, the solution will have entered the invariant interval since ym ≥ L′ implies thath3ym ≥ h3L′.Lemma 4.4. If h2 ≥ 1 and h3 < 1L′ , then all positive solutions of System (4.1) eventuallyenter the invariant interval [h3L′,max{h1, h2L′}).Proof. We start by showing that [h3L′,max{h1, h2L′}) is invariant. Let h3L′ ≤ ym <max{h1, h2L′} for some m > 0. If ym ∈ [h3L′, 1), then ym+1 = h1ym > ym ≥ h3L′since h1 > 1. Also, ym+1 < h1 ≤ max{h1, h2L′} since ym < 1. If ym ∈ [1, L′) thenym+1 = h2ym ≥ ym ≥ 1 > h3L′ since h2 ≥ 1. Also, ym+1 = h2ym < h2L′ ≤ max{h1, h2L′}254.2. PSEUDO-CYCLES WITH TWO THRESHOLDSsince ym < L′. If ym ∈ [L′,max{h1, h2L′}), then ym+1 = h3ym < ym < max{h1, h2L′} sinceh3 < 1. Also, ym+1 = h3ym ≥ h3L′ since ym ≥ L′. Hence, the interval is invariant.Next, we show that all solutions eventually enter the invariant interval. We only needto show this for initial conditions outside the interval. Suppose y0 < h3L′. Then ym+1 =h1ym > ym as long as ym < 1 because h1 > 1. Since the sequence hm1 → ∞, there is asmallest m such that hm1 y0 < 1 ≤ hm+11 y0. At that point, the solution has entered theinvariant interval since ym < 1 implies that h1ym < h1.Suppose y0 ≥ max{h1, h2L′}. Then ym+1 = h3ym < ym as long as ym ≥ L′ because 0 <h3 < 1. Since the sequence hm3 → 0, there is a smallest m such that hm3 y0 ≥ L′ > hm+13 y0.Furthermore, ym+1 = h3ym ≥ h3L′ since ym ≥ L′. Hence, at that point, the solution hasentered the invariant interval.Of the two possible cases yielding two-threshold cycles (see above), we choose the onewhere h2 ≥ 1 and h3 < 1/L′. Our choice is inspired by the system on Bylot Island whereit is observed that predation by the arctic fox alone is not sufficient to reduce lemmingpopulation growth below unity (hence h2 ≥ 1) while the addition of owl predation doesresult in a net decrease of the lemming population (hence h3 < 1).We now proceed to a more detailed asymptotic analysis of our chosen type of two-threshold cycles. Our aim is to characterize the qualitative behaviour of solutions in a waythat is similar to our approach in Theorems 3.12 and 3.13 for the one-threshold case. Forthe long-term dynamics, we may assume that solutions have entered the invariant interval.By the preceding lemmas, this invariant interval is [h3L′,max{h1, h2L′}) for our parameterchoices. Since our model has two thresholds, we expect to obtain more complex behaviourand to have to distinguish more case than in the preceding chapter. We formulate ourresults in several lemmas below and summarize them in Table 4.1 at the end.Lemma 4.5. Consider h1 > h2 > 1 > 1/L′ > h3 > 0. Suppose h2h3 < 1 and h1h3 < L′.Then a positive solution cannot contain two or more consecutive decreases.Proof. Assume that there exist two consecutive decreases. Then there exists y0 ≥ L′ suchthat h3y0 ≥ L′. We have h3y0 < max{h1h3, h2h3L′} since y0 < max{h1, h2L′}. Thus, weget h3y0 < L′ since h1h3 < L′ and h2h3 < 1. This contradicts h3y0 ≥ L′. Therefore, apositive solution cannot contain two or more consecutive decreases.Lemma 4.6. Consider h1 > h2 > 1 > 1/L′ > h3 > 0. Suppose that h2h3 < 1 and takek ∈ N such that 1 ≤ hk+12 h3. Then a positive solution cannot contain more than k + 1consecutive increases.Proof. Assume that there are more than k + 1 consecutive increases. Then there existsh3L′ ≤ y0 < L′ and positive integers m, l such that yk+1 = hl1hm2 y0 < L′ where m + l =k + 1. We have, yk+1 = hl1hm2 y0 ≥ hl1hm2 h3L′ since y0 ≥ h3L′. We also have hl1hm2 h3L′ >hl+m2 h3L′ = hk+12 h3L′ since h1 > h2 > 1 and k + 1 = m + l. Together, we get yk+1 ≥hk+12 h3L′ ≥ L′ since hk+12 h3 ≥ 1. This contradicts yk+1 < L′. Therefore, a positive solutioncannot contain more than k + 1 consecutive increases.264.3. ILLUSTRATION OF THE DYNAMICS WITH TWO THRESHOLDSLemma 4.7. Consider h1 > h2 > 1 > 1/L′ > h3 > 0. Suppose that h2h3 < 1 andh1h3 > L′, and take k ∈ N such that h1hk+13 ≤ L′. Then a positive solution cannot containmore than k + 1 consecutive decreases.Proof. Assume that there are more than k+1 consecutive decreases. Then there exists L′ ≤y0 < max{h1, h2L′} such that hk+13 y0 ≥ L′. We have, hk+13 y0 < max{h1hk+13 , h2h3 · hk3L′}since y0 < max{h1, h2L′}. The first term of the maximum is bounded by L′ by the definitionof k. The second term is also bounded by L′ because h2h3 < 1 and h3 < 1. Hence, wefind hk+13 y0 < L′, which contradicts the assumption. Therefore, a positive solution cannotcontain more than k + 1 consecutive decreases.Lemma 4.8. Consider h1 > h2 > 1 > 1/L′ > h3 > 0. Suppose that h2h3 > 1. Then therecannot be two or more consecutive increases. Furthermore, there exist integers k, l suchthat h2hk+13 < 1 and h1hl+13 < L′. Let m = max{k, l}. There cannot be more than m + 1consecutive decreases.Proof. Assume that there are two consecutive increases. Then there exists h3L′ ≤ y0 < L′such that h1y0 < L′ or h2y0 < L′. Since y0 ≥ h3L′ and h1 > h2, we have h1y0 > h2y0 ≥h2h3L′ > L′ by assumption. Hence, we have a contradiction.Next, assume that there are more than m+ 1 consecutive decreases. Then there existsL′ ≤ y0 < max{h1, h2L′} such that ym+1 = hm+13 y0 ≥ L′. Since y0 < max{h1, h2L′},we have ym+1 = hm+13 y0 < max{h1hm+13 , h2hm+13 L′}. The first term of the maximum isbounded by L′ since m ≥ l. The second term is bounded by L′ since m ≥ k. Altogether,we find that ym+1 = hm+13 y0 < L′, which is a contradiction. Therefore, a positive solutioncannot contain more than m+ 1 consecutive decreases.4.3 Illustration of the dynamics with two thresholds4.3.1 Shape of the cyclesGiven the above lemmas, we now know what complex behaviours to expect from ourtwo-threshold model under certain parameter conditions. To illustrate these behaviours, weprovide several simulation results for certain interesting dynamics of this system.We begin with the maximum number of consecutive increases (see Lemma 4.6). We fixparameters h1, h3, and L′ and vary h2. The lemma states that there cannot be more thank + 1 consecutive increases when 1 ≤ hk+12 h3 or, equivalently,k + 1 ≥ − lnh3lnh2. (4.4)This expression is positive since we have h3 < 1 < h2. We note that the maximum numberof consecutive increase increases as h2 decreases to 1. We give three examples in Figure 4.3.For the chosen parameter values (see figure caption), the formula in (4.4) gives k+1 = 2.32274.3. ILLUSTRATION OF THE DYNAMICS WITH TWO THRESHOLDSParameter Conditions Solution Behaviour Lemmah2h3 < 1h1h3 < L′ No more than oneconsecutive decrease4.5h1h3 > L′,h1hk+13 ≤ L′No more than k + 1consecutive decreases4.71 ≤ hk+12 h3 No more than k + 1consecutive increases4.6h2h3 > 1No more than oneconsecutive increase 4.8h2hk+13 < 1,h1hl+13 < L′No more than m+ 1consecutive decreasesTable 4.1: A summary of parameter conditions and their corresponding solution behavioursfor positive solution sequences in the model with two thresholds. We have k, l ∈ N, m =max{k, l}. Furthermore, parameters satisfy the relations h1 > h2 > 1 > 1L′ > h3 > 0.in (a), k + 1 = 4.78 in (b) and k + 1 = 16.88 in (c). Since k has to be an integer, we canhave no more than 3, 5, and 17 consecutive increases, respectively.Figure 4.3 shows the maximum number of consecutive increases that we found numeri-cally by varying the initial condition. Only in Figure 4.3(a) do we find as many consecutiveincreases as the theoretical threshold predicts. In the two other cases, we find fewer. Thisdiscrepancy is to be expected since in the proof of Lemma 4.6, we use the upper boundh2 < h1, which is somewhat close in case (a) but not in (b) and (c).Similarly, we illustrate the maximum number of consecutive decreases (see Lemma 4.7).We choose h3 as our free parameter and fix all others. From the condition in the lemma,we obtaink + 1 =ln L′h1lnh3. (4.5)Given the chosen parameter values in Figure 4.4, we find the maximum number of decreasesto be k + 1 = 2.19 in (a), k + 1 = 1.67 in (b), and k + 1 = 1.25 in (c). Since k has to be aninteger, we cannot have more than 3 consecutive decreases in (a) and no more than 2 in (b)and (c). From the plots in Figure 4.4 that in this case, the theoretically possible maximumnumber of consecutive decreases is actually achieved numerically as well.4.3.2 Coexistence of m and n cyclesFinally, we return to the question of cycles that we studied in the one-threshold model.There, with only two parameters (denoted r and R), we found that we obtained cycles oflength k = l + m if rlRm = 1. Here, we have three parameters, namely h1, h2 and h3.284.3. ILLUSTRATION OF THE DYNAMICS WITH TWO THRESHOLDSWe can choose these parameters such that conditions for cycles of two different lengths aresatisfied simultaneously. For example, we can choose h1h3 = 1 so that a 2-cycle can existand h22h3 = 1 so that a 3-cycle can exist. The 2-cycle will take values alternatingly below 1and above L′. The three cycle will have two values in [1, L′] and one above L′. Dependingon initial conditions, the system approaches either a 2-cycle or a 3-cycle (see Figure 4.5).As in the previous chapter, these cycles are not unique but there exists a continuum of eachof them.Bounds on the initial conditionsWe can analytically determine which initial conditions give rise to a 2-cycle, and whichto a 3-cycle. First, we explain the steps intuitively, and then mathematically.In order to develop our intuition around this calculation, we begin with the 2-cycle.We know that the 2-cycle alternates between points in the first and third regions. Hence,any first-region point in the two cycle must be large enough so that h1y0 falls above thelower threshold of the third region, L′. The lowest possible initial point is the one for whichh1y0 = L′, and the highest possible initial point is the one for which y0 = 1. Hence, theinitial points L′/h1 < y0 ≤ 1 all give rise to a 2-cycle, and the corresponding points in thethird region are L′ < y1 ≤ h1. Once these two regions have been identified, one can findthe initial values outside of the 2-cycle which eventually yield an iterate within a 2-cycle byworking backwards. For example, if we consider the 2-cycle points in region 1, and take onebackwards step, we find that the points L′/h21 < y0 < 1/h1 yield first iterate points insidethe lower 2-cycle region. The other subregions are obtained in an analogous manner. Thisprocess is illustrated in Figure 4.6.We now proceed to the mathematical calculation. We find the bounds of these regions asfollows. In the two threshold model, given the assumption h1 > L′ > h2 > 1 > L′ > h3 > 0,first consider a 2-cycle with the condition h1h3 = 1.Let yl ∈ N be a positive solution sequence of this system that starts below 1 andy0h1 ≥ L′. Then, L′h1 ≤ y0 < 1. We write L′ ≤ y1 < h1 since y1 = y0h1 and h1 > 1. Sinceh1h3 = 1, we get y0 = y2 = y4 = . . . and y1 = y3 = y5 = . . .. Furthermore, a solutionsequence that begins in L′h21≤ y0 < 1h1 hits a 2-cycle after an increase in (0, 1) since y0 < 1and h1 > 1. Similarly, we write, for any n ∈ Z+0 , a solution sequence that begins from[L′hn+11, 1hn1)hits a 2-cycle after n consecutive increases in (0, 1) (see Figure 4.6).Next, let yl ∈ N be a positive solution sequence of this system that starts above L′.Then a solution sequence with initial condition L′h3≤ y0 < h1h3 hits a 2-cycle after a decreasein [L′,∞) since y0 > L′ and 1 > h3. Similarly, we write, for any n ∈ Z+0 , a solution sequencethat starts from[L′hn3, h1hn3)hits a 2-cycle after n consecutive decreases in [L′,∞).Let yl ∈ N be a solution sequence of this system that starts in [1, L′) and hk2y0 ≥ L′ >hk−12 y0 for k ∈ N. If a sequence starts in [1, L′), then the sequence increases by k numberof consecutive increases until yk ≥ L′ > yk−1 since h2 > 1. Subsequently, if yk ∈ [ L′hn+13 ,h1hn3),then the positive solution sequence hits a 2-cycle after n consecutive decreases.294.3. ILLUSTRATION OF THE DYNAMICS WITH TWO THRESHOLDSHence, for any n ∈ N, a solution sequence that starts from [ L′hk2hn+13, h1hk2hn3) hits a 2-cycleafter k consecutive increases in [1, L′) and n consecutive decreases in [L′,∞).Similarly, we can find the bounds of initial value intervals of 3-cycle. Here, the only3-cycles that can occur in the presence of the above 2-cycle have the parameter conditionh22h3 = 1. For k-cycles, where k ≥ 3, the calculation is analogous.Sensitivity to Initial ConditionsWe note that the cycles obtained here are structurally unstable, in that a small pertur-bation in the population numbers will result in the system shifting to a new cycle. Thisbehaviour is also observed in the classic Lotka-Volterra predator-prey model (Kot, 2001).The Rosenzweig-MacArthur and Leslie-Gower-May models, however, exhibit stable limitcycles when they exist (Strohm and Tyson, 2009). The same observation was also made forthe one-threshold model (see Figure 3.6, and accompanying text).304.3. ILLUSTRATION OF THE DYNAMICS WITH TWO THRESHOLDS(a) (b)(c)Figure 4.3: Sample cycles illustrating the maximum number of consecutive increases(Lemma 4.6). Note that for all of the pseudo-cycles shown, there are sometimes multi-ple consecutive increase steps, but never more than one decrease step. According to (4.4),there cannot be more than 3 (a), 5 (b), and 17 (c) consecutive increases. This requirementis fulfilled by all of the pseudo-cycles shown. Parameter values are h2 = 2 (a), h2 = 1.4 (b),and h2 = 1.1 (c). Fixed parameters are h1 = 2.3, h3 = 0.2 and L′ = 3. Black dashed linesindicate the thresholds. The pink dashed and solid lines give the upper and lower bound ofthe invariant interval.314.3. ILLUSTRATION OF THE DYNAMICS WITH TWO THRESHOLDS(a) (b)(c)Figure 4.4: Sample cycles illustrating the maximum number of consecutive decreases(Lemma 4.7). Note that for all of the pseudo-cycles shown, there are sometimes multi-ple consecutive decrease steps, but never more than one increase step. According to (4.5),for the cases shown there cannot be more than 3 (a), 2 (b), and 2 (c) consecutive de-creases. This requirement is fulfilled by all of the pseudo-cycles shown. Parameter valuesare h3 = 0.4 (a), h3 = 0.3 (b), and h3 = 0.2 (c). Fixed parameters are h1 = 15, h2 = 2 andL′ = 2. Black dashed lines indicate the thresholds. The pink dashed and solid lines givethe upper and lower bound of the invariant interval.324.3. ILLUSTRATION OF THE DYNAMICS WITH TWO THRESHOLDS(a) (b)(c)Figure 4.5: Illustration of coexisting cycles. (a) Solutions form a 2-cycle from an initialvalue at 0.4. (b) Solutions form a 3-cycle from an initial value at 0.7. (c) Asymptoticoutcome depending on initial condition. Parameters are h1 = 4, h2 = 2, h3 = 0.25 andL′ = 3. Black dashed lines indicate the thresholds. Pink dashed and solid lines correspondto the upper and lower bound of the invariant interval.334.3. ILLUSTRATION OF THE DYNAMICS WITH TWO THRESHOLDSFigure 4.6: Illustration of the calculations to determine the initial values that give rise to2-cycles for the following parameter values: h1 = 4, h2 = 2, h3 = 0.25 and L′ = 3. Bluedots: Initial values that give rise to 2-cycles. Red dots: Initial values that give rise to3-cycles.34Chapter 5Discussion and conclusionsThe mathematical models in our analysis are inspired by the lemming-arctic fox-snowyowl predator-prey system on Bylot Island. We assume a constant adult fox population dueto the approximately fixed number of available territories. Similarly, we assume a constantsnowy owl population, if present, because it is also a territorial predator. Their presenceon the island depends on sufficient available prey. In the study, the model parameterssuch as birth, death and predation remain fixed during the winter and summer seasons.However, the predation pressure changes when the prey population crosses the thresholdvalues. The interaction terms between the species of this model follow the non-linear massaction dynamics in the Lotka-Volterra model. Our assumptions allow us to reduce thesystem of nonlinear ordinary differential equations to first-order difference equations. Thisdiscrete time representation of the system allows us to analyse the lemming cycles in detailand find new insights regarding the cycles in this predator-prey system.Our first result is for the solution behaviour of the model with one threshold. To arriveat this model, we assume that the predation pressure of the resident predator (the arcticfox) is approximately similar in the summer and winter seasons. It is clear that steady statesdo not exist except if R = r = 1. For all other cases, the solutions are cyclic; either theyare periodic or quasi-periodic (i.e., pseudo-cycles). We explore the mechanisms behind thenumber of increases and decreases in the periodic oscillations. Next, we provide a theory forthe asymptotic behaviour of the solutions. We find that there is variation in the structureof the increases and decreases within each period; the specific structure observed dependson precise conditions involving the prey growth rates and threshold values. Furthermore,we demonstrate that these cyclic populations eventually enter an invariant interval, and wedetermine the upper and lower bounds of this interval.Subsequently, we allow the arctic fox predation pressure to change seasonally. We findnew dynamical behaviour in the solutions of this second model, e.g., coexistence of twostable cycles. As the initial condition changes, solutions of the system with two thresholdscan approach, for example, a 2-cycle or a 3-cycle.We demonstrate that our simplified mathematical model is capable of exhibiting complexbehaviour while maintaining some defined system structure. Specifically, we present lemmasthat determine the maximum number of increases/decreases in a population pattern in termsof the parameter values. For instance, we present that, given growth rates and thresholdvalues, there cannot be more than one population increase in a solution sequence. Inaddition, we find the maximum number of consecutive decreases in the solution sequence.This pattern is distinct from some classical model population patterns. The Rosenzweig-MacArthur model exhibits cycles where the prey population has gradual increases followed35by gradual decreases (Rosenzweig and MacArthur, 1963). On the other hand, the Rickermodel exhibits cycles of alternating population increase and decrease (Kot, 2001).This study advances scientific knowledge by providing a mathematical analysis of athree-species predator-prey system in which the structure of the food web is determinedby prey density. This work is inspired by the work presented in Hutchison et al. (2020),which is a study of the interactions between resident (lemming and arctic fox) and migrant(snowy owl and snow goose) species on Bylot Island. According to those authors, specifi-cally, there is a need for theory regarding predator-prey systems in seasonal food webs. Infact, arctic food webs experience drastic seasonal changes that have strong effects on preyabundance, predation pressure, and migration patterns (Therrien et al., 2014). In addition,seasonal changes in, for example, snowy owl migration patterns affect other food webs inwhich they are involved. Therefore, analyzing the seasonal changes among predator andprey communities on Bylot Island helps us identify the mechanisms influencing populationdynamics and their adaptation to those changes.We also ask what these results mean in terms of Bylot Island. The communities onBylot Island, are subjected to regular population fluctuations rather than reaching a sta-ble equilibrium and the food web on Bylot Island experiences structural changes while thespecies adaptto these changes. The bounds on the system solutions indicate that the lem-ming population does not go extinct or continue to grow over time, even though lemminggrowth is exponential at low density. More importantly, the coexistence of cycles such as2-cycle and 3-cycle in the model with two thresholds implies that the species in this foodweb may be sensitive to the initial population densities on Bylot Island. This explains theLFO system’s asymptotic behaviour. Moreover, this may explain why a seasonal systemlike arctic system tends to show predictable cyclic patterns in the long run depending onthe initial population.Our work leads to several open questions which would expand our study in this area. Tobegin with, we envision using this approach to understand the 4 species LFO-snow goosesystem studied in Hutchison et al. (2020). With this shift in the ecological scenario, anadditional threshold might appear in a corresponding model. This will induce changes inthe periodic oscillations of the solutions. Another future interest of this study is allowingvariable fox and/or owl population densities. Studying varying predator populations wouldgive us some insights into the factors that regulate predator population sizes. Anotherfuture work includes stochasticity in mathematical models as such models would be morerealistic. For instance, protection of lemmings from predators due to the existence of snowlayers (Bilodeau et al., 2013) could be presented as a noise term of the lemming growthrate in the winter. When snow layers are present, lemmings construct tunnels under them,making predation harder. This protection allows lemming populations to grow faster thanin the summer, as predators like arctic foxes struggle to predate them. Overall, our studysets the stage for further exploration of predator-prey systems to provide more theory tobetter understand seasonal predator-prey systems and their stability.36BibliographyAbbott, K. C. (2007). Does the pattern of population synchrony through space reveal ifthe moran effect is acting? Oikos, 116(6):903–912. 1Adams, L. (1959). An analysis of a population of snowshoe hares in northwestern Montana.Ecological Monographs, 29(2):141–170. 1Alebraheem, J. and Hasan, Y. A. (2013). Dynamics of a two predator–one prey system.Computational and Applied Mathematics, 33(3):767–780. 2, 3Allen, L. J. S. (2007). An introduction to mathematical biology. Pearson/Prentice Hall. 10Angerbjo¨rn, A., Tannerfeldt, M., and Lundberg, H. (2001). Geographical and temporalpatterns of lemming population dynamics in Fennoscandia. Ecography, 24(3):298–308.1, 4Barraquand, F., Louca, S., Abbott, K. C., Cobbold, C. A., Cordoleani, F., DeAngelis, D. L.,Elderd, B. D., Fox, J. W., Greenwood, P., Hilker, F. M., Murray, D. L., Stieha, C. R.,Taylor, R. A., Vitense, K., Wolkowicz, G. S., and Tyson, R. C. (2017). Moving forwardin circles: challenges and opportunities in modelling population cycles. Ecology Letters,20(8):1074–1092. 1Bauer, S. and Hoye, B. J. (2014). Migratory animals couple biodiversity and ecosystemfunctioning worldwide. Science, 344(6179). 3Berlow, E. L., Neutel, A., Cohen, J. E., De Ruiter, P. C., Ebenman, B., Emmerson, M.,Fox, J. W., Jansen, V. A. A., Iwan Jones, J., Kokkoris, G. D., Logofet, D. O., McKane,A. J., Montoya, J. M., and Petchey, O. (2004). Interaction strengths in food webs: issuesand opportunities. Journal of Animal Ecology, 73(3):585–598. 3Bilodeau, F., Kenney, A. J., Gilbert, B. S., Hofer, E., Gauthier, G., Reid, D. G., Berteaux,D., and Krebs, C. J. (2013). Evaluation of a technique to trap lemmings under the snow.ARCTIC, 66(1). 36Bolduc, D., Fauteux, D., Gagnon, C. A., Gauthier, G., Beˆty, J., and Legagneux, P. (2023).Testimonials to reconstruct past abundances of wildlife populations. Basic and AppliedEcology, 68:23–34. 3Boxall, P. C. and Lein, M. R. (1982). Territoriality and habitat selection of female snowyowls (nyctea scandiaca) in winter. Canadian Journal of Zoology, 60(10):2344–2350. 437BibliographyChen, Y. and Zhang, F. (2013). Dynamics of a delayed predator–prey model with predatormigration. Applied Mathematical Modelling, 37(3):1400–1412. 3Debnath, S., Majumdar, P., Sarkar, S., and Ghosh, U. (2023). Memory effect onprey–predator dynamics: Exploring the role of fear effect, additional food and anti-predator behaviour of prey. Journal of Computational Science, 66:101929. 3Dostalkova, I., Kindlmann, P., and Dixon, A. F. (2002). Are classical predator–prey modelsrelevant to the real world? Journal of Theoretical Biology, 218(3):323–330. 2Duchesne, D., Gauthier, G., and Berteaux, D. (2011). Evaluation of a method to deter-mine the breeding activity of lemmings in their winter nests. Journal of Mammalogy,92(3):511–516. 4Elton, C. S. (1924). Periodic fluctuations in the numbers of animals: Their causes andeffects. The Journal of Experimental Biology, 2:119–163. 1, 3, 4Farhan, A. G., Khalaf, N. S., and Aldhlki, T. J. (2023). Dynamical behaviors of fourdimensional prey-predator model. European Journal of Pure and Applied Mathematics,16(2):899–918. 3Fauteux, D., Gauthier, G., and Berteaux, D. (2015). Seasonal demography of a cycliclemming population in the Canadian Arctic. Journal of Animal Ecology, 84(5):1412–1422. 1, 4Freedman, H. and Waltman, P. (1984). Persistence in models of three interacting predator-prey populations. Mathematical Biosciences, 68(2):213–231. 2, 3Gakkhar, S., Negi, K., and Sahani, S. K. (2009). Effects of seasonal growth on ratiodependent delayed prey predator system. Communications in Nonlinear Science andNumerical Simulation, 14(3):850–862. 3Gao, X., Pan, Q., He, M., and Kang, Y. (2013). A predator–prey model with diseasesin both prey and predator. Physica A: Statistical Mechanics and its Applications,392(23):5898–5906. 3Gauthier, G. (2004). Trophic interactions in a high arctic snow goose colony. Integrativeand Comparative Biology, 44(2):119–129. 2, 3, 4Gauthier, G., Berteaux, D., Beˆty, J., Tarroux, A., Therrien, J.-F., McKinnon, L., Legag-neux, P., and Cadieux, M.-C. (2011). The tundra food web of Bylot Island in a changingclimate and the role of exchanges between ecosystems. E´coscience, 18(3):223–235. 1, 4Gilg, O., Hanski, I., and Sittler, B. (2003). Cyclic dynamics in a simple vertebrate predator-prey community. Science, 302(5646):866–868. 3, 4Gilg, O., Sittler, B., and Hanski, l. (2009). Climate change and cyclic predator–prey popu-lation dynamics in the high arctic. Global Change Biology, 15(11):2634–2652. 338BibliographyGoltsman, M., Kruchenkova, E. P., Sergeev, S., Johnson, P. J., and Macdonald, D. W.(2005). Effects of food availability on dispersal and cub sex ratio in the Mednyi Arcticfox. Behavioral Ecology and Sociobiology, 59(2):198–206. 8Gravel, R., Lai, S., and Berteaux, D. (2023). Long-term satellite tracking reveals patterns oflong-distance dispersal in juvenile and adult arctic foxes (vulpes lagopus). Royal SocietyOpen Science, 10(2). 4Gruyer, N., Gauthier, G., and Berteaux, D. (2009). Demography of two lemming specieson Bylot Island, Nunavut, Canada. Polar Biology, 33(6):725–736. 1, 4Hang, L., Zhang, L., Wang, X., Li, H., and Teng, Z. (2020). A hybrid predator–preymodel with general functional responses under seasonal succession alternating betweenGompertz and logistic growth. Advances in Difference Equations, 2020(1). 3, 4Hanski, I., Hansson, L., and Henttonen, H. (1991). Specialist predators, generalist predators,and the microtine rodent cycle. The Journal of Animal Ecology, 60(1):353. 1Hanski, I., Henttonen, H., Korpima¨ki, E., Oksanen, L., and Turchin, P. (2001). Small-rodentdynamics and predation. Ecology, 82(6):1505–1520. 1, 4Hanski, I. and Korpima¨ki, E. (1995). Microtine rodent dynamics in northern Europe:Parameterized models for the predator-prey interaction. Ecology, 76(3):840–850. 1Hansson, L. and Henttonen, H. (1988). Rodent dynamics as community processes. Trendsin Ecology and Evolution, 3(8):195–200. 1Harvell, C. D., Mitchell, C. E., Ward, J. R., Altizer, S., Dobson, A. P., Ostfeld, R. S., andSamuel, M. D. (2002). Climate warming and disease risks for terrestrial and marine biota.Science, 296(5576):2158–2162. 3Hastings, A. and Powell, T. (1991). Chaos in a three-species food chain. Ecology,72(3):896–903. 2, 4Holt, R. D. and Polis, G. A. (1997). A theoretical framework for intraguild predation. TheAmerican Naturalist, 149(4):745–764. 2, 3Ho¨rnfeldt, B. (1978). Synchronous population fluctuations in voles, small game, owls, andtularemia in northern Sweden. Oecologia, 32(2):141–152. 1Hutchison, C., Guichard, F., Legagneux, P., Gauthier, G., Beˆty, J., Berteaux, D., Fau-teux, D., and Gravel, D. (2020). Seasonal food webs with migrations: multi-seasonmodels reveal indirect species interactions in the Canadian Arctic tundra. PhilosophicalTransactions of the Royal Society A: Mathematical, Physical and Engineering Sciences,378(2181):20190354. 3, 4, 5, 36Ims, R., Henden, J., and Killengreen, S. (2008). Collapsing population cycles. Trends inEcology and Evolution, 23(2):79–86. 1, 239BibliographyIms, R. A. and Fuglei, E. (2005). Trophic interaction cycles in tundra ecosystems and theimpact of climate change. BioScience, 55(4):311. 1, 3, 4Ji, J. and Wang, L. (2022). Competitive exclusion and coexistence in an intraguild predationmodel with Beddington–DeAngelis functional response. Communications in NonlinearScience and Numerical Simulation, 107:106192. 3Jin, M., Xu, F., Shen, C.-S., Zhang, J.-Q., and Wang, C.-Y. (2019). Nontrivial effect oftime-varying migration on the three species prey-predator system. Communications inTheoretical Physics, 71(1):127. 3Keith, C. and Rongstad, B. M. C. (1984). Demography and ecology of a declining snowshoehare population. Wildlife Monographs, 90:3–43. 1Khaliq, A., Ibrahim, T. F., Alotaibi, A. M., Shoaib, M., and El-Moneam, M. A.(2022). Dynamical analysis of discrete-time two-predators one-prey Lotka–Volterramodel. Mathematics, 10(21):4015. 3Kharitonov, S. P., Volkov, A. E., Willems, F., van Kleef, H., Klaassen, R. H. G., Nowak,D. J., Nowak, A. I., and Bublichenko, A. G. (2008). Brent goose colonies near snowyowls: Internest distances in relation to abundance of lemmings and arctic foxes. BiologyBulletin, 35(3):270–278. 3Klemola, T., Koivula, M., Korpima¨ki, E., and Norrdahl, K. (1997). Small mustelid predationslows population growth of microtus voles: A predator reduction experiment. The Journalof Animal Ecology, 66(5):607. 1Kot, M. (2001). Elements of Mathematical Ecology. Cambridge University Press. 2, 10,30, 36Krebs, C. J., Boonstra, R., and Boutin, S. (2017). Using experimentation to understandthe 10-year snowshoe hare cycle in the boreal forest of North America. Journal of AnimalEcology, 87(1):87–100. 1, 3Krebs, C. J., Boonstra, R., Boutin, S., and Sinclair, A. (2001). What drives the 10-yearcycle of snowshoe hares? BioScience, 51(1):25. 1Kumar, A., Jani, K., Khatri, D., Patel, H., Sharma, A. K., and Sahu, N. K. (2021).A computer-based simulation showing balance of the population of predator and preyand the effects of human intervention. IOP Conference Series: Materials Science andEngineering, 1031(1):012049. 3Lai, S., Beˆty, J., and Berteaux, D. (2016). Movement tactics of a mobile predator ina meta-ecosystem with fluctuating resources: the arctic fox in the high arctic. Oikos,126(7):937–947. 840BibliographyLatimer, C. E. and Zuckerberg, B. (2019). How extreme is extreme? Demographic ap-proaches inform the occurrence and ecological relevance of extreme events. EcologicalMonographs, 89(4). 4Legagneux, P., Gauthier, G., Berteaux, D., Beˆty, J., Cadieux, M.-C., Bilodeau, F., Bolduc,E., McKinnon, L., Tarroux, A., Therrien, J.-F., Morissette, L., and Krebs, C. J. (2012).Disentangling trophic relationships in a high arctic tundra ecosystem through food webmodeling. Ecology, 93(7):1707–1716. 4Lin, J. and Kahn, P. B. (1978). Qualitative dynamics of three species predator-prey systems.Journal of Mathematical Biology, 5(3):257–268. 2Lindstro¨m, J., Ranta, E., Kokko, H., Lundberg, P., and Kaitala, V. (2001). From arcticlemmings to adaptive dynamics: Charles Elton’s legacy in population ecology. BiologicalReviews, 76(1):129–158. 1, 3Montoya, J. M., Rodr´ıguez, M. A., and Hawkins, B. A. (2003). Food web complexity andhigher-level ecosystem services. Ecology Letters, 6(7):587–593. 3Murdoch, W. W., Kendall, B. E., Nisbet, R. M., Briggs, C. J., McCauley, E., and Bolser,R. (2002). Single-species models for many-species food webs. Nature, 417(6887):541–543.1Murray (2002). Mathematical Biology: I. An Introduction. Springer New York. 2Murray, J. D. (1993). Mathematical Biology. Springer Berlin Heidelberg. 2Myers, J. H. (2018). Population cycles: generalities, exceptions and remaining mysteries.Proceedings of the Royal Society B: Biological Sciences, 285(1875):20172841. 1Norrdahl, K. and Korpima¨ki, E. (1995). Effects of predator removal on vertebrate preypopulations: birds of prey and small mammals. Oecologia, 103(2):241–248. 4O’Donoghue, M., Boutin, S., Krebs, C. J., and Hofer, E. J. (1997). Numerical responses ofcoyotes and lynx to the snowshoe hare cycle. Oikos, 80(1):150. 1Pattanayak, D., Mishra, A., Dana, S. K., and Bairagi, N. (2021). Bistability in a tri-trophicfood chain model: Basin stability perspective. Chaos: An Interdisciplinary Journal ofNonlinear Science, 31(7). 2Post, E., Forchhammer, M. C., Bret-Harte, M. S., Callaghan, T. V., Christensen, T. R.,Elberling, B., Fox, A. D., Gilg, O., Hik, D. S., Høye, T. T., Ims, R. A., Jeppesen, E.,Klein, D. R., Madsen, J., McGuire, A. D., Rysgaard, S., Schindler, D. E., Stirling, I.,Tamstorf, M. P., Tyler, N. J., van der Wal, R., Welker, J., Wookey, P. A., Schmidt,N. M., and Aastrup, P. (2009). Ecological dynamics across the arctic associated withrecent climate change. Science, 325(5946):1355–1358. 441BibliographyReid, D. G., Krebs, C. J., and Kenney, A. (1995). Limitation of collared lemming populationgrowth at low densities by predation mortality. Oikos, 73(3):387. 4Rosenzweig, M. L. (1973). Exploitation in three trophic levels. The American Naturalist,107(954):275–294. 2Rosenzweig, M. L. and MacArthur, R. H. (1963). Graphical representation and stabilityconditions of predator-prey interactions. The American Naturalist, 97(895):209–223. 2,36Samelius, G. and Alisauskas, R. T. (2017). Components of population growth for arctic foxesat a large arctic goose colony: the relative contributions of adult survival and recruitment.Polar Research, 36(sup1):6. 4Smith, M. J., White, A., Sherratt, J. A., Telfer, S., Begon, M., and Lambin, X. (2007).Disease effects on reproduction can cause population cycles in seasonal environments.Journal of Animal Ecology, 77(2):378–389. 1Stenseth, N. C. (1999). Population cycles in voles and lemmings: Density dependence andphase dependence in a stochastic world. Oikos, 87(3):427. 1, 3, 4Strohm, S. and Tyson, R. (2009). The effect of habitat fragmentation on cyclic populationdynamics: A numerical study. Bulletin of Mathematical Biology, pages 1323–1348. 2, 30Sun, G.-Q., Jin, Z., Liu, Q.-X., and Li, L. (2008). Dynamical complexity of a spatialpredator–prey model with migration. Ecological Modelling, 219(1–2):248–255. 3Tannerfeldt, M., Angerbjo¨rn, A., and Angerbjorn, A. (1998). Fluctuating resources andthe evolution of litter size in the arctic fox. Oikos, 83(3):545. 4Teramoto, E., Kawasaki, K., and Shigesada, N. (1979). Switching effect of predation oncompetitive prey species. Journal of Theoretical Biology, 79(3):303–315. 2Therrien, J.-F., Gauthier, G., Korpima¨ki, E., and Beˆty, J. (2014). Predation pressure byavian predators suggests summer limitation of small-mammal populations in the Cana-dian Arctic. Ecology, 95(1):56–67. 4, 36Turchin, P. (2013). Complex Population Dynamics: A Theoretical/Empirical Synthesis(MPB-35). Princeton University Press. 2, 4Tyson, R., Haines, S., and Hodges, K. (2010). Modelling the Canada lynx and snowshoehare population cycle: the role of specialist predators. Theoretical Ecology, 3:97–111. 3Tyson, R. and Lutscher, F. (2016). Seasonally varying predation behavior and climate shiftsare predicted to affect predator-prey cycles. The American Naturalist, 188(5):539–553.3, 442BibliographyVance-Chalcraft, H. D., Rosenheim, J. A., Vonesh, J. R., Osenberg, C. W., and Sih, A.(2007). The influence of intraguild predation on prey suppression and prey release: Ameta-analysis. Ecology, 88(11):2689–2696. 3Wang, Q. (2019). Qualitative analysis of a Lotka–Volterra predator–prey system withmigration. Journal of Mathematical Analysis and Applications, 472(1):421–431. 3Wang, X. and Zou, X. (2017). Modeling the fear effect in predator–prey interactions withadaptive avoidance of predators. Bulletin of Mathematical Biology, 79(6):1325–1359. 343
Thesis/Dissertation
eng
Mathematics
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Tyson, Rebecca
Lutscher, Frithjof
Predator-prey interactions in Arctic Canada : a discrete-time model of lemmings, foxes and snowy owls on Bylot Island
Text