## Abstract

We propose that delayed predator–prey models may provide superficially acceptable predictions for spurious reasons. Through experimentation and modelling, we offer a new approach: using a model experimental predator–prey system (the ciliates *Didinium* and *Paramecium*), we determine the influence of past-prey abundance at a fixed delay (approx. one generation) on both functional and numerical responses (i.e. the influence of present : past-prey abundance on ingestion and growth, respectively). We reveal a nonlinear influence of past-prey abundance on both responses, with the two responding differently. Including these responses in a model indicated that delay in the numerical response drives population oscillations, supporting the accepted (but untested) notion that reproduction, not feeding, is highly dependent on the past. We next indicate how delays impact short- and long-term population dynamics. Critically, we show that although superficially the standard (parsimonious) approach to modelling can reasonably fit independently obtained time-series data, it does so by relying on biologically unrealistic parameter values. By contrast, including our fully parametrized delayed density dependence provides a better fit, offering insights into underlying mechanisms. We therefore present a new approach to explore time-series data and a revised framework for further theoretical studies.

## 1. Introduction

Most predator–prey models assume that predators respond to present abundance of prey; for example, growth and ingestion rates are influenced by immediately available prey levels. We might, however, expect that a predator will be influenced (behaviourally, physiologically) by prey levels at some past time, possibly a day or a year, depending on its life history, storage ability, and ability to perceive and respond to changing environments. For instance, a predator may increase its fitness when prey are abundant, allowing it to be more resilient in the future when prey are scarce. Including such nutritional history or ‘delayed responses’ in models seems intuitive. Theory, in fact, indicates that allowing predators to respond to past, rather than to present, prey levels can drive systems away from equilibrium, towards cycles [1–3]. These and earlier analyses have led to a body of research on how ‘delayed density dependence’ alters population cycles, primarily by combining observations of predator–prey cycles with theory, and to a lesser extent experimentation, to fit models to data and to explore dynamics [4–12]. Critically, a common assumption in this field is that past-prey availability influences the predator's reproductive rate, rather than its ingestion rate.

Consequently, studies tend to assume the numerical response (change in predator numbers versus prey abundance) depends on past conditions, whereas the functional response (predator ingestion rate versus prey abundance) does not. Concomitantly, virtually all investigations adhere to a Lotka–Volterra-based structure that explicitly assumes the predator (*C*, consumer) population growth is directly proportional to ingestion (*I*) of prey (*P*), as in equations (1.1) and (1.2):
1.1and
1.2where *f*(*P _{t}*) is prey-specific growth rate at time

*t*;

*I*(

*P*) is the functional response;

_{t}*e*is assimilation efficiency (the fraction of ingested prey converted to predator);

*δ*is the predator mortality; and

*τ*is a time delay,

*P*

_{t −τ}being prey abundance

*τ*time steps previously. This formulation undoubtedly stems from an expectation that feeding is influenced by capture and handling of presently available prey, whereas past-prey regimes impose physiological, behavioural or perhaps maternal effects on the present reproduction. Our first concern is that this assumption has not been adequately supported though experimentation and may be invalid, as there may be an impact on prey capture and handling. Second, existing studies have not independently examined the direct influence of past-prey levels on the numerical response; rather they follow the above-mentioned logic (equations (1.1) and (1.2)), assuming the influence on reproduction is directly proportional to past feeding, which may not be so, as

*e*and

*δ*will vary with prey abundance [13,14]. Finally, many studies allow the time delay (

*τ*) to vary, to fit time-series data, following an information criterion approach [15].

Although the approach described by equations (1.1) and (1.2) can be insightful, it will potentially provide best-fit parameter values that are not biologically meaningful; for example, estimates of *τ*, *e* and *δ* may be unrealistic. If this occurs, then models may be generated that are phenomenologically adequate to fit time-series but are mechanistically uninformative, and therefore inadequate for providing general insights regarding specific processes or applying the model in a more general manner.

To expand on and complement past studies, a more mechanistic approach is required: one that directly quantifies the influence of past-prey levels separately on the functional and numerical responses, and compares their predictions with independently derived time-series data (i.e. in-sample fitting with out-of-sample validation). This view was stimulated by fieldwork and associated models on lemmings and their predators [16,17]. We therefore offer a new strategy: independent assessment of the functional and numerical responses coupled with the impact of an experimentally dictated time-delay on a range of past-prey availabilities. Fenton *et al.* [13] explicitly formalized a model structure that offers a framework amenable to this approach. Their *independent response* model incorporates independently derived functional and numerical responses that intrinsically allow *e* and *δ* to vary with prey abundance, providing a mechanistically more realistic model [14,15,18,19]. This offers an empirically motivated counterpoint to the generally applied approach to biomass conversion that arises within the standard models, where growth is a linear function of consumption [20]. Previous applications of the independent response model, however, have not considered past-prey abundance. Here, we take this framework and experimentally impose the effect of past-prey abundance on the functional and numerical responses. Then, by comparing model predictions, based on the parametrized responses, with independently obtained time-series observations, we provide insights into the process driving predator–prey cycles.

Specifically, we determine the influence of past-prey abundance on the functional and numerical responses by imposing a fixed delay in the feeding and reproductive rates. This differs from previous studies that fit responses to existing cycles; rather than assessing the influence of delay duration (parametrizing *τ*, equation (1.2)), we experimentally fix the delay at a value of the order of one generation (a period expected to influence the present generation; for support of this approach, see [8]) and assess the magnitude of the numerical and functional responses to that delay. By including these comprehensive responses in a population model, we then provide direct, controlled evidence that functional and numerical responses that depend on delayed prey densities promote population cycles. Thus, we first rigorously support the poorly evaluated notion that reproduction, not feeding, is highly dependent on the past food regime. We then indicate that this delay in reproduction affects both short-term (transient) and long-term (equilibrium) population dynamics. Critically, we show that although superficially the standard approach (equations (1.1) and (1.2)) may provide a parsimonious fit to time-series data, it does so by applying unrealistic parameter estimates (i.e. a reasonable conclusion is obtained from poor assumptions). By contrast, our mechanistic approach provides a better fit to the time-series data and insights into underlying mechanisms. This new approach offers a means for future exploration of existing oscillating datasets and a revised framework for consideration in theoretical studies.

## 2. Material and methods

### (a) Experimental organisms and medium

We develop and parametrize models using an experimental system that has elucidated population dynamics for almost a century (see electronic supplementary material S1): the ciliates *Didinium nasutum* (predator) and *Paramecium caudatum* (prey). Prey were maintained in Volvic water containing 0.55 g l^{−1} protozoan pellets (Carolina Biological Supply, USA); crushed pellets were boiled for 30 min; cooled medium was inoculated with the bacterium *Serratia marcescens*. Predators were maintained in this medium and fed prey.

Experiments used this medium thickened with methyl cellulose (M_{n} 40 kDa, Sigma-Aldrich, St Louis, MO) with viscosity of 3.5 cP; this slows organisms (allowing accurate counts), increases predator–prey cycles persistence and parallels other work on this system [21–23]. High bacteria levels were maintained, ensuring maximum prey growth (data not shown). Experiments were conducted at 20°C in 6 ml of medium. To ensure random distribution (verified visually), containers were constantly but gently agitated.

### (b) Laboratory analyses

#### (i) Time series and prey logistic growth

Time series (*n* = 5) were initiated at one predator and 10 prey ml^{−1}, and followed for 12 days, enumerating microscopically thrice daily. To obtain prey logistic growth parameter estimates, replicate cultures (*n* = 3) were enumerated twice daily until abundance remained at the carrying capacity. For both of these, to reduce toxin accumulation, every second day 50% of the medium was replaced [23]. The solution to equation (2.1) was fitted to prey abundance (*P,* ml^{−1}) data by the Marquardt–Levenberg least-squares algorithm (SigmaPlot v. 11, SPSS Inc., Chicago, IL),
2.1where *μ* (d^{−1}) is prey-specific growth rate and *K* (ml^{−1}) is the carrying capacity; mean values (±s.e.) were determined from replicates.

#### (ii) Effect of past-prey abundance on predator ingestion and growth rates

The influence of past-prey abundance (*A* = past : present prey, from 0.25 to 2) on predator-specific growth (*r*, d^{−1}) and ingestion rate (*I*) was determined as follows: predators experienced growth-saturating levels (more than 80 prey ml^{−1}) for 2 days, to ensure they were similar prior to manipulations. Then, predators experienced a fraction of the present-day-prey level for 2 days. Finally, predators experienced a present-day-prey level, at which specific growth and ingestion were determined (e.g. a final present-day level of 100 prey ml^{−1} previously experiencing 200 prey ml^{−1} for 2 days provided conditions for measuring growth or grazing at 100 prey ml^{−1}; for these measurements *A* = 2).

#### (iii) Baseline functional response

Prey were labelled by feeding them 3 μm fluorescent microspheres (Sigma-Aldrich) for 90 min at approximately 10^{5} microspheres ml^{−1}, allowing visualization of prey within predators [24]. Forty predators (conditioned to a past-prey abundance) were inoculated into medium containing a defined number of labelled prey. After incubation (time-dependent on prey abundance but ensuring less than 30% of predators fed), treatments were glutaraldehyde-fixed (2%) and examined (epifluorescent microscopy) to detect ingested prey. Ingestion rate was determined as the fraction of predators that had consumed prey over the incubation time, assuming a predator ingested no more than one prey during the incubation (30 min to 2 h), which is supported by our results. Prey ranged from 1 to 150 ml^{−1} (prey carrying capacity, approx. 110 ml^{−1}) with many concentrations (see §3) and no replicates [25]; average prey and predator abundance, used to determine responses over the incubation, were obtained from initial and final estimates, following standard methods [26].

The response of predator ingestion rate (*I, P* d^{−1}) to prey abundance (*P,* ml^{−1}) at each level of *A* was determined by fitting a type II functional response (equation (2.2)) to the data,
2.2where *I*_{max} = maximum ingestion rate (*P,* d^{−1}), and *k* (ml^{−1}) is the half-saturation constant. As the data reflected a type II response (see §3), no others were considered. Responses were fit using the Marquardt–Levenberg least-squares algorithm, appropriate for fitting models to such data [27].

#### (iv) Baseline numerical response

Reproductive rate at different prey concentrations was determined for predators that had experienced past-prey levels (as above), following standard methods [28], with average prey abundance determined [26]. Predators experienced prey ranging from 1 to 150 ml^{−1}. Predator-specific growth rate (*r*, d^{−1}) was determined over 24 h, following: *r* = ln(*N _{t}*/

*N*

_{o})/

*t*, where

*N*and

_{t}*N*

_{o}are, respectively, final and initial abundance, and

*t*is time (days).

The response of specific growth rate to prey abundance (*P,* ml^{−1}) at each level of *A* was determined by fitting a modified type II response (equation (2.3)) to the data,
2.3where *r*_{max} is maximum-specific growth rate (d^{−1}), *k*_{1} is a constant (*P,* ml^{−1}) and *P′* is threshold prey level (where *r* = 0). This semi-mechanistic function captures predator reproduction and mortality rates, and implicitly incorporates variable assimilation efficiencies and variable mortality rate with prey abundance [13,14].

### (c) Data analyses

#### (i) Incorporating past-prey abundance into functional and numerical responses

For modelling, functions were established to predict predator ingestion and growth rates in response to varying present : past-prey abundance (*A*). Functions were fitted to the *A*-influenced response data. The asymptotic maximum rates were assumed to not be directly influenced by *A*; rather, the half-saturation constant (*k*) for the functional response and the constant *k*_{1} and the threshold concentration (*P′*) of the numerical response were considered to be influenced by *A* (see equations (2.4) and (2.5)); these alter the curvature of the responses, and the rate at which the theoretical maximum is approached. Several variants of the general structure of the responses were examined that included nonlinear elements (power, exponential, polynomial functions), and the most parsimonious yielding the best fit (judged by the adjusted *R*^{2}) was applied. This produced semi-mechanistic responses that reveal trends and reflect biological phenomena (i.e. nonlinear responses of growth and ingestion rates to changing levels of *A*; see §3). The *A*-influenced functions for predator ingestion (*I _{A}* (equation (2.4)) and growth (

*r*, equation (2.5)) include the constants

_{A}*α*,

*β*,

*ϕ*,

*λ, a, b, c, d*and

*f*. Adjusted

*R*

^{2}-values for the responses and standard errors of the estimates were determined (SigmaPlot), to indicate their goodness of fit. 2.4and 2.5

#### (ii) Combining data and theory: comparison of the two approaches

To reiterate, we are exploring two fundamentally different approaches to incorporate delays in predator–prey models: (i) the standard (parsimonious, biomass conservation [20]) approach (equations (1.1) and (1.2)), with the numerical response directly proportional to the functional response at some previous time (*τ*); and (ii) the independent response (less parsimonious but more mechanistic [13]) approach, with independently derived functional and numerical responses that incorporate empirically derived functions of past prey at a defined time via the acclimation function *A* (past : present-prey ratio).

Recognizing the difference between the model parametrization and structure, we compare these approaches in terms of fits with time-series data and evaluate the long-term population dynamic behaviour as predicted by each. We adopted a trajectory-fitting approach, using the R package FME [29], to fit the predator and prey time-series predicted by the different models to the data (for details, see electronic supplementary material S2). In brief, independent fits were made for each replicate time series (*n* = 5); best fits of each parameter and initial densities were calculated across these, and the two approaches were compared based on *R*^{2}-values between the best-fit models and data (paired *t*-test, *α* = 0.05). The best-fit versions of the models were then simulated over 200 days to explore long-term dynamical behaviour (i.e. whether stable equilibria, damped oscillations or sustained oscillations were predicted).

#### (iii) Exploring the standard delay equations

To assess the efficacy of the standard approach (equations (1.1) and (1.2)), we examined its ability to mimic the time-series data, where and *I*(*P _{t}*) is as in equation (2.2). Parameter values for equation (2.2) were obtained from the fully parametrized equation (2.4), by setting

*A*= 1. This model was fitted to the data under two scenarios: (i) no time delay (

*τ*= 0) and (ii)

*τ*estimated through fitting. In both cases, assimilation efficiency (

*e*), mortality rate (

*δ*), and initial prey and predator abundance were allowed to vary, to provide best estimates; the other parameter values in the model (e.g. those in equation (2.4) and prey growth equation, above) were held constant at their estimated values.

#### (iv) Exploration of *A*-influenced functional and numerical response

To examine the influence of *A* on population dynamics and assess if its inclusion improved the model's predictive ability, we examined four scenarios following the independent response structure [13]. In all of these, prey grew following parameters derived from equation (2.1) and were consumed based on parameters derived from equation (2.4); predator growth rate was prey-dependent with negative growth below threshold levels, predicted from parameters derived from equation (2.5). Equations (2.6) and (2.7) describe the model, with parameters and variables as above:
2.6and
2.7

Four scenarios were examined to assess whether *A*-influenced ingestion, growth or both need be considered in models: *A* = 1 for both functional and numerical responses (equations (2.4) and (2.5)); *A* = 1 for the numerical response (equation (2.5)), with the fully parametrized functional response following equation (2.4); *A* = 1 for the functional response (equation (2.4)), with the fully parametrized numerical response following equation (2.5); and *A* varied according to fully parametrized functional and numerical responses (equations (2.4) and (2.5)). As including an *A*-influenced numerical response was required to best mimic the time-series data (see §3), we explored the best fit of the fully parametrized model to the time-series data by allowing parameters that altered the curvature of the numerical response (*b*, *c*, *d*, *f*; equation (2.5)) and the initial prey and predator abundance to vary in the fitting process.

## 3. Results

### (a) Time-series and prey logistic growth

Time-series incubations followed a replicable pattern (figure 1), prey populations becoming extinct within one cycle (fewer than 8 days), with predator extinction following. For prey logistic growth (equation (2.1)), growth rate (*μ*) = 0.70 ± 0.04 d^{−1} and carrying capacity (*K*) = 109 ± 0.7 ml^{−1}.

### (b) Functional and numerical responses

Although ingestion and specific growth rate both followed rectangular hyperbolic responses to prey abundance, the predicted ingestion rate responses did not appear asymptotic over the range examined (figure 2); the two responses were distinctly different in shape (figures 2–4).

Both responses were influenced by the past : present-prey ratio (*A*), but in distinctly different ways. For the functional response, the half-saturation constant (*k*) did not vary with *A*, except at low levels, where it was unrealistically high, as no asymptote was reached (figures 2 and 3*a*). The maximum observed ingestion rate (predicted at the highest prey level) initially rapidly increased with *A* and then gradually increased as *A* increased from 0.5 to 2.0 (figure 3*b*). There was a stronger effect of *A* on the shape of the numerical response (figures 2 and 3): *k*_{1} monotonically decreased with *A*; *P′* declined in an exponential manner with increasing *A*; and the maximum observed growth rate increased with *A* (figure 3*c–e*). Predictive functions for the interactive effects of prey abundance (*P,* ml^{−1}) and *A* on ingestion and growth rate (equations (2.4) and (2.5)) provided good fits to the data (see figure 4; electronic supplementary material S3).

### (c) Model analysis

Model scenarios (figure 5*c–k*) were compared with time-series data (figure 5*a,b*). Model simulations were in an infinite volume, while empirical data were in 6 ml; this allowed simulations to reach lower, but reasonable, numbers (approx. 1 l^{−1}), and thus recover, potentially producing cycles.

#### (i) Exploration of the standard delay equations

The best fit of equations (1.1) and (1.2) to the data with *τ* = 0 estimated assimilation efficiency (*e*) = 0.25 ± 0.01 and mortality rate (*δ*, d^{−1}) = 0.68 ± 0.14. Long-term dynamics of this model resulted in stable equilibrium at approximately 5 predators ml^{−1} and approximately 25 prey ml^{−1} (figure 5*c*). This prediction is unprecedented, as no incubation that we have conducted over several years has yielded stable equilibrium at any levels (figure 1; D.J.S.M. 2010–2013, unpublished data).

When *τ* was allowed to vary, the best fit of equations (1.1) and (1.2) estimated *e* = 0.22 ± 0.01, *δ* = 0.28 ± 0.01 d^{−1} and *τ* = 0.94 ± 0.03 days; these predicted long-term sustained predator–prey cycles (figure 5*d*). In the short term, the predicted dynamics differed from the empirical data in terms of the magnitude of predator and prey abundance and the timing of cycles (figure 5*b,e*). The mean *R*^{2}-value for the model fit to the data was 0.609 ± 0.027.

#### (ii) Exploration of past- and present-prey abundance

Simulations with no influence of past-prey abundance imposed on both the functional and numerical response (*A* = 1), or with *A* = 1 for the numerical response (equation (2.5)) but varying in the functional response (equation (2.4)), exhibited virtually identical behaviour (figure 5*f,g*). Two features not exhibited by the empirical data appeared: prey cycle amplitude was approximately 50% lower, and cycles converged at approximately 5 predators ml^{−1} and approximately 10 prey ml^{−1} (figure 5*a,f,g*); again this has never been empirically observed.

Simulations with past-prey abundance imposed only on the numerical response (*A* = 1 for equation (2.4) and varying according to equation (2.5)), or on both responses (*A* varying according to equations (2.4) and (2.5)), predicted similar, large oscillations with sustained predator–prey cycles (figure 5*h,i*). Fitting the latter model to the data and allowing the parameters of the numerical response, *b*, *c*, *d* and *f* (equation (2.5)), to vary (best-fit values: *b* = 11.3 ± 0.46, *c* = −0.55 ± 0.04, *d* = 33.8 ± 3.7, *f* = 0.44 ± 0.10), the model achieved the closest fit to the empirical data (figure 5*a,j,k*); it is important to note that these values are similar to those obtained through experimentation (see electronic supplementary material S3), suggesting our empirically derived estimates are appropriate. Although the model underestimated peak prey abundance, it provided a better estimate of predator abundance than the standard delay model, and the timing of population cycles was slightly improved over that of the standard delay model (cf. figure 5*b,e,k*). The mean *R*^{2}-value for this model fit was 0.749 ± 0.026, significantly higher than that obtained for the standard structure that included *τ*; *t* = 5.595, *p* = 0.0005, with the caveat that it contains more parameters. This suggests that the independent response model provides a better (albeit less parsimonious; see §4) predictor than the standard model.

## 4. Discussion

Time-delay effects on either prey or predators can influence population dynamics, with the potential to destabilize systems causing cycles [1,10,16,17,30,31]. Often, however, delays are incorporated phenomenologically by introducing a lag (*τ*, equation (1.2)) that governs predator-specific growth rate according to past-prey levels. This flexible approach allows a range of dynamic behaviours by varying *τ*, often enabling population models to fit time-series data. However, it is not necessarily clear which biological process is being described by the lag or whether the fitted values are biologically meaningful. For example, using our model system (*Paramecium*, *Didinium*; electronic supplementary material S1), Jost & Ellner [8] significantly extended the general concepts addressed by Harrison's work on the same system [4] and concluded, among other aspects, that imposing lags on *Didinium* growth rate significantly improved model fits to time-series data from a MSc thesis by B. G. Veilleux (see [7,23] and http://robjhyndman.com/tsdldata/data/veilleux.dat). These studies, however, lacked data on the influence of the past-prey regime; rather they explored a range of scenarios to explain a best fit (e.g. type II and III functional responses).

In the decade since Jost & Ellner's work [8], there have been numerous theoretical explorations inspired by and directly evaluating *Paramecium–Didinium* time-series data [32–38], but these and similar works tend to lack mechanistic underpinning. They therefore may arrive at what seem as sensible conclusions for the ‘wrong’ reasons (i.e. they may obtain biologically unrealistic parameter values or include unrealistically parametrized functions). Applying such incorrect values or functions would be counterproductive to understanding the system, making it unlikely that the resulting model would correctly predict the outcome of perturbations to the system. By contrast, a more mechanistic, complex model may allow this (see [39] for similar arguments). Here, we adopted a novel approach, using experimentally derived relationships that directly capture the nonlinear impact of past-prey abundance on the functional and numerical responses and assess the ability of a model, based on these, to fit predator–prey population dynamics.

We observed that both the functional and numerical responses were affected by past-prey levels, and therefore suggest that both should be independently considered in delay-dependent models. This differs from many works that ignore past-prey abundance on feeding. However, the numerical response had a greater dependence on past-prey levels than the functional response (figures 2 and 3), supporting the notion that growth is strongly dependent on the past, whereas feeding depends on the present-prey regime. It also superficially supports the long-standing approach to modelling that does not independently examine the numerical response (equations (1.1) and (1.2)).

However, we then provide strong evidence that this standard structure (equations (1.1) and (1.2)) requires reconsideration. First, past-prey abundance influences the functional and numerical responses in distinctly different ways (figures 3 and 4). Consequently, it may be inappropriate to apply the logic of equation (1.2), where the numerical response is directly derived from the functional response; instead, the numerical response should be independently assessed [13,14,18,19].

Clearly, including no time delay in either the functional or numerical response produces erroneous behaviour, predicting that the predator–prey system should reach a stable equilibrium (figure 5*a–c*). Note that, because our time-series data only lasted for one cycle, it is impossible to categorically state whether these data represent transient, damped oscillations or sustained limit cycles; however, stable equilibrial dynamics have never previously been observed in longer incubations of this system, and we have no reason to expect that would be the case here. By including a delay into the standard framework (equation (1.2)), predator–prey dynamics that more closely mimicked the empirical data were obtained (figure 5*a,b,d*), but predator and prey abundance (and to some extent the timing of the cycles) differed from the empirical data; furthermore, the parameter values that yield these cycles were unrealistic. Both assimilation efficiency (*e*) and mortality rate (*δ*) vary with prey abundance [13,14,40]; specifically, *δ*, for which we have *Didinium-*based data, may range from less than 0.1 to 1.6 d^{−1} [14]. Consequently, the values that provided the best fit (*e* = 0.22 and *δ* = 0.27 d^{−1}) are likely to be biologically unrealistic. This illustrates our main concern: that due to the flexibility of the standard, parsimonious model, it is possible to obtain the right answer (a good fit to the data) for the wrong reasons (inappropriate model parametrization).

We next explored the influence of past-prey abundance on population dynamics using the independent response model that accounts for variables *e* and *δ* [13]. This indicated that past-prey abundance should be considered, and it is the influence of past-prey abundance on the numerical response that leads the model to best represent the empirical data (figure 5); this better fit (than the standard framework) undoubtedly arises from an increase in parameters. However, it has the advantage that the parameter values, and the more complex model structure, are biologically meaningful, reflecting large, nonlinear trends in the experimentally obtained data (figures 2⇑–4). This illustrates the need to respect parsimony in model construction only where it is appropriate and suggests that more complex and mechanistically realistic models must be applied when they are justified by empirical work, as we have done here. By developing and experimentally parametrizing these models, it is now possible to obtain insights into, and biological understanding of, the process driving the observed predator–prey dynamics, thus getting the right answer for the right reasons. Given that, we might then explore those reasons: delay-driven responses could include reserves of well-fed individuals [4], changes in swimming behaviour, predator : prey ratio dependence through a range of processes [8,41], the role of demographic or environmental stochasticity in driving those dynamics [42], or changes in the life cycle, all of which would require careful experimentation to assess.

While it will be impossible to obtain such data for all predator species, the framework developed and principles revealed here should be applicable to a wider range of predator–prey systems. Finally, we suggest that to evaluate predator–prey dynamics, exploring the interaction between past- and present-prey abundance on both growth and ingestion at carefully chosen levels of *τ* may be as fruitful as other presently applied approaches.

## Funding statement

This study was supported by a grant from a British Ecological Society, Small Project and the National Nature Science Foundation of China (no. 31071893).

## Acknowledgements

We thank reviewers and colleagues who provided comments.

- Received June 3, 2013.
- Accepted July 15, 2013.

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