## Abstract

Seasonality is an important component in many population systems, and factors such as latitude, altitude and proximity to the coastline affect the extent of the seasonal fluctuations. In this paper, we ask how changes in seasonal fluctuations impact on the population cycles. We use the Fennoscandian vole system as a case study, focusing on variations in the length of the breeding season. We use a predator–prey model that includes generalist and specialist predation alongside seasonal forcing. Using a combination of bifurcation analysis and direct simulations, we consider the effects of varying both the level of generalist predation and the length of the breeding season; these are the main changes that occur over a latitudinal gradient in Fennoscandia. We predict that varying the breeding season length leads to changes in the period of the multi-year cycles, with a higher period for shorter breeding season lengths. This concurs with the gradient of periodicity found in Fennoscandia. The Fennoscandian vole system is only one of many populations that are affected by geographical and temporal changes in seasonality; thus our results highlight the importance of considering these changes in other population systems.

## 1. Introduction

Seasonal forcing represents a pervasive source of environmental variability in natural systems, with many species exhibiting seasonal changes in their life-history parameters [1]. This can affect both epidemiological and interacting population dynamics and has been proposed as a cause of multi-year cycles in these systems [2,3]. Considerable focus has been placed on the impact of seasonality in predator–prey systems, and theory has been developed both for general model frameworks [4,5] and for specific wildlife populations [6–8]. This has shown that seasonal forcing, alongside the nonlinear dynamics, can lead to complicated population dynamics including multi-year cycles, quasi-periodic solutions and chaos.

Many populations span large geographical areas over which the characteristic features of the seasonal fluctuations can vary considerably [9–12]. This variation is, for instance, dependent on latitude, altitude, proximity to the coastline and prevailing weather patterns. There is also evidence that climate change has led to changes in seasonal patterns such as the earlier onset of breeding and increased variability in climatic conditions [13]. It is, therefore, important to extend the analysis of the effects of seasonality to assess the impact of changes in the amplitude and functional form of seasonal forcing on the underlying population dynamics. In particular, to what extent do modifications in seasonal patterns lead to shifts in population dynamics?

While numerous examples of widespread populations exist in which seasonal forcing and predator–prey dynamics combine to produce cyclic population dynamics, in order to examine the impact of changes in seasonality, we will focus on a well-known case study—namely the Fennoscandian vole system [14–18]. This system spans a large geographical area and undergoes high-amplitude seasonal forcing that varies with location (typically, the length of the breeding season decreases with latitude). The important features of this vole system are the large amplitude, multi-annual cycles north of 60° N, synchronous over large spatial scales with a period of 3–5 years, with the 5 year cycles occurring furthest north. In comparison, the southern Fennoscandian region has low-amplitude ‘non-cyclic’ seasonal fluctuations, effectively an annual cycle. Therefore, period and amplitude decrease as one moves along a north–south gradient in Fennoscandia [9,14,15,19].

It has been suggested that changes in generalist predation with latitude explain this gradient within Fennoscandia [17,20,21] and this is supported by theory [6,7]. However, the seasonal forcing term remains fixed in these models. This contrasts with the pronounced differences in the pattern of seasonality between the southern and northern parts of Fennoscandia, and the effect of this on the predation hypothesis has received almost no attention (a notable exception being Dalkvist *et al*. [22] who consider three different breeding season lengths). In southern Fennoscandia, the breeding season can be over seven months long, whereas in northern areas it is usually three to four months [16,22,23]. Furthermore, there is evidence to suggest that the vole cycles have been reducing in period length because of climate change [24,25]. The vole system can be used to see whether geographical variations in seasonal forcing can help explain the gradient in cyclicity which is currently found in Fennoscandia, as well as how changes in the pattern of seasonality will alter the population dynamics. To capture the changing seasonality, we will vary breeding season length. Using the vole system as a case study, we will highlight the wide range of complex population dynamics that can arise owing to seasonal forcing, and indicate how variation in seasonality can generate shifts in population dynamics.

## 2. Model

A wide range of intrinsic and extrinsic factors have been proposed to explain the geographical gradient in cyclicity found in Fennoscandian voles, but the most widely supported is the predation hypothesis, first formulated by Hanski *et al*. [20] and further discussed in the earlier studies [7,17,19]. This states that the northern cycles are caused by specialist predation, usually by the least weasel. As the weasels have little alternative prey, they have delayed reactions in their reproduction in response to changes in vole numbers, generating the cycles (although food limitation may play a role in determining the size of the prey peak [26]). Both the protective layer of snow cover, which hinders predators unable to tunnel under the snow, and the low diversity of prey in these northerly locations mean that there are few generalist predators. In the south there are numerous generalist predators, such as birds, badgers and foxes. In fact, 80 per cent of rodent mortality by predation in the south is by generalists [15]. Because generalist predators exhibit rapid prey switching when vole numbers drop, this type of predation has a stabilizing effect on the vole population. The predation hypothesis effectively states that specialist predators induce cycles, while the generalist predators in the south act to dampen these fluctuations. This is backed up by predator exclusion experiments [27,28] as well as through analysis of many of the proposed reasons for the cycles (such as the predation hypothesis, landscape fragmentation, maternal effects) using statistical methods to assess the extent of their agreement with the vole data collected [9,21,29].

### (a) Model details

Our work is based on the model of Turchin & Hanski [7], which incorporates both specialist and generalist predation. Importantly, we introduce seasonal forcing to the growth rates and implicitly also the carrying capacities of both prey and predator through a seasonal term *S*(*t*), which includes breeding season length (see below for details). When suitably non-dimensionalized (see the electronic supplementary material for details) the model is
2.1Here, *x*(*t*) and *y*(*t*) are the densities of prey and specialist predator at time *t*, respectively. The prey undergo logistic growth. They are affected by two predation terms, the first owing to generalist and the second specialist predation. Generalist predation is a Holling type III functional form, because generalists will switch to other prey items when vole numbers are low. The specialist predation is the Holling type II functional form which incorporates handling time of prey. The predators have a logistic growth with growth rate *s* and a carrying capacity which is determined by the density of prey. We use the parameter values of Turchin & Hanski [7] namely *r* = 6, *s* = 1.25, *d* = 0.04, *a* = 15 and *h* = 0.1. These values were chosen based on a combination of field data, previous literature and time-series analysis. The parameter *g* represents the level of generalist predation, which changes along a north–south gradient. We vary *g* between 0 and 1, because Turchin & Hanski [7] estimate the value of *g* to be in the range 0.933–0.066 at four vole localities, ranging from Revinge, Sweden at 56° N to Kilpisjärvi, Finland at 69° N. Note that small *g* indicates less generalist predation, corresponding to a more northerly location.

The main aim of this paper is to consider how geographical changes in seasonal forcing affect the population dynamics and we do this by varying the breeding season length. We introduce a parameter *l* as part of the seasonal term *S*(*t*) which is given by the following equation (see the electronic supplementary material for more details):
2.2The parameter *l* determines the length of the breeding season, defined as when the forced growth term is above its unforced value: for the prey this is when *rS*(*t*) > 6. Figure 1 shows how our forcing term changes depending on the value of *l*. As *l* increases, then the breeding season length decreases. Although we vary *l* from 0 to 5, the most biologically relevant part of the range is 0.5–3.9, which gives a breeding season between three and eight months long; the data discussed above shows that this is the appropriate variation across Fennoscandia.

In addition, figure 1 shows that as *l* increases, the mean value of the growth rate decreases. This reflects the fact that if the breeding season is longer, voles will produce more litters and female young from the first one or two litters have time to mature and also are able to breed [30]. Turchin & Ostfeld [31] suggest 14 as a maximum value for the growth rate and our values are below this limit. However, this does mean our results are based on both breeding season length and mean value of the growth rate varying as we vary *l*, and we do not determine how each component affects the results separately.

### (b) Model analysis

We studied the model using two different methods—a bifurcation and a simulation approach, following the procedure in Taylor *et al*. [32]. The bifurcation results were calculated using essentially standard numerical continuation techniques implemented via the software Auto [33]. This method determines parameter ranges in which different multi-year cycles exist, and when they are stable. The simulation results were produced by solving the model equations (2.1) using Matlab (ode15s) for 2000 years and the solutions were tested to determine whether there existed a periodic solution of 1–9 years; if not, the solution was labelled as quasi-periodic, meaning that the solution was approximately periodic, but does not have a finite period. We also used fast Fourier transforms to determine the period most closely reflected in the solution (‘dominant period’); this was especially helpful when the cycle was quasi-periodic. We performed our simulations over a grid of points in parameter space, with 50 simulations at each point to illustrate the relative sizes of the basins of attraction of the different solutions. See electronic supplementary material for further details.

## 3. Results

Hanski *et al*. [6] and Turchin & Hanski [7] focused on the case where breeding season is exactly six months long (*l* = 1), and considered the effects of varying the extent of generalist predation (*g*). This revealed a change from annual to multi-annual cycles (period 3–4 years), as the level of generalist predation is reduced. We illustrate this in figure 2 which shows simulations for *l* = 1 and various values of *g*, for two sets of initial conditions. Note that for *g* = 0.5 3 year cycles occur for one set of initial conditions and annual cycles for the other. These results highlight the existence of multiple solutions, with strong sensitivity to initial conditions. Clearly bifurcation analysis and/or extensive simulation are necessary to reveal the full range of model predictions.

We used the bifurcation method to investigate potential dynamics as both season length *l* and generalist predation level *g* are varied (figure 3*a*). This divides the *l*–*g* parameter plane into different regions denoting the different possible dynamics. This shows that for moderate or high *l* and *g* values only annual cycles are possible, with multi-year and quasi-periodic solutions for lower *l* and *g*. Figure 3*a* therefore shows that when the breeding season is longer there is a larger range of generalist predation levels giving multi-year cycles. However, this figure does not show the period of the cycles as a function of parameters.

To examine the solutions in more detail, we show parameter regions giving both existence and stability for 2, 3, 4 and 5 year cycles (figure 4). For cycles with higher period (6–9 years) the results are similar to figure 4*d*, although the regions tend to get smaller with increasing period. The cycles exist inside the boundary curves, with stability/instability shown by the thick black/thin grey lines, respectively. The coloured curves represent different types of bifurcation: period-doubling bifurcation (blue), saddle–node (aka fold or tangent; red) and Neimark–Sacker (green); a Neimark–Sacker bifurcation is a discrete version of a Hopf bifurcation. (More detail and a full bifurcation diagram are provided in the electronic supplementary material.)

Each of the regions in figure 4 has a complicated form; in particular for the 3 and 4 year cycles, there are two separate regions of stability. From these diagrams we can see that cycles of longer period occur only for lower values of *g*. There are many areas where the cycles are unstable and here we expect either quasi-periodic or annual solutions to be stable. The cyclic regions also intersect (see the electronic supplementary material, figure S3.2) and so there are parameter regions that give rise to coexisting stable multi-year cycles. Therefore, it is important to consider which cycles actually occur for different values of *l* and *g*.

Simulation results displaying the population dynamical outcome in parameter space (figure 3*b*) provide a clearer view of how the multi-year solutions overlap and indicate the relative sizes of the basins of attraction when multiple solutions coexist. The results compare favourably with the bifurcation picture in figure 3*a*, with a relatively clear split between the yearly and quasi-periodic solutions interspersed with multi-year solutions. Owing to the large regions of instability within the multi-year cycle regions shown in figure 4, there were not many regions in figure 3*b* where two multi-year cycles coexisted and were both stable. Thus there were hardly any pie charts which showed multiple multi-year stable solutions, although in one case both 4 and 9 year cycles occur (where *l* = 1.5, *g* = 0.2). However, for a number of parameter pairs multi-year solutions coexist with annual or quasi-periodic solutions.

Figure 3*a* shows that when the breeding season is longer, there is a larger range of generalist predation levels giving multi-year cycles. But the multi-year cyclic dynamics in Fennoscandia occur in the areas where the breeding season is shorter. To produce further detail on the nature of the population dynamics in regions where quasi-periodicity is predicted, we use fast Fourier transforms to analyse each simulation, revealing the dominant cycle period (figure 3*c*). As *l* is increased for low *g*, i.e. as the breeding season length is reduced, the periodicity of the solutions increases. This concurs with the gradient of periodicity that has been found in Fennoscandia where in mid-Fennoscandia, there are 3 year cycles and further north (with shorter season lengths) there are 5 year cycles [15,25].

## 4. Discussion

The basic conclusion from previous modelling of the predation hypothesis for Fennoscandian vole cycles is that cycle period increases as the generalist predation level decreases. Our work shows that this is an over-simplification. Even with fixed breeding season length, the existence of multiple stable solutions implies that switches between, say, annual and 3 year cycles are possible. When variations in both season length and generalist predation level are considered, there is a rich variety of potential dynamics, with coexisting multi-year solutions and numerous possible cycle periods. This includes parameter values showing both 4 and 9 year cycles, an interesting phenomenon considering that short-term and long-term cycles are often seen as either–or dynamics, rather than both being possible for the same system [19].

Empirical support for the dependence of vole cycle period on breeding season length is provided by the work by Strann *et al*. [34] on voles in Kirkesdalen, Norway. Here, the population undergoes cycles with a period of 3 years, which contrasts with the 5 year cycles at other sites with a similar latitude. Strann *et al*. [34] attributed this difference to climatic differences, owing to Kirkesdalen's proximity to the coast. Our results suggest that a key implication of this climatic difference might be a longer breeding season; our model predicts that this would indeed cause lower period cycles.

The variation in vole breeding season length in Fennoscandia is considerable: from as little as three months in the most northern locations to eight months in the south [16,23]. This wide range is important for our conclusions. In particular, if one restricts attention to breeding season lengths of six to eight months, the results in figure 3 predict little variation in the dynamics. This is consistent with the results of Dalkvist *et al*. [22]. They studied data on Fennoscandian voles with breeding season lengths of six to eight months, and found that the difference between these season lengths was not a statistically significant factor in differences in cyclicity.

To explore our results further, we have plotted parameter planes in which the generalist predation level varies on one axis, with the breeding season length on the other. The north–south gradient in Fennoscandia corresponds to a curve in this parameter plane. There is currently insufficient data to enable a detailed formulation of this curve, but it is likely to be nonlinear, especially when factors such as landscape fragmentation [22,35,36] are taken into account. Our results suggest that more information on the form of this nonlinear curve would be a key step in understanding the complexities in vole dynamics, such as the sudden change from annual to multi-year cycles that occurs in mid-Fennoscandia [7].

In the model, we forced both the prey and predator growth rates with the same forcing term. Although weasels are able to breed throughout the whole year, this happens only in peak years of vole populations [37]. Therefore, including a weasel breeding season is valid, but it is unclear whether this should follow the same pattern as the vole breeding season, on which the forcing term was based. In order to understand the effects of different forcing in the predator growth rate, we considered the scenario of no forcing for the weasels (see the electronic supplementary material). This yielded similar results, showing the complexity of moving on a north–south gradient, the change to annual cycles as generalist levels increase and the possibility of switching between different multi-year cycles. However, the pattern of decreasing breeding season length leading to increased period was less clearly defined.

Significant changes in vole population dynamics have been observed in recent years in Fennoscandia, with the boundary between annual and multi-year cycles moving further north, and very few locations still showing 5 year cycles [24]. These changes are widely attributed to climate change. This is affecting the Fennoscandia system in many different ways and the implications for vole dynamics are consequently complex. However, two effects of climate change are increases in both breeding season length and generalist predation level [38,39]. Our results (figure 3) predict that these two effects could together cause decreases in cycle length and a northerly movement of the location of the switch between annual and multi-year cycles. Because annual cycles typically have a much lower amplitude than multi-year solutions (figure 2), and because voles are a keystone species in Fennoscandia [39], a switch to annual cycles is highly significant for the whole ecosystem.

Fennoscandia is not the only region in which there are established geographical variations in vole cyclicity. For example, vole cycles occur over large regions in central Europe [40] and in Hokkaido, Japan [10]. Notably, Stenseth *et al*. [10] showed the importance of season length in driving the gradient of cyclicity in Hokkaido. In central Europe, both the period and amplitude of vole cycles increase from north to south [40]; this is in the opposite direction to the trend in Fennoscandia. Data on breeding season length and level of generalist predation are currently very limited for this region, and our work argues for the importance of such data in understanding this trend.

Large-scale geographical variations in population dynamics also occur in other populations: Canadian muskrats are another classic example [11]. Our results argue for a detailed appraisal of the role of breeding season variability on the population dynamics in such cases. Although we expect similar trends to those found in the present study, they may be quenched by the effects of seasonality in parameters other than simply the birth rate; specific case studies are essential.

Seasonality varies geographically owing to a multitude of factors including latitude, altitude, proximity to a coastline and general weather patterns. Our work argues that this cannot be ignored as variations in seasonality can be an important driver of the observed population dynamics.

- Received November 15, 2012.
- Accepted December 20, 2012.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.