{'fraction_non_alphanumeric': 0.041160818227174455, 'fraction_numerical': 0.013354398802594378, 'mean_word_length': 4.672735849056604, 'pattern_counts': {'":': 0, '<': 22, '<?xml version=': 0, '>': 22, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Indirect effects of primary prey population dynamics on alternative prey

We develop a theory of generalist predation showing how alternative prey species are affected by changes in both mean abundance and variability (coefficient of variation) of their predator's primary prey. The theory is motivated by the indirect effects of cyclic rodent populations on ground-breeding birds, and developed through progressive analytic simplifications of an empirically-based model. It applies nonetheless to many other systems where primary prey have fast life-histories and can become locally superabundant, which facilitates impact on alternative prey species. In contrast to classic apparent competition theory based on symmetric interactions, our results suggest that predator effects on alternative prey should generally decrease with mean primary prey abundance, and increase with primary prey variability (low to high CV) - unless predators have strong aggregative responses, in which case these results can be reversed. Approximations of models including predator dynamics (general numerical response with possible delays) confirm these results but further suggest that negative temporal correlation between predator and primary prey is harmful to alternative prey. We find in general that predator numerical responses are crucial to predict the response of ecosystems to changes in key prey species exhibiting outbreaks, and extend the apparent competition/mutualism theory to asymmetric interactions.


Introduction
Most predators switch to alternative prey species when their primary, favourite prey becomes scarce. The classical prey model of optimal foraging theory reflects such a switch, showing how the diet breadth of consumers increases when consumers' favourite food items are at low density (Emlen, 1966;Schoener, 1971;Pulliam, 1974;Charnov, 1976;Pyke et al., 1977;Pyke, 1984;Stephens and Krebs, 1986), and shrinks again at high favourite food density. For community dynamics, this implies that alternative prey species are influenced by the dynamics of the primary (i.e. preferred) prey through their influence on the predator's functional (Oaten and Murdoch, 1975) and numerical (e.g. Holt and Kotler, 1987;Wilson and Bromley, 2001) response to shifts in primary prey abundance. Predator numerical responses can be aggregative through movements (Turchin and Hanski, 1997), or demographic through reproductive changes; in most cases a mixture of both. There is, however, an essential tension between the consequences of functional and aggregative/numerical responses for alternative prey. When the primary prey increases in numbers, individual predators tend to eat less of alternative prey species, but predators also tend to be more numerous, which potentially increases the overall impact of predators on alternative prey populations. The impact of predators on alternative prey results from a balance between the numerical and functional components of the predator response, which determines whether the primary prey has a positive or negative impact on the alternative prey. Apparent interactions can therefore take the form of competition or mutualism (Abrams et al., 1998;Bêty et al., 2002;Brassil and Abrams, 2004), or even amensalism versus commensalism, if the primary prey has greater impact on the secondary prey than the reverse, which often seems to be the case because of differences in overall biomass and densities (Sinclair, 2003). Here, we will consider these asymmetric interactions, i.e. prey 1 indirectly affects prey 2 but not the other way around. It seems that such asymmetric subsets of the food web (trophic modules) are not only quite widespread (Stouffer et al., 2012), but also important to consider from a functional or conservation perspective (DeCesare et al., 2010;Wittmer et al., 2012).
Predators generally prefer to feed on fast-reproducing species, r-strategists that invest heavily on reproduction and less on survival, if only because these prey are easier to catch and have large maximal densities. Such fast-reproducing species are strongly influenced by environmental variability, and therefore tend to have dynamics that are both highly variable (Saether et al., 2002;Sinclair, 2003) and often nonstationary (e.g. Angerbjörn et al., 2001). The influence of variability and nonstationarity of primary prey dynamics on alternative predation is not yet well developed in the otherwise abundant apparent competition/mutualism literature (Holt and Lawton, 1994;Abrams et al., 1998;Brassil and Abrams, 2004;Brassil, 2006;Schmidt and Ostfeld, 2008), despite its relevance to a number of species of scientific and conservation importance (see Schmidt and Ostfeld, 2008 for examples). In this paper we investigate how the interplay between the mean and variability of primary prey abundance affect alternative prey species demography. We restrict the definition of alternative prey in two ways in this article. First, it is a species that -unlike the primary prey -is either not available or not nutritious enough for the predator to specialize on it year-round. Second, which might be a corollary of the first, the alternative prey cannot drive the aggregative/numerical response of the predator, which is mostly influenced by the primary prey (e.g. New et al., 2011New et al., , 2012. These restrictions make the asymmetric nature of the interaction all the more likely, and correspond well to the examples described below. A typical terrestrial example of highly variable primary prey are rodents such as voles and lemmings (subfamily Arvicolinae), that exhibit large-amplitude cycles, especially in boreal and arctic regions. These cycles are often non-stationary, because rodent vital rates react to trends in climatic variables (Kausrud et al., 2008). As Lack (1946) remarked early on, rodents are preyed upon by an important guild of avian and mammalian predators. This guild includes foxes and mustelids, on the mammalian side, and on the avian side, raptors and some other birds such as skuas and corvids. All these rodent-eating predators have the habit of switching to alternative food sources when the rodent cycle is at a trough; often such alternative prey consists of eggs and juveniles of gamebirds, waterfowl, etc. (see Valkama et al., 2005, for more details). Thus, the breeding success of many bird species is severely impacted by the population dynamics of rodents (Lack, 1946;Summers et al., 1998;Wilson and Bromley, 2001;Blomqvist et al., 2002;Bêty et al., 2002;Valkama et al., 2005;Schmidt and Ostfeld, 2008). It has even been hypothetized that bird breeding habitat in the Arctic is partly determined by association with cylic rodents (Blomqvist et al., 2002;Gilg and Yoccoz, 2010).
Taking a larger view of the food web, such variable primary prey are often key species in their ecosystems, with a large number of links to other species (Sinclair, 2003;Jordán, 2009;Krebs, 2011). Examples of key prey whose dynamics can influence alternative prey species range from lemmings and hares (Lepus sp.) in boreal and arctic landscapes (Krebs, 2011) to wildebeest (Connochaetes sp.) in the Serengeti (Sinclair, 2003); or in marine ecosystems, from sandeels (family Ammodytidae; Matthiopoulos et al., 2008) to small, overabundant pelagic fishes that create so-called "wasp-waist" ecosystems (Cury et al., 2000).
Though such species are sometimes referred to as "keystone" (e.g. Cornulier et al., 2013), this is a misnomer, since their importance is largely due to their large maximal biomass at peak densities (unlike keystone species which have a disproportionate effect on the ecosystem per unit of biomass, Power et al., 1996). Hence, we refer to rodents and their r-strategist counterparts in other ecosystems simply as "key" prey species.
In a previous paper, Matthiopoulos et al. (2007) developed a model for predation by harriers (Circus cyaneus) on red grouse chicks (Lagopus l. scoticus, an alternative prey to voles Microtus agrestis) in Scottish moors, for different levels of abundance of voles and pipits (Anthus pratensis), both of which are important in the harriers' diet. This was done using empirically measured multispecies functional and aggregative responses. The study was largely motivated by management of a human -wildlife conflict; hunters typically want to shoot more grouse and see less harriers, while conservationists wish to protect the hen harrier, an endangered raptor in the UK (Thirgood and Redpath, 2008). The impact of harriers on grouse is managed through a number of techniques, that range from the illegal killing of raptors to diversionary feeding (Redpath et al., 2001). A modelling assessment of the various management techniques is given in New et al. (2012). We use this system as a key empirical example to motivate general theory rather than the object of study per se. Red grouse corresponds well to our abovementioned definition of alternative prey.
Field vole numbers -the primary prey of harriers -were assumed for simplicity to be a constant in Matthiopoulos et al. (2007). However, vole abundance can vary greatly from year to year, and harriers react numerically to these variations (Redpath et al., 2002). The model we develop here relaxes the assumption of constant vole primary prey availability. The way we represent primary prey dynamics is akin to a resource pulse (Holt, 2008;Schmidt and Ostfeld, 2008). Through a series of progressive approximations of the detailed, empirically-based model of vole-harrier-grouse dynamics, we formulate a simplified model that can apply to a number of other prey communities that share the same basic features with the rodent-raptor-gamebird system, which are: (1) strongly fluctuating primary prey and (2) asymmetric interactions.
We initially extend the vole-harrier-grouse model to account for vole variability, and show how grouse average numbers critically depend on the temporal average of grouse breeding success and therefore, temporal average of the total response (total number of grouse chicks killed = aggregative × functional response). This then leads us to a more detailed examination of how the average total response (number of grouse chicks killed) depends on primary prey variability (i.e. coefficient of variation), and how the effect of primary prey variability interacts with that of mean primary prey density. As previously remarked by Schmidt and Ostfeld (2008), variability in primary prey density can affect a predator's average consumption, because their total response is a nonlinear function of primary prey abundance (a consequence of Jensen's inequality). We show that the direction of such effects is driven by mean primary prey abundance. We first consider the case of fast aggregative responses to primary prey density (section 2 and 3), but also extend some of these results to delayed numerical responses (section 4).

The model
The model we present here is derived from Matthiopoulos et al. (2007), with the addition of fluctuating dynamics for voles. In Matthiopoulos et al. (2007), a third prey, pipits, was considered. However, for simplicitly, in most of this paper, we fix the density of the third prey, pipits, to a constant. This assumption (relaxed in Appendix A) is justified because pipit numbers vary usually less than vole numbers.

Building blocks 2.2.1. The vole model
Throughout the manuscript, we use the letters V to denote vole or more generally primary prey density, and G for grouse or alternative prey. We use two alternative descriptions of vole population dynamics, a deterministic model producing limit cycles (based on Maynard Smith and Slatkin, 1973) or a simple probability distribution of values (either lognormal or gamma, because abundance distributions are non-negative and skewed to the left). These contrasting formulations are justified biologically because observed temporal fluctuations in voles range from deterministic-looking cycles to lognormal noise driven by abiotic variables (Stenseth, 1999), and mathematically because the periodicity of vole dynamics does not influence the main results which are driven by the mean and coefficient of variation.
The Maynard-Smith and Slatkin (MSS) model is given by eq. 1, in which R V is the maximal vole growth rate, K V is the threshold density marking the onset of densitydependence, k is the time lag of density-dependence (one or two timesteps), and γ is an exponent characterizing the abruptness of density dependence.
Alternative dynamical models such as autoregressive models on a linear scale (Royama, 1992) provide similar results for the predation models considered below, so long as we consider aggregative predator responses.

The grouse model and harrier predation
The grouse model is based on Matthiopoulos et al. (2007), and is essentially a delayed version of the Hassell model for population cycles, including some predation on the maximal birth rate (eq. 2 below) Here, S is the grouse survival rate, K G is a threshold for density-dependence, and β controls the shape of density-dependence. B(G t , V t ) is the birth rate, modelled as a function of both prey densities H t is harrier density and f (B 0 G t , V t ) is the raptor's functional response (equation 5) to the densities of grouse chicks and voles (see Table 1). Assumptions about H t are either H t = H 0 (a constant) or H t = α 1 V t + α 2 P where P is pipit density and α 1 , α 2 relate vole and pipit density to H t (this is a linear aggregative response, as in Matthiopoulos et al., 2007). Later, in section 4, we consider population dynamics for H t . Note that we can also write the model as

Harrier functional response
The functional response relates the intake rate of harrier pairs to the numbers of voles, grouse young, and pipits (considered constant here). Assuming the density of the i th prey species is denoted by N i , the functional response on species i can be written where parameters a i , m i , h i have been estimated empirically in Smout et al. (2010), with values given in Table 1 (for consistency with Smout et al., 2010, i = 1 for grouse, 2 for pipit, 3 for vole). We use this specific parametrization of the model so as to consider a biologically relevant scenario (illustrated in Fig. 1); more general functional forms are considered later on. , while in C we see the functional response on grouse (prey 2) given vole density (prey 1). Three kinds of aggregative responses (D)constant (plain line), linear (dashed) and sigmoid (dash-dotted) -eventually lead to three kinds of total response in E (total intake on prey 2 as a function of prey 1= aggregative response to prey 1 × functional response on prey 2 as a function of prey 1). All functional response parameters are identical to those estimated in Smout et al. (2010).

Effect of vole variability on the full community model
Vole population dynamics (Fig. 2) modify grouse population cycles through harrier predation, though not to a level generating synchrony between vole and grouse  Fig. 2, as well as more simplifed probabilistic representations (lognormal noise), we found that the temporal average of grouse population size, the temporal average of per capita grouse breeding productivity, and the temporal average of the total response are all strongly correlated (Fig. 3), and respond similarly to variations in mean and variability (i.e. coefficient of variation) in vole density. In Appendix C, we show why such correlations are expected in models for regulated (albeit fluctuating) populations. Therefore, in the following section, we focus only on the temporal average of the total response, which is more general and easier to compute.  Table 1 for other parameters. , and average number of chicks predated (C, i.e. average total response). Vole population dynamics is given by a simple lognormal probability distribution, whose mean and CV are varied in the x and y axes. Parameters are otherwise equal to those in Table 1 (with empirically measured aggregative response instead of constant predator numbers).

Average total response on alternative prey: effect of mean density and variability in primary prey
In the preceding section, we showed how the total response, averaged over the years, provides a good indicator of predator impact on secondary prey, when the predator aggregates to the primary prey. This will be the case in situations where predators mostly impact the juvenile component of the alternative prey population (see Appendix C). Here, we therefore analyze the temporal averageT of the total response T (V t ), In the long run, this value converges to its expectation over realizations, given a probability distribution ρ(v) for voles densities. Because of the ergodic theorem, this is true irrespective of possible autocorrelation in the V i values and leads to 3.1. General formulation using gamma functions Using Taylor approximations, it is possible to derive insights on how the average total response is affected by the mean and variability of primary prey density (this works best for CV < 0.3). However, because the primary prey has usually very variable dynamics (i.e. CV ≥ 1), it is worthwhile to derive results for "large noise" conditions. Using relatively mild assumptions for the functional and aggregative responses, it is possible to derive such results. A typical form for the functional response on prey 1 as a function of prey 2 could be f (V ) = f 0 e −λV (grouse chicks killed as a function of rodent density V ), and for the aggregative response, some power function H(V ) = αV l (number of predators purely driven by the primary prey density). More complicated functions are possible, but do not change the main argument. In general, the functional response f depends on both G and V ; the reason why we approximate f (G, V ) by f (V ) is because of the stronger dependency on V (Fig. 1), especially when the grouse population is regulated. Similar functions are presented in Schmidt and Ostfeld (2008), and the Appendix B describes how a nearly exponential decrease of the functional response can be obtained even in cases of optimal predator foraging, where we would usually expect a function with both concave and convex parts.
This combination of functional forms yields a total response T (V ) of the kind A l V l e −λV (Fig. 4A, with the exponent l changing how the total response accelerates at low densities, and a total response maximum at l/λ). The total response, for efficient predators that forage optimally, might drop faster than exponentially, e.g. steeper declines in the functional response, which complicates the mathematics but does not change the qualitative results. Because it will be important to the computations, let us mention that AV l e −λV has the same functional form as the probability distribution of a gamma random variable (presented in Appendix D). Then, using a gamma random variable to model the probability distribution of primary prey values, we reduce the problem of computing the average total response to computing particular integrals of gamma functions (Appendix D).
We found that in the case l = 1 (linear aggregative response H(V ) = αV ), increasing the mean rodent density m first increases predation on prey 2, up to a maximum after which it starts decreasing (Appendix D and Fig. 4).
This means that the deterministic total response curve T (V ) is a good indicator of the effect of the average rodent density (denoted here m) on E(T (V )). The relationship may seem intuitive -the effect of the mean on the average total response is the same as that of that of the deterministic vole density on the deterministic total response -but verification is important for two reasons. First, there is substantial variability around the mean (e.g. CV = 1, see Fig. 4), and second, the variability is asymmetrically distributed (due to the gamma distribution, which corresponds well to the shape of real abundance values).
In contrast to the case with linear aggregative response, when l = 0 (contant predator numbers, H(V ) = H 0 ), we have a constant decrease of average predationT with mean rodent density m, again predictable from the deterministic T (V ) curve.
The effect of variability is not so intuitive, though the curvature of the total response gives some indication. For l = 1 (H(V ) = αV ), increasing variability slightly increases predation for low primary prey mean m, then decreases predation for intermediate m values, and for very large m, increasing CV increasesT (Appendix D for a general proof, similar results for l = 5 in Fig. 4). The latter case of large m is unlikely because it would correspond to overabundant key herbivores all the time. For a large range of mean values m, in the presence of a strong aggregative response, more primary prey variability will therefore decrease predation. In contrast, for l = 0 (H(V ) = H 0 ), increased CV always generates increased predation (Appendix D). Table 2

Vole-Grouse-Harrier example: numerical integration
In section 3.1, we used theoretical models without measured parameters. In this section, we present an example with the empirically measured parameters for the voleharrier-grouse system. Fig. 5 shows how increasing pipit density (which is assumed constant over time) makes the aggregative response switch from a linear function of vole density to a near-constant, and modifies the effect of vole variability on grouse. In all cases, however, it appears that the total response is mostly decreasing and convex, and therefore increased mean rodent density decreases predation on grouse while increased rodent density CV increases predation.  Table 1.

Including delays in the predator response: Taylor approximation and predator-prey covariance
Delays in the predators' numerical response can arise either in the aggregative response (due to slow movements) or in the predators demography (due to the delayed effects of food on reproductive success). Delays in the aggregative response can happen when the environment is characterized by strong seasonality, a situation covered extensively in Appendix E that yields results qualitatively similar to those presented below. The consequences of a delay in the predators' reproductive response are difficult to study fully with a parametric model when predation is important, because of feedbacks between predator and prey populations. We plan to perform detailed numerical simulations in a future publication. However, we consider here approximations that take into account the possible effect of the predator-prey covariance generated by delays in the predator numerical or aggregative responses, assuming predation pressure is relatively low, so that prey species 1 and 2 cannot be considered to be fully synchronized (in which case the average total response is not correlated to average population size of the alternative prey). For small predation pressure then, we can focus on the average total responseT to find except now the predator density H i has its own dynamics, correlated to some extent with V i . Simplications are obtained by remarking that We can derive results using classical small-noise approximations to nonlinear func- (Hilborn and Mangel, 1997). Using the previous formula f (V ) = f 0 exp(−λV ), and dropping indexes, we obtain .
We could also make use of a first order approximation for the functional response ) for small deviations (using interchangeably f (V ) and f (E(V )) to simplify notations). The first order approximation helps simplifying the covariance Cov H). This finally leads to the approximate formula for the average number of grouse killed (=average total response) In this formula, we see • As before, increase in primary prey variability, for a relatively constant predator number, can increase the predation rate • A positive covariance between predator and primary prey density will decrease the average number of alternative prey killed. A positive correlation in numbers -as happens in the aggregative response case -implies that when the individual kill rate on alternative prey is high (e.g. low rodent density), there are also fewer predators.
• In contrast, a negative covariance between predator and primary prey density does not benefit the alternative prey. This might typically happen with a delayed predator numerical response to the primary prey (with some degree of phase opposition in the cycles).
The results above are corroborated by those for delayed aggregative responses, in Appendix E.

Discussion
We modelled the effect of generalist predation on alternative prey as a function of the mean and variability of primary prey numbers. Our work was motivated by an empirically-driven model for the interaction between voles and grouse through harrier predation, and recents reports of changes in the mean and sometimes variability of vole or lemming numbers (Ims et al., 2008;Schmidt et al., 2012;Cornulier et al., 2013). However, other key herbivores (e.g., hares, wildebeest), with large maximal growth rates compared to other prey species and large maximal densities, have rather variable and non-stationary dynamics, and have the capacity to indirectly influence the population dynamics of alternative prey species.
We showed that the temporal average of the total response (the total number of alternative prey killed each year) is a good index of predator impact for this model. It was indeed strongly correlated with both alternative prey breeding productivity and average abundance. Our results apply to a large set of asymmetric apparent interaction modules, where prey 1 influences prey 2 but not the other way around.
low V + -Note: First column: By Cov(H,V) we denote the covariance (at temporal lag zero) between predator and primary prey. Cov(H,V)>0 corresponds to the aggregative response scenario, Cov(H,V)=0 to constant or randomly varying predator scenario and Cov(H,V)<0 to a numerical response scenario (negative correlation is generated by a short demographic or movement delay). All of our results assume Cov(G,V) ≈ 0 (unsynchronised G and V). Situations with synchronised G and V, corresponding to high predation pressure on alternative prey, would warrant further analyses (see Appendix F). The three other columns show, from left to right (1) the phase of the cycle during which the alternative prey suffers the most predation, i.e. the cycle phase in which management actions would be best implemented to protect the alternative prey (2) the effect of increasing the temporal mean of the primary prey density on the alternative prey (3) the effect of increasing primary prey variability, as measured by the coefficient of variation (CV). * If the mean primary prey density is very low (very skewed distribution), intermediate V values from the prey perspective might not represent intermediate values in the total response of predators. * * Negative effect possible at very large mean primary prey density, which is unlikely (it would require a permanent overabundance of key herbivores).

Effect on alternative prey of mean and variability in primary prey abundance
The total response on prey 2 as a function of prey 1 density is the functional response times the aggregative response. Here, the term functional response differs from its classical use, i.e. it is the per predator kill rate of prey 2 as a function of the density of prey 1. Such a functional response is derived from the multispecies functional response, as opposed to Holling's classic single-species functional responses. The problem of evaluating the effect of primary prey abundance and variability on alternative predation was simplified by considering only variation in the primary prey, given it is the prey that dominates both the functional and aggregative responses. Therefore, in our model, the total response for prey 2 depended only on prey 1.
The key quantity in the model was the temporal average of the total response. The total response is generally a nonlinear function, hence its temporal average obeys Jensen's inequality, E(T (V )) > T (E(V )) if T is convex, which generates possible effects of prey 1 variability. Usually, approximations of nonlinear functions are used when such averaging is needed. However, the primary prey (e.g. rodents) in many real systems exhibits considerable variability in abundance (i.e. CV ≥ 1) which implies that results from linearized models may have limited applicability. We circumvented this problem using a more general approach, specifying functional forms that allowed us to represent primary prey interannual variation as gamma distributions. Mathematical and computational results using different probability distributions for primary prey variability agree, and we summarize these results below and in Table 2.
We found that when predators are in relatively high and constant numbers, the overall kill rates of alternative prey decline almost exponentially with respect to primary prey density (such in the vole-harrier-grouse case when pipits are numerous, Fig. 5). Therefore, predation on alternative prey decreases when mean primary prey increases, but predation increases when primary prey variability increases.
With respect to the effect of the mean, it is an instance of commensalism (0/+) between prey, and our finding add to those of Abrams et al. (1998) who found the possibility of mutualism (+/+) in interactive systems. With respect to the effect of variability, Schmidt and Ostfeld (2008) suggested a similar outcome based on Jensen's inequality, i.e. E(T (V )) > T (E(V )) because T is convex. In that case, primary prey variability is inherently detrimental to alternative prey survival.
However, the literature shows cases where predators can be sparse or even absent when primary prey is low (Gilg et al., 2006). In these cases, our theoretical modeling shows that such strong predator aggregative responses may generate the opposite pattern: 1) Increases in mean primary prey are detrimental to secondary prey but 2) temporal variability in primary prey can decrease predation on secondary prey. The occurrence of such cases might be restricted to high-arctic latitudes however, and calls for more empirical research.
More generally, our results suggest that empirical measurements of predators' aggregative responses to primary prey are at least as important as measurements of their functional responses when predicting predation on alternative prey species.

Limitations and extensions 5.2.1. Can variability and mean really be separated?
So far we equated primary prey variability to its interannual coefficient of variation (CV). This choice needs to be discussed, especially in the light of changes in small mammal cycles currently being observed (e.g. Cornulier et al., 2013). First, our decomposition separates the effect of the mean primary prey density/abundance from the effect of the CV. We chose CV because it is well-known that real data of population abundance show a scaling between the mean and variance of abundance. In the absence of demographic stochasticity, Taylor's law suggests a log-variance/log-mean regression with a coefficient of 2, which is equivalent to a CV independent of the mean. Thus, over short enough timescales for the dynamics to be approximated by a stationary process, and for large populations in which demographic stochasticity has limited influence, we can assume independence of mean and CV.
However, representing changes in variability through the CV does not solve all problems when considering nonstationary dynamics. Recent research on cyclic populations, for example, suggest that mean and CV covary to some extent Barraquand and Husek, 2014). Changes in vole dynamics in e.g. the Kielder forest suggest predominantly decreases in the mean spring vole abundance Cornulier et al., 2013. If changes in primary prey dynamics are primarily driven by a decrease in the mean, cycle "loss" (absence of peak years) will be associated to more predation on alternative prey. If changes are mostly a decrease in CV, then cycle "loss" will usually be associated in contrast to less impact on alternative prey (unless there are strong predator aggregative responses).
Changes in the temporal dynamics of key herbivores might be represented in the mean-CV plane . This might help deciphering the possible covariation of mean and variability over long timescales, for nonstationary dynamics, and predict the effect of changes in (cyclic) primary prey dynamics.

Primary prey temporal autocorrelation
The results we derived here, assuming mostly simple predator aggregative responses, are valid irrespective of temporal autocorrelation in the primary prey abundance. This follows from Birkoff's ergodic theorem, the generalisation of the law of large numbers for autocorrelated sequences. For example, even though rodent cycles have a famous autocorrelation structure (Turchin, 2003), Birkoff's ergodic theorem implies that the temporal average of a nonlinear function of rodent densities will be equal to its expected value (average over realizations), calculated according to the marginal probability distribution of rodent densities. For another use in ecology, see Benaïm and Schreiber (2009); Barraquand and Yoccoz (2013).
For the simple model with aggregative predator response, use of Birkoff's ergodic theorem is correct, but our ability to ignore the autocorrelation critically depends on how we relate predator numbers to primary prey densities. When delays are added, predator numbers depend on past prey densities, and the autocorrelation structure between primary prey densities and predator numbers determines the average predation rate. Seasonality can generate similar temporal patterns (Appendix E). This can be dealt with qualitatively using simple approximations based on Taylor expansions. Eq. 11 clearly shows that positive temporal correlation between predator and prey suppresses the kill rate on alternative prey, while negative temporal correlation increases it. Thus, longer predator reproductive delays increase alternative predation. But such approximations can break down when primary prey variability is large. To a certain extent, delays might be present in the vole-grouse-harrier example, because predators settle as a function of current vole density (aggregative response) but have a tendency to stay longer (New et al., 2011(New et al., , 2012. This happens because of other alternative prey that sustain predators in the absence of their primary prey (as pipits do for harriers), and it can take time to deplete these other food sources. We do not think this should alter the qualitative results in the vole-harrier-grouse case, but are currently considering the consequences of such scenarios, that tend to favour synchrony between prey species when predator abundance is large (Barraquand et al. unpublished data).

Synchronization of alternative prey dynamics
We assume throughout the paper that although breeding success and numbers of prey 2 are affected by predation, often resulting in a cyclic breeding success, predation is not enough to fully synchronize prey 2 abundance with that of prey 1. This is correct for the empirical values we possess for vole-harrier-grouse in Scotland and many other bird-rodent systems (e.g. Blomqvist et al., 2002), but might not be true everywhere. In contrast, one hypothesis to explain relatively short cycles of Fennoscandian grouse is predator sharing with cyclic rodents (the Alternative Predation Hypothesis, Angelstam et al., 1984). Larger numbers of predators than what we considered here can lead, in the harrier predation model, to situations where the vole and grouse cycles synchronize (with a one-year timelag, Appendix F). Grouse periodicity then matches vole periodicity, and this frequency entrainment is stronger when the vole cycle exhibits strong asymmetry (as with the MSS model). In the harrier case, densities of raptors approximately above 1 or 2 pairs per km 2 are needed to generate such high predation pressure, which might be difficult to obtain at large spatial scales in most habitats. However, recent empirical work suggests that communities of rodents are for instance synchronized by the conjugated action of the whole predator guild Carslake et al., 2011), which makes large predator densities feasible.

More feedbacks in prey 1 -predator -prey 2 modules
The longer snowshoe hare cycle, which is thought to be generated by predators (Krebs et al., 2001) provides yet another nice empirical motivation. Numerous predators such as lynx and coyotes react demographically to snowshoe hare densities (O'Donoghue et al., 1997(O'Donoghue et al., , 1998a and inflict losses to alternative prey (e.g. squirrels, sheep, O'Donoghue et al., 1998b;Arthur and Prugh, 2010). Hence it would be worthwhile to continue the work on asymmetric apparent interaction modules to investigate synchronizing effects, now modeling feedbacks between hares and predators (see King and Schaffer, 2001), which should allow large numbers of predators to influence the dynamics of prey 2 when prey 1 is low. This situation can be represented prey 1↔ predator → prey 2. This would be a logical extension to the models presented here -preliminary simulations suggest that such delayed predator responses favour synchronization between prey species.
Another logical extension, to investigate apparent interactions and synchronization between similar prey species (as in Korpimäki et al., 2005;Carslake et al., 2011), would be a different model structure, with also feedbacks between prey 2 and predator dynamics (prey 1↔ predator ↔ prey 2). Abrams et al. (1998) considered such a model with predator-generated cycles, and they show that mutualism is possible between prey species in the cyclic regime, which goes well with our findings (increased mean prey abundance and decreased prey variability are beneficial to alternative prey). However, Hanski and Henttonen (1996), again with a differential equation model, inspired by interactions between rodents and mustelids, show that apparent competition can mediate competitive coexistence, and that a rodent species not likely to cycle can be entrained by a competitor more susceptible to predation. Induced cyclicity is then helping coexistence, but the presence of prey 1 is hindering coexistence. In the case where the alternative prey can also influence predator dynamics, our results therefore will not always hold, so this warrants further research.

Conclusion
Under the assumption of both asymmetric interactions (prey 1→ predator → prey 2) and absence of synchronization of prey abundances, increased mean abundance of primary prey will usually help alternative prey species, and increased variability of primary prey will be detrimental to alternative prey species. Unless predator aggregative response is strong, in which case both results can actually be reversed (Table 2). Northern rodents such as voles and lemmings, that are primary prey in many such apparent interaction modules with ground-nesting birds, have both large maximal abundances and large variability (CV close to 1). Their presence often protects birds, as increases in the mean primary prey density decreases predation on alternative prey, although in a suboptimal way, because a less variable primary prey dynamics would usually offer better protection. It is likely that such suboptimal protection by primary prey is felt by other alternative prey species in a variety of ecosystems.

Appendices Appendix A. Dealing with variation in alternative prey species
Appendix A.1. Variation in the focal alternative prey So far, our 2-prey total response models are assuming prey 2 density is roughly constant, i.e. E(T (V )) approximates If G is not fixed but fluctuating instead, E(T (V, G)|G) is a random variable. We can remark that in general, with ρ V,G (v, g) the joint probability density of V and G. And in the simple case where they are independent (no synchrony) we arrive at So we just need to integrate once more the average total response, this time with respect to grouse variation. This shows why our simple approximation is working, despite the variation in grouse density in the numerical models.
We might even take into account synchrony between grouse and vole using some functional form G = G 0 +βV . That is a deterministic link expressing an extraordinarily strong level of synchrony. A more realistic model could assume that grouse is linearly related to vole density on average, e.g. G|V ∼ N(G 0 + βV, σ 2 G|V ), assuming here a normal probability distribution conditional to vole density (we recover independence between V and G when β = 0). This would result in the conditional p.d.f. ρ G|V =v (g; v) being known, and then we can write where all the elements are known, or again, where the quantity between parentheses is the expectation of the functional response on the grouse as a function of grouse density, conditional on vole density. That expectation is possible to compute, at least by numerical integration. Note: Similar decompositions using conditional expectations could help dealing with parameter uncertainty. A difficulty not addressed here is that synchrony between V is G usually working with a lag (Appendix F), thus the arguments above might need to apply to some function of G.

Appendix A.2. Variation in another prey species
The above arguments can be extended for variation in other prey species than the primary prey or the alternative prey of interest. Thinking of the harrier example, let us consider a density of pipits P t , that varies over time, independently of V t . Assuming G t is constant (e.g. applying the arguments of the previous subsection), the total response writes Due to independence, we arrive at which is then equal to E(T (V, P )) =ˆ∞ 0 E(T (V, P |P = p))ρ P (p)dp (A.9) We therefore show that considering pipit density, provided marginal probability distributions of prey densities can be considered independent, only requires one additional integration over the whole pipit probability distribution.
Appendix B. Inclusion of optimal foraging behaviour in the functional response Here we modified a classic type-II multispecies functional response model to allow for optimal foraging -this is generally thought to stabilize such models and enlarge the spectrum for prey species coexistence (Krivan, 1996;Krivan and Sikder, 1999). We follow closely Fryxell and Lundberg (1998, chap. 2) and the initial type II response, similar to eq. (5) in the main text, comes from Hanski and Henttonen (1996), hereafter referred to as HH96.
The basic model of foraging theory, i.e. the diet choice model (Stephens and Krebs, 1986), predicts that prey items should be ranked according to their relative profitabilities e i /h i (energetic content/handling time) and the least profitable prey is ignored until the most profitable type drops below some critical density threshold. For convenience here we assume the most profitable prey type is 1, i.e. e 1 /h 1 > e 2 /h 2 .
In comparison, the type II functional response (sometimes called "shared predation" response, Norrdahl and Korpimäki, 2000) assumes a heavier predation on the less profitable prey when the main prey is abundant. Let us assume that γ i is the probability of accepting the prey item i upon encounter and modify the functional response. Using an encounter (discovery) rate parameter a, the functional response (predator intake rate) now reads classically where V 1 and V 2 are the densities of the prey one and two, respectively. We have γ 1 = 1 and γ 2 = 1 if V 1 > η 1 , 0 if V 1 < η 1 . In the classic theory Stephens and Krebs (1986); Fryxell and Lundberg (1998); Krivan (1996), the long-term individual intake is written Through the derivation of w w.r.t. γ 1 and γ 2 , we find the critical main prey density threshold is η 1 = e 2 a(e 1 h 2 −e 2 h 1 ) Stephens and Krebs (1986); Fryxell and Lundberg (1998); Krivan (1996). This step-like formulation is not practical and is not likely to be observed in nature (partial preferences, Berec and Křivan, 2000) so that Fryxell and Lundberg (1998, chap. 2) propose a sigmoid function instead, γ 2 = exp(zη 1 ) exp(zV 1 )+exp(zη 1 ) (z is a smoothing parameter controlling the steepness of the sigmoid). We therefore included optimal diet choice, at the expense of two additional parameters only 1. The ratio k of energetic content of two prey items 2. The steepness z of the partial preference curve (γ 2 = f (V 1 )) The results in the figure show that the inclusion of optimal foraging tend to produce exponential-like functional response on prey 2 as a function of prey 1 density (bottom right panel).

Appendix C. Correlations between average population size and average breeding success / natural mortality in alternative prey
Let us assume the alternative prey has population dynamics of the form G t+1 = G t R(V t )/k(G t ), with a maximum per capita growth rate R(V t ) = B(V t ) + D(V t ) (either maximum breeding success or natural mortality can be affected by primary prey density) and the form of regulation k(G t ) = e λGt (Ricker), or (1 + G t /K) β (Hassell), or (1+(G t /K) γ ) (Maynard-Smith). By definition, in a regulated population the long-term growth rate has to be zero. The only assumption made here is that primary prey density affects mostly the maximal per capita growth rate, and not the density-dependent structure (this is only an approximation, but makes sense if one thinks density-dependence is mostly set by resources not predators). Illustration of how the probability of attack of alternative prey should vary with main prey density, for an optimal predator preying mainly on rodents (prey 1) and alternatively on gamebird (prey 2). There is a slight threshold (inflection point) but it occurs at quite low rodent densities. Resulting functional response on the right (adaptive predator) and comparison to a type II functional response (HH96, from Hanski and Henttonen (1996)'s apparent competition model). Parameters: D 1 = 5.0, D 2 = 10.0, h = 1/600, a = 1 Dh , k = e 1 /e 2 = 2, z = 0.5 Because ln(G t+1 /G t ) = ln(R(V t ))−ln(k(G t )), we arrive at the equation E ln(R(V )) = E ln(k(G)) (some uses of this technique in Barraquand and Yoccoz, 2013). In the case of the Ricker model, this even transforms into a more direct relationship between population size and the expected log-max pop growth rate E(G) = E ln(R(V )) λ . If either B(V t ) or D(V t ) is close to 0.5 -let us assume D(V t ) is, as in the grouse case -and B(V t ) has small fluctuations, we simplify ln(1+B(V t )−0.5) ≈ B(V t )−0.5, and thus approximately E(G) ∝ E(B(V ) + constant). If we also have E(B(V )) ∝ E(T (V )) (proportionality of the breeding success to the expected total response), which should be the case whenever G does not fluctuate too much, we then obtain E(G) ∝ E(T (V )). This is what we found with the vole-harrier-grouse model. This relationship works in a number of models for k(G t ) (e.g. Ricker, Hassell, Maynard-Smith).
Appendix D. Sensitivity analysis of the average total response to mean and coefficient of variation in primary prey density Gamma distribution reminder. The probability density function is ρ( with (k, β) = (k, 1/θ) and E(X) = kθ, V(X) = kθ 2 . Therefore, the coefficient of variation is c = 1/ √ k, and the shape parameter is k = 1/c 2 .
Appendix D.1. Computation of the general integral for the total response The integration uses the fact that´∞ 0 ρ(x)dx = 1 because ρ is a pdf on [0, +∞) to simplify the integral. Let's call the integralT =´∞ 0 ρ(x)T (x)dx. We havē We therefore have in the general casē Let's consider a simple case Case l = 1, H(V ) = αV : linear aggregative response. Because Γ(k + 1) = kΓ(k), we obtainT With the other formulation for the gamma distribution, denoting λ = 1/η (η is the vole density at which the individual predation rate on the alternative prey drops to e −1 = 37% of its initial value), we obtain NB: β λ+β = η η+θ with β = 1/θ. Or, writing the mean rodent density m = kθ and c = 1/ √ k the coefficient of variation of rodent density (the relevant variable for rodent variability) Or again, since k = 1/c 2 , a more practical formulation for computations is We use only the last version, in the following calculations so we drop the subscript 3.
Case l = 0: constant predator density. We then have H(V ) = H 0 and A = H 0 f 0 in the integral.T Of course, we can have any combination of the two previous cases, such as H(V ) = H 0 + αV ; due to the additivity of the integral this leads tō Below, both cases are interpreted separately, while B11 will be evaluated at a later date.

Appendix D.2. Sensitivity analysis of the averaged total responseT
We computed partial derivatives with respect to the average rodent density m and its coefficient of variation c (through parameter k = 1/c 2 ). Finally, The mean rodent density m has a positive effect on the averaged total responsē T = ϕ(m, k) when m is small, and a negative effect when m is large (maximum effect at m = η for the averaged total response, η being also the maximum the non-averaged total response T (V ) here).
Appendix D.3.2. Case l = 0 (≈ without aggregative response) We need to compute ∂ψ ∂m . Let's write A = H 0 f 0 . We have , we obtain ∂ψ ∂m < 0. So the effect of increasing mean rodent density is to decrease the average kill rate when l = 0 (no aggregative response), and corresponds to the monotonically decreasing total response in that case.
Let's assume that m kη << 1, which seems likely (m < η implies mean rodent density is lower than the threshold at the functional response drops to e −1 = 37% of its maximum, and k > 1). In that case, we have ln(1 + m kη ) ≈ m kη which then leads to So u ′ (k) < 0 iff m > η (assuming the approximation works). We know ∂ϕ ∂c = Amu ′ (k)e u(k) ∂k ∂c and ∂k ∂c = −2 Or, in other words, there is a positive effect of rodent variability on the kill rate of grouse at large average rodent density. The results presented here are actually similar to those that can be obtained using "small-noise" approximations (i.e. with Taylor expansions of nonlinear functions). Now the question that remain is: can there be other 'negative' effects of variability at not only very high, but also very low rodent density? This may occur when the aggregative response is accelerating, i.e. H(V ) = α l V l ⇒ T (V ) = A l V l exp(−λV ). To investigate this question, we need to redo these computations with a power l > 1 in the integral, and this is done in Appendix D.5.
kη . x > 0 so ln(1 + x) > 0 which means u ′ (k) < 0 and ∂ψ ∂k = ∂ ∂k e u(k) = u ′ (k)e u(k) < 0. Consequently, ∂ψ ∂c > 0, and increase in rodent variability, for a constant mean, will increase the kill rate on the grouse. We recover the negative effects on grouse of rodent variability that can be found using simply Jensen's inequality, in the case of constant predator numbers.
Appendix D.5. Effect of variability, case l > 1, accelerating total response at low rodent densities Here we show that positive effects of vole variability on grouse kill rates are unlikely at very low mean rodent density, despite the convex shape of the total response curve (l > 1) suggesting so. We rewrite the general eq. D.4 for the averaged total response when assuming the flexible v l e −λv form: we have the relations β = 1 θ = k m , λ = 1 η . Hence, using also Π(k, l) = Γ(k+l) Γ(k) = k(k + 1)...(k + l − 1), we arrive at Taking the derivative with respect to k here is not easy. Because we are concerned with very low m values, we will make the approximation that k/m >> λ. This is equivalent to the approximation m kη << 1 done for the case l = 1. There we arrive at There, we have Therefore, Thus for low m, ∂φ ∂k > 0 ⇐⇒ ∂φ ∂c < 0, increased prey variability decreases the total response on the alternative prey, which is "good" from the alternative prey point of view. Confirmed with numerical simulations (Fig. 4 in the main text).
Appendix E. Seasonality in the cycle, or how to introduce delays in the aggregative response A case where the aggregative response assumption (predator numbers are a function of the current primary density) is not possible is obviously when seasonality implies delays. In the case of rodent/raptor/ground-nesting bird systems, the (cyclic) autocorrelation structure in rodent dynamics can then matter. We did not consider the fact that predators settle in early spring (as in Gilg et al., 2003, for lemming predators); when rodents decrease rather than increase over the summer, one can sometimes observe heavy predation pressure on alternative prey at the end of the summer. So far, all density values had to be thought of as end-of-spring/early-summer values. Now let us consider a two-season model with censuses for a rodent (e.g. vole) 1. At the end of winter V i,1 2. At the end of summer V i,2 The functional response is not affected by any delay, i.e. f (V i,s ) = f 0 e −λV i,s . In contrast, the aggregative response is assumed to depend on spring density only: This is simplest case of aggregative response (more complicated functions would not affect the main conclusions). In that case, we havē We have then two different cases to examine with respect to the seasonality in vole dynamics.
Appendix E.1. Rodent crashes occur mostly in winter, but declines are nonetheless possible in summer We model the dynamics of V i,1 = V i (winter-driven population cycles) and write V i,2 = R S V i,1 with R S the summer growth rate. We obtain the same model as before exceptT with a corrected functional response on the grouse that can typically be larger due to the decline in vole availability during the summer: If R S is small (e.g. 1/4), for most V values that are not too large, we typically have e −λR S V ≈ 1 >> e −λV (simply because the negative exponential function is decreasing, and relatively steep at the beginning). Thus, seasonality can have an important impact and bring the alternative prey population to levels of predation that would be unexpected given the yearly average of the vole density. However, the following case generates even stronger predation.
Appendix E.2. Rodent crashes always occur during summer: delays appear and strongly increase the level of predation Let us consider now that rodent population crashes occur during the course of the summer (as in collared lemmings, Gilg et al. 2003), which means early summer densities can reach temporarily high levels, but if end-of-winter densities are large, then end-ofsummer densities are low. Winter density is then treated as a fixed factor (or random growth rate), R W multiplied by the summer density. V i,2 = V i becomes the dynamic variable, and we have V i,1 = R W V i−1,2 = R W V i−1 . Then we obtain This model is interesting only in the case of wide amplitude fluctuations, with peaks followed by a severe crash. In this case, we can have a very high beforecrash autumn density V i−1 , generating very high abundance in the next spring, R W V i−1 . Let's assume for the sake of simplicity R W = 1 (or larger).
Following that logic, f (V i−1 ) at the end of winter/spring beginning is negligible (because f is decreasing). But then, at the end of summer the predation rate is f (V i ), which can be quite high (V i << V i−1 if the crash happens in between) and this is multiplied by the number of predators proportional to V i−1 , which is still very high.
Therefore, a delay in the aggregative response will generate much more predation when crashes occur in summer than in winter. To see how much, we have to use a model where we can manipulate amplitude easily and have control of the periodicity.
A possible candidate is the sequence V = (v 1 , v 2 , v 3 , v 1 , v 2 , v 3 , v 1 , v 2 , v 3 , ...) where v 1 = v, v 2 = Rv, v 3 = R 2 v (geometric growth with annual growth rate R combined with periodic crashes). This would lead for the average kill ratē (E.5) instead of the winter-crash case described in the previous section (E.6) We can just compare the two averages vf (Rv) + Rvf (R 2 v) + R 2 vf (v) and vf (v) + Rvf (Rv) + R 2 vf (R 2 v) • Summer crash case: the average is roughly T SC = vf (Rv)+Rvf (R 2 v)+R 2 vf (v) ≈ R 2 vf (v) (other values are negligible provided R is quite large, and most predation occurs in the year with the crash) • Winter crash case: T W C = vf (v) + Rvf (Rv) + R 2 vf (R 2 v) ≈ Rvf (Rv) (other values are negligible because v is close to zero and f (R 2 v) as well, most predation occurs in the intermediate year: remember the total response in that case has a maximum for intermediate vole values) • We know f (Rv) < f (v) because R > 1 and f is decreasing. We can even assume f (Rv) << f (v) if R >> 1. We therefore obtain T W C ≈ Rvf (Rv) << Rvf (v). Then, T W C ≈ Rvf (Rv) << R 2 vf (v) ≈ T W S (since R is in general > 10 for these species of rodents or very variable key herbivore). We can generalize the result to cycles longer than 3 years.
Conclusion: rodent dynamics with a summer crash is worse for alternative prey than rodent dynamics with a winter crash. Crash being understood here as a sharp reduction in abundance, visible in the annual time series (either spring-to-spring or autumn-to-autumn).
Considering rodent dynamics with a winter crash (the classic case when thinking about voles, but not lemmings), maximal predation will likely happen during the intermediate year: there are not enough predators at low vole density to do any damage, and the numerous predators in the peak year have very low predation on alternative prey because they prey mainly on abundant rodents. However, when instead there is a summer crash, most predation occurs right after the crash; when there are both numerous predators and no rodents, which means there is extremely heavy alternative predation.
This suggests that nest predation might be worse for ground-nesting birds when the main rodent is a lemming rather than a vole, because lemmings often crash during the summer while they increase in winter. Note: all this is conditional to predators having a strong aggregative response, although this is usually the case when it comes to lemmings in the high North (Gilg et al., 2003(Gilg et al., , 2006. Similar effects will occur when we allow for a predator numerical response, but on a multiannual rather than seasonal scale. The basic mechanism is the same: a negative temporal correlation between primary prey and predator densities. Most of the predation will then occur in the crash phase, because the predators stay quite numerous and are without anything else to eat. Unfortunately, because predation can be very strong in this case, this implies that the simplifying assumption that vole and grouse cycles are not synchronized might break down (see Discussion).

Appendix F. Synchronization of prey dynamics
Using larger densities of predators than what we observed in Scotland (Fig. E1  below), it is possible to synchronize the two prey species in our vole-harrier-grouse model, so that they fluctuate with the same period, the grouse cycle lagging however one year behind that of the voles (thus we do not really observe synchrony in a strict ecological sense, but rather "lagged synchrony"). The typical 6-year period of the grouse cycle (using parameter values of Table 1) is then reduced to the 4-year period of the vole cycle. There is change in the period of the grouse cycle generated by external forcing, as is observed for forced oscillators in physics. This phenomenon is usually referred to as frequency locking or frequency entrainment, as opposed to "phase locking" that refers to perfectly in-phase, fully synchronous dynamics (Pikovsky and Rosenblum 2007). The  Figure E1: Effect of large number of predators (H=0.5 ind.km −2 ) on the grouse periodicity (middle row), when the system is forced by a strong vole cycle (top row). The bottom row shows the grouse cycle in the absence of predation on grouse dynamical behaviour observed here is reminiscent of that of simple nonlinear forced oscillators (e.g. the standard circle map, see Weisstein).