Fitting stochastic predator-prey models using both population density and kill rate data

04/03/2019
by   Frederic Barraquand, et al.
Université de Bordeaux
0

Most mechanistic predator-prey modelling has involved either parameterization from process rate data or inverse modelling. Here, we take a median road: we aim at identifying the potential benefits of combining datasets, when both population growth and predation processes are viewed as stochastic. We fit a discrete-time, stochastic predator-prey model of the Leslie type to simulated time series of densities and kill rate data. Our model has both environmental stochasticity in the growth rates and interaction stochasticity, i.e., a stochastic functional response. We examine what the kill rate data brings to the quality of the estimates, and whether estimation is possible (for various time series lengths) solely with time series of counts or biomass data. Both Bayesian and frequentist estimation are performed. The Fisher Information Matrix suggests that models with and without kill rate data are all identifiable, although correlations remain between parameters that belong to the same functional form. However, our results show that if the attractor is a fixed point in the absence of stochasticity, identifying parameters in practice requires in fact kill rate data as a complement to the time series of population densities, due to the relatively flat likelihood in this case. Only noisy limit cycle attractors can be identified directly from ecological count data (as in inverse modelling), although even in this case, adding kill rate data - including in small amounts - can make the estimates much more precise. Our framework highlights how to combine data types in dynamical models for multiple species, and may be extended to other biotic interactions for which data on both interaction rates and population counts are available.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 15

page 29

01/05/2022

The effect of biologically mediated decay rates on modelling soil carbon sequestration in agricultural settings

Microbial biomass carbon (MBC), a crucial soil labile carbon fraction, i...
04/28/2021

Estimation of Poisson Autoregressive Model for Multiple Time Series

A Poisson autoregressive (PAR) model accounting for discreteness and aut...
06/13/2018

PoARX Modelling for Multivariate Count Time Series

This paper introduces multivariate Poisson autoregressive models with ex...
08/08/2019

Application of Levy Processes in Modelling (Geodetic) Time Series With Mixed Spectra

Recently, various models have been developed, including the fractional B...
10/31/2018

A Sequential Design Approach for Calibrating a Dynamic Population Growth Model

A comprehensive understanding of the population growth of a variety of p...
03/09/2020

Strong mixing condition for Hawkes processes and application to Whittle estimation from count data

This paper focuses on the time series generated by the event counts of s...
11/01/2020

A general modelling framework for open wildlife populations based on the Polya Tree prior

Wildlife monitoring for open populations can be performed using a number...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

Introduction

The parameterization of dynamical systems for interacting species usually involves either independent data on process rates (e.g., reproduction and kill rates, Turchin and Hanski, 1997) or some form of inverse modelling (Turchin and Ellner, 2000; Jost and Ellner, 2000; Jost and Arditi, 2001 or more recently, DeLong et al., 2018; Rosenbaum et al., 2019; Curtsdotter et al., 2019). In inverse modelling approaches, a dynamical system (often deterministic) is adjusted to time series of population counts, densities or biomasses (Stouffer, 2019). However, even in a perfect deterministic world, inverse modelling from time series is challenging due to numerous identifiability issues (Raue et al., 2009; Eisenberg and Hayashi, 2014; Kao and Eisenberg, 2018), that is, several sets of parameters can produce identical time series. For some models, even with an unlimited amount of data, two parameters cannot be separated from each other because they contain redundant information (Cole et al., 2010; Little et al., 2010), a problem that has long been known to ecologists working on capture-recapture models (Catchpole and Morgan, 1997).

An additional concern when modelling ecological communities over time is that environmental stochasticity usually has a pervading influence on their dynamics (Lande et al., 2003; Mutshinda et al., 2009; Mac Nally et al., 2010). It makes therefore great sense to fit community-level models that allow for such environmental stochasticity () to perturb the dynamics, now creating a stochastic mapping between state variables at to those at . However, this increases rather than decreases the complexity of the inverse problem, since environmental stochasticity can then combine in non-intuitive ways in nonlinear models (Greenman and Benton, 2005). While in an observational stochasticity setting (e.g., Rosenbaum et al., 2019) effects of stochasticity do not accumulate, with environmental stochasticity, perturbations to the growth rates of species within the community accumulate over time (Jost and Arditi, 2001; Turchin and Ellner, 2000), creating even more parameter configurations possibly leading to similar time series. Even in simple, phenomenological statistical models for two species in interaction, there can therefore remain considerable uncertainty about the model structure and parameters, especially for short ecological time series of length time steps (Barraquand and Nielsen, 2018).

How to decrease the uncertainty on parameter estimates in stochastic community models? One strategy that has enjoyed a great success in population ecology is to use a combination of data sets (Besbeas et al., 2002; Maunder, 2004; Schaub and Abadi, 2011). This allows to increase the precision of parameters that could be estimated (but with a large uncertainty) and often to estimate parameters that were not identifiable when the models were fitted to one data type only (Besbeas et al., 2002). In the case of predator-prey modelling, such strategy would be equivalent to combining process rate data with population densities, which is what we attempt here. Ours is not, of course, the first study to consider the relative benefits of forward modelling (process rate modelling) vs inverse modelling (e.g., Turchin and Ellner, 2000). The novel perspective here is on merging different sources of data on process rates (i.e., the kill rate or biomass flow) and densities (or biomasses) into a fully stochastic framework.

Indeed, even in models that include environmental stochasticity on the population growth rates of predators and prey (e.g., Karban and De Valpine, 2010; Ives et al., 2008), the functional response (Solomon, 1949; Holling, 1959) describing the relationship between predator kill rate and population densities is almost always viewed as a deterministic object (outside of corrections for sampling variation, e.g., Fenlon and Faddy, 2006). Exceptions to this rule are to be found in continuous-time ODE models where stochasticity has been added to all parameters (Turchin and Hanski, 1997; Turchin and Ellner, 2000) and in stochastic individual-based predator-prey models (e.g., de Roos et al., 1991; Wilson, 1998; Hosseini, 2003; Murrell, 2005), in which interaction stochasticity naturally arises from the encounter process. By contrast, in most models using stochastic difference or differential equations, kill rates of individual predators are the deterministic outcome of prey population densities. However, in addition to the vagaries of the interaction process (interaction stochasticity per se, due to the randomness of encounters of pairs of individuals), the predator kill rate typically depends on multiple quantities that vary over time (and space, though we will not discuss this here): the densities of several species and not only the main prey (Abrams, 1982, 2010), abiotic factors such as temperature (Rall et al., 2012), traits such as body size (Vonesh and Bolker, 2005), and the list goes on. With such an imperfect knowledge of all factors affecting the kill rate of individual predators, and measures of population densities restricted to two species (the predator and its prey), it can therefore make more sense to view the functional response as a stochastic object. There are also empirical arguments for this view: point clouds surrounding the functional response curve usually display considerable variation (e.g., Fig. 3 of Kalinkat et al., 2013, Fig. 3 of Pritchard et al., 2017, Fig. 5 in Aldebert and Stouffer, 2018, Fig. 5 of Rosenbaum and Rall, 2018). Considering the functional response (FR) as a stochastic rather than deterministic object is therefore both biologically realistic and statistically convenient. In fact, some researchers have chosen to move away from the classical food web modelling paradigm precisely because of its reliance on deterministic functional relationships (Planque et al., 2014). To our knowledge, embracing fully the stochasticity of the functional response in a classical parametric inference context has only been done by Gilioli et al. (2008, 2012). However, they did not use kill rate data to inform their stochastic functional response as we do here. Note well that by stochasticity of the functional response, we mean that the kill rate can change from one time step to the next in the absence of changes in population densities, and not simply that there is a large uncertainty over parameter values (on the latter aspect, see e.g., Smout et al., 2010; Aldebert and Stouffer, 2018).

In light of the above considerations, to integrate community-level time series of densities and kill rate data, we consider a ground truth dynamical model for predator-prey systems with stochasticity in both growth rates and interactions parameters, in discrete time. We first simulate datasets using the model, and then fit this model and a variant to the simulated data. We use both frequentist and Bayesian estimation paradigms, to check that our results are congruent between the two frameworks and therefore fairly robust. Identifiability is examined in both frameworks, using the Fisher Information Matrix in a frequentist setting and the posterior-prior overlap in a Bayesian setting, with or without kill rate data. We consider both long time series (T = 100) by ecological standards and then time series of realistic ecological length (T = 50 or 25). Finally, we vary the percentage of data points for which kill rate data are available, which can contribute to a more optimal allocation of time and effort in the field, when designing surveys of community dynamics.

Models and statistical methods

Predator-prey model in discrete time

We chose a model with a numerical response of the Leslie type, but similar analyses are performed for Rosenzweig MacArthur models in Supplement S1, to ensure generality. A Beverton-Holt function for density-dependence in the prey growth rate was chosen to avoid cycles in the prey in absence of the predator, so that the model behaviour is more reminiscent of its model continuous-time counterpart (see Weide et al., 2018 on connecting discrete-time to continuous-time predator-prey models). Our model writes

(1)
(2)

The roots of this discrete-time formulation can be traced back to Leslie (1948); Leslie and Gower (1960) who included a Beverton-Holt regulation for the prey and predator to make discrete-time models more similar to their continous-time counterparts. We consider additionally a saturating functional response here, which makes our model resembles the continuous time models of Tanner (1975) or May (1973) and later Turchin and Hanski (1997), whose notations we have kept. Except that now, the functional response is actually , with a stochastic term included.

One of the advantages of the Leslie-type parameterization over the Rosenzweig-MacArthur linear numerical response is that the former can be parameterized using the observed predator maximal intrinsic growth rates as Tanner (1975) did. Parameter values were inspired by small mammals (Turchin and Hanski, 1997), which we modified slightly to be able to get both stable and cyclic dynamics. The division by in expresses the fact that all terms within the exponential are on the prey fitness scale (per capita mortality), a very classic representation (e.g., Ives, 1995). This is similar to assuming a Nicholson-Bailey type predation term (Weide et al., 2018).

Until now, we have not specified a model for the functional response . With a deterministic functional response (FR), we have a classic stochastic predator-prey model with log-normal environmental noise but an otherwise ‘deterministic skeleton’. A stability analysis of the deterministic skeleton of the model is performed in Appendix 1, and used to determine which parameters lead to a perturbed fixed focus or node, or a noisy limit cycles in the stochastic model. Here we consider a data-generating process where the functional response is not deterministic but itself stochastic, as in the following equation for a stochastic Holling type II functional response:

(3)

This corresponds to a case where there is mild Gaussian fluctuations around the functional response. Because there can be substantial noise on the FR, we also consider more complicated cases where the parameters and are themselves allowed to vary in time in Supplement S2.

Statistical methodology

Although we apply here data integration to a predator-prey case, the methodology is more general and can in principle be applied to any addition of auxiliary information (interaction rate, demographic rate) which is available over time and not simply produced by the count data. The approach has similarities to Integrated Population Modelling (Besbeas et al., 2002; Maunder, 2004; Péron and Koons, 2012) and data fusion/assimilation approaches (Niu et al., 2014). In a predator-prey context, we can cite the seal predation model of Cook et al. (2015) or the recent endeavour of Ferguson et al. (2018) to combine count and isotopic (diet) data, but their models differ in that they do not view interaction rate as a fully stochastic process.

We illustrate the method with the predator-prey case where log-densities for both predator and prey are gathered in a log-count vector

. Auxiliary information on the functional response, or rather the observed kill rate per predator is denoted . To this can be added other demographic (vital) rates that are stacked in a vector as well, . Currently we use but it may useful to add more information in other applications; hence the derivation is kept general.

We consider a discrete-time dynamical system (nonlinear difference equation). The population dynamics part of the model gives us the probability law of

, since the counts at time depend both on past abundances and the interaction or demographic data based on the chosen mechanistic model. We also know the probability law of (in our simple predator-prey case, the functional response). We can therefore write down easily the joint likelihood for both data sources in quite general terms, denoting and :

(4)

where and are continuous probability densities for the matrix

and its conditional pdf, respectively. The conditional pdf can be further decomposed using the chain rule of probability (i.e., relying on conditional independence)

(5)

where represents the dynamical system, given auxiliary information on process rates (here, the functional response), and is the functional response model (or a demographic model), which is conditional only to the current set of abundances. We therefore end up with a closed form expression for the likelihood

(6)

where we swapped and to get the dynamical system model first.

Application to the stochastic predator-prey model

In the simplest case highlighted by our discrete-time dynamical systems of the sections above, is the product of the two gaussian pdf for log-densities conditional on past densities and auxiliary information. Using the equations (1) and (2), we obtain the following difference equation on the log-scale:

(7)

and

(8)

with a functional response model ( in eq. 6) also given by a Gaussian pdf (this is the simplest case where we assume near Gaussian noise on the functional response, see Discussion)

We considered two parameter sets for the model, given in Table 1.

Parameter Perturbed FP Noisy LC
C 2.5 15
D 1.0 0.25
Table 1: Parameter values for both types of attractors considered (Perturbed FP: Perturbed fixed point; noisy LC: limit cycle perturbed by noise). The rest of the parameters are .

The first parameter set (FP) corresponds to a forced fixed point (focus), and although the eigenvalues are imaginary (Appendix 

1) in this case there are no discernable cyclic oscillations (i.e., no quasi-cycles sensu Nisbet and Gurney, 1982). This is a crucial example because not all fluctuating predator-prey systems give rise to limit cycle oscillations. We also consider a noisy limit cycle (LC), i.e., parameters that give rise to a limit cycle without the noise, so that the cycle has a broad amplitude and well-defined attractor even as noise is added, but an irregular periodicity due to the random perturbations. We use here the wording ‘limit cycle’ by analogy to the continuous time theory, although quite formally, this parameter set gives rise to an invariant loop rather than a discrete-time limit cycle (see, e.g., Caswell, 2001, chapter 16). Representative time series of the model for both parameter sets are provided in Fig. 1.

Figure 1: System dynamics for T=100 (1 simulation). In panel (A), we show densities of prey in blue and predator in red , for the perturbed FP case, (B) corresponding trajectories in phase plane and (C) functional response: kill rate of individual predator as a function of prey density. In the second row, identical panels (D-E-F) for the noisy LC simulation.

These two sets of parameters are crossed with several modalities of data availability:

  • time series of length , and

  • we consider that kill rate (KR) data is available along a fraction , , or of the time series data. This is meant to emulate common scenarios in which the kill rates are not monitored over the whole time series, and to quantify the benefits of adding small amounts kill rate data.

Formal identifiability results involving the FIM or prior-posterior overlap are reported for and or , and will be presented first. Throughout the manuscript, KR data will refer only to the predator intake rates whilst functional response (FR) data will refer to the pair .

Note that in the case without KR data (), we fit the model without noise on the functional response, since there is no data to fuel the statistical model of the functional response. This is what has usually been done in this case in ecology, thus we wanted to have a meaningful comparison to the rest of the literature. Having a latent, completely unobserved stochastic state variable for the kill rate was considered in early simulations but typically leads to even worse estimates, and is more complicated to fit (the case does consider an unobserved latent state of 75% of the time series though).

For each parameter data availability scenario, we fit the models in a Bayesian framework. A Bayesian framework is here more convenient as requires to estimate a latent state for unobserved kill rate data (it is doable in a frequentist setting using an EM algorithm, with more effort). Because we consider 100 simulations for each parameter x data availability scenarios, we have distributions even for the point estimates (means over the posteriors), whose upper and lower bounds reveal the precision of the estimator in a frequentist sense.

Model fitting implementation

We fitted the model by MCMC in JAGS (Plummer et al., 2003), using 3 chains and vague priors (for instance for parameters and ), full code provided in a GitHub repository111https://github.com/fbarraquand/predatorpreyDynamics_wFRdata/figures_article, made public upon acceptance will all details of the implementation. We also derived mathematically the full model likelihood, which we then implemented in a function to plot likelihood surfaces and ascertain identifiability from a frequentist viewpoint. Having a well-defined likelihood function also allows to find point estimates of the parameters through quasi-Newton methods (hill-climbing algorithms), which we have done using the BFGS algorithm with optim() in R. These matched qualitatively the results obtained with MCMC (not shown).

Identifiability was inspected in different ways. In a frequentist setting, we computed the Fisher Information Matrix (FIM) from the likelihood, defined as

(9)

for a dataset of length , . The FIM is obtained by taking expectations over many time series of length . In practice, we use and simulations; we previously computed the Hessian matrix for very large , with similar results. A parameter set is identifiable whenever the FIM is non-singular (i.e., no zero eigenvalue) because this condition is equivalent to the matrix being positive definite, the FIM being real and symmetric. Positive definiteness, in turn, is essential to having a unique minimum of the negative log-likelihood (Rothenberg et al., 1971). It is therefore possible to quantify identifiability by assessing whether the FIM has some zero eigenvalues. These arguments may appear a bit dull to ecologists, so we give a geometric explanation that could be more intuitive: the FIM is the expected Hessian of , i.e., the expected second derivative of the negative log-likelihood, or in other words, its curvature. If the vector of parameters , where is some subset of parameter space where all maximize the likelihood (e.g., there is a true ridge on the likelihood surface), then both the slope and the curvature of must be zero for along some direction, which is equivalent to having a zero FIM eigenvalue.

Moreover, as classical statistical theory dictates, asymptotically (as the time series length gets large or as we average of over many datasets), the vector of parameters

. We therefore also computed the expected pairwise correlation between the parameters: the expected variance-covariance matrix of the parameters is defined by

so that we get easily the expected pairwise parameter correlation matrix , with each element defined as . This analysis was done here for the FIM, or expected Fisher Information, but as a check we have also performed similar inspections of the variance-covariance matrix derived from the Hessian at the maximum likelihood estimate (which is the observed FIM for one simulation, Viallefont et al., 1998) with similar results.

In a Bayesian setting, we inspected the correlations in the Markov chains for pairs of parameters, which translates into pair posterior distributions of parameters. Parameters whose chains were too positively or negatively correlated were considered to be unidentifiable separately. We also examined whether the prior and posterior distributions overlapped

(Gimenez et al., 2009), which is a way to check that the model is identifiable in practice, in the sense that the likelihood brings enough information to the posteriors.

Results

Identifiability diagnostics

Frequentist analysis
Fisher Information Matrices

We consider here two contrasted cases, for : , we have data on kill rates for the whole time series vs , no data on kill rates for the whole time series.

Figure 2: Distribution of eigenvalues of the FIM, perturbed fixed point case. First row, eigenvalues of the FIM for the model with kill rate data, second row, eigenvalues for the model without kill rate data. Third and fourth rows, eigenvalues for Hilbert matrices of similar size for comparison; firth and sixth rows, eigenvalues for Wishart matrices of similar size for comparison.

We have reported the eigenvalues of the FIM with and without kill rate data for the fixed point in Fig 2 (first and second row, respectively). It is challenging to interpret whether eigenvalues are truly zero (Viallefont et al., 1998; Gimenez et al., 2004), thus the condition number, the ratio between the largest and smallest eigenvalue is usually given. As eigenvalues are reported in log-scale in Fig 2, the distance between the smallest and largest eigenvalue on the graph represents the condition number. To aid the interpretation, we have plotted in Fig 2 the eigenvalue spectrum of two types of matrices of similar size: Hilbert matrices (with elements ), notorious for being on the border of invertibility, and random Wishart matrices, generated as with a matrix with i.i.d. entries. The structure of the FIM being block-diagonal, which helps invertibility, we also constructed two equivalent block-diagonal Hilbert and Wishart matrices, whose elements are set to zero whenever FIM elements are zero. The comparison of the FIMs to these matrices show that while the FIMs with and without kill rate data have lower condition number that random Wishart matrices (which should come close to what the FIM eigenvalues can look like under ideal circumstances), their condition numbers are also larger than those of ill-behaved Hilbert matrices. FIMs therefore be considered to have non-zero eigenvalues with and without kill rate data. Similar results are presented in Appendix 2 for the noisy limit cycle parameter set.

Variance-Covariance Matrices

The variance-covariance matrices display considerable correlation between pairs of parameters that belong to the same functional form of the model ( and , and , and ), and null to weak correlation between parameters that belong to different functional forms. As we shall show below, these within-functional forms correlations are very strong but not necessarily too problematic. Fig. 3 shows that the within-functions correlations emerge with and without kill rate data, and for both parameter sets, the perturbed fixed point as well as the limit cycle. However, without the kill rate data, and for both parameter sets, additional between-functions correlations between the functional response and prey growth parameters () start to emerge.

(a) With kill rate data, FP
(b) Without kill rate data, FP
(c) With kill rate data, LC
(d) Without kill rate data, LC
Figure 3: Correlation matrix between parameters at the true parameter value (FP: perturbed fixed point; LC: noisy limit cycle)
Likelihood surfaces

Based on the variance-covariance matrix results, we output now the (negative) log-likelihood surfaces for the pairs of parameters that were correlated , and . Very similar results can be obtained using directly the residual sum of squares, without taking into account the residual variances ; we present only the negative log-likelihood in Fig. 4 for consistency and simplicity.

(a) Perturbed fixed point parameter set
(b) Noisy limit cycle parameter set
Figure 4: Likelihood surfaces for pairs of parameters. In columns, pairs of parameters , and ; upper row: with kill rate data (panels A-B-C); bottom row: without kill rate data (panels B-C-D). Both parameter sets are considered.

In spite of the correlation between parameters, we see in Fig. 3(a) that a maximum can be found on the likelihood when kill rate data is provided, as the eigenvalues of the FIM suggested. However, when the kill rate data is not included (D and F in Fig. 3(a)), the prey growth rate and functional response parameter values do not correspond clearly to minima of the negative log-likelihood. Fig. 3(b) shows the same plots for the noisy limit cycles parameter set, where with or without the kill rate data, the minima are very well-defined.

Bayesian analysis
Prior-posterior overlap

For one simulation, we now plot the overlap between prior and posterior distributions in Fig. 4(a). For parameters and , the near-absence of overlap without kill rate data and near-perfect prior-posterior overlap in the case without kill rate data demonstrates the practical unidentifiability of the functional response parameters without kill rate (KR) data, for the perturbed fixed point parameter set.

(a) Perturbed fixed point dataset
(b) Noisy limit cycle dataset
Figure 5: Prior-posterior overlap with kill rate (KR) data (left) and without kill rate data (right). Top row, parameter ; bottom row, parameter . True parameter values are red vertical lines.

These results hold for all simulations, as exemplified by the reported estimates for the perturbed fixed point case without KR data: the posterior mass concentrates in this case at the mode of the prior (see section Data availability scenarios). Only in the noisy LC case can we identify correctly in practice and (See Fig. 4(b)) without the use of KR data. And, even though all panels of Fig. 4(b) demonstrate low prior-posterior overlap, the posterior mass for both and parameters is much more narrowly distributed with KR data included, so that the final precision of the estimates is much improved with KR data no matter the parameter set considered.

Pairwise correlations in the joint posterior distribution
(a) With kill rate data
(b) Without kill rate data
Figure 6: Biplots of MCMC samples per pair of parameters, perturbed fixed point (FP) parameter set.

The correlation in the posteriors mirrors the frequentist results obtained on the variance-covariance matrix, and is shown in Fig. 6 for the fixed point parameter set. Pairs of parameters belonging to the same functional form are not estimated as independent. The correlation between parameters typically becomes worse without the kill rate data (Fig. 5(b)). Very similar results are presented in Appendix 2 for the noisy limit cycle parameter set.

Estimation of curves rather than single parameters
Figure 7: Estimated average functional response with vs without correlation between parameters, with (top row) and without (bottom row) kill rate data. The red line is the true, simulated average functional response (note that Gaussian noise is added to that curve in the stochastic predator-prey model, eq. 3).
Figure 8: Estimated average prey growth rate-density curve with vs without correlation between parameters, with (top row) and without (bottom row) kill rate data. The true curve is drawn in red.

Examination of the shape of the growth rate-density curves for both predator and prey, as well as the average functional response, reveals interesting model properties. Namely, the correlation between parameters belonging to the same function that is estimated serves as to maximize the precision of the estimation of the function as a whole. In other words, we estimate well the function and the parameter pair but not the and parameters separately. This can be seen by overlay of the functions for each iteration in the Markov chain, with vs without the correlations between parameters (Fig. 7). We see that the cloud of lines that corresponds to the posterior distribution of the average functional response itself is much thinner when the correlations are included (when the functional response can be estimated, that is, with kill rate data in Fig. 7). The same holds for the prey growth rate-density curve in Fig. 8, which can be well estimated with or without the kill rate data. Therefore, even though correlations in the posterior between parameters seem to make the model unidentifiable, the opposite is actually true: correlation between parameters are instrumental to make the functional forms of the model more precise. A reparameterization proposed in the next section allows to get rid of most of those correlations and to estimate the parameters independently, but does not markedly improve the estimation of the functional curves when compared to the present parameterization with posterior correlations included.

We show here plots for the prey density-dependence curve but similar results can be obtained for the predator density-dependence curve. Similar plots are also provided in Appendix 2 regarding the noisy LC parameter set.

Re-parameterization of the model

In this section, we consider a re-parameterized form of the model given by a prey carrying capacity as in the classical Beverton-Holt model, and a transformation of to following the same logic. The functional response has for its part being changed from (classic Monod or Michaelis-Menten formulation) to the more mechanistic Holling parameterization . We then arrive at the equations:

(10)
(11)

The correlations between pairs of parameters have markedly decreased with the reparameterized model, as exemplified by the correlation matrices derived from the new Fisher Information Matrix at the true parameter value shown in Fig. 9 (similar results have been obtained in a Bayesian setting, Appendix 3).

(a) With kill rate data, FP
(b) Without kill rate data, FP
(c) With kill rate data, LC
(d) Without kill rate data, LC
Figure 9: Correlation matrix between parameters at the true parameter value (FP: perturbed fixed point; LC: noisy limit cycle); here.

In this case, the reparameterization in terms of carrying capacities worked well, especially for the predator, and to some degree for the Holling-type parameterization of the functional response. The decrease in identifiability without the kill rate data is also robust to the change in parameterization (Appendix 3). However, it is important to realize that while improving the identifiability of individual parameters, such reparameterization does not markedly improve the quality of the estimation of the functional response and growth rate curve per se (see Appendix 3). We therefore have no reason to strongly favor the reparameterized form over the original form of the model given by eqs. 1-2.

Data availability scenarios

Our numerical experiment crosses 3 time series length and 3 data availability scenarios (fraction of time units for which the kill rate is available, 0%, 25% and 100%). Figs. 10 and 11 illustrate the bias and precision for the and parameters of the functional response, for the fixed point (FP) parameter set. The previous Bayesian results using a single simulated dataset are confirmed: total absence of kill rate data prevents proper estimation of parameters and (for all 100 simulated FP datasets). However, adding kill rate data for only one fourth of the duration of the time series () gives estimates of that are almost as precise as with kill rate data for the whole of the time series. The half-saturation constant benefits as well from the addition of kill rate data, though less than when looking at the marginal posterior distribution. But it should be kept it mind that as these parameters co-vary strongly, if is better estimated, must be too.

Results regarding and for the noisy LC parameter set are shown in Appendix 2. For the noisy limit cycle (LC) parameter set, estimation of and is possible in absence of kill rate data but adding kill rate data greatly decreases bias and increases precision. Especially, adding kill rate data for one fourth of the time series in the case yields bias and precision that are barely different from . The benefits are still very great for , and then for even with both density and kill rate data, there remains considerable uncertainty about the parameter values, with absence of identifiability for some simulations.

Figure 10: Posterior probability densities for the parameter , for 100 simulations (each black line is one posterior probability density). The vertical red line materializes the true parameter value.
Figure 11: Posterior probability densities for the parameter , for 100 simulations (each black line is one posterior probability density). The vertical red line materializes the true parameter value.

Discussion

We have simulated data from a stochastic predator-prey model with a stochastic functional response, for parameters where the deterministic skeleton of the model would predict either convergence towards a stable fixed point or a limit cycle. These dynamic models, when environmental and interaction stochasticity were included, exhibited considerable variation in both population densities and the functional response, as could be expected from theory (Nisbet and Gurney, 1982; Greenman and Benton, 2005).

We then studied the identifiability of parameters in both theory (sensu Fisher Information Matrix, averaged over realizations of a stochastic process) and practice (for a given dataset of finite time series length). Our model examined the relative contributions to identifiability of density time series and kill rate data. We have found that kill rate data is essential to identifiability in the case of a fixed point equilibrium perturbed by stochasticity, and that it greatly improves bias and precision of functional response parameters in the case of a noisy limit cycle (formally, a noisy invariant loop). Finally, we have shown that a small amount of data on kill rates might go a long way towards improving the practical identifiability for both parameter sets. Now we comment these results in greater detail, as well as their possible extensions and limitations.

Identifiability in theory and in practice of the model with kill rate data

The model is always identifiable in theory (often called ‘structural identifiability’), in the sense of having a nonsingular Fisher Information Matrix. It is also always identifiable in practice when kill rate data is included, in the sense of a low prior-posterior overlap when starting with vague priors.

However, for the initial parameterization of the model, pairs of parameters belonging to the same functional forms of the models are highly correlated. These pairs are, respectively, , , and . These substantial correlations occur in a congruent manner in both the (joint) posterior distributions of Bayesian estimates and the variance-covariance matrix in a frequentist setting.

Reparameterizing using carrying capacities substantially decreased parameter correlations within density-dependence functions. Regarding the functional response (FR), whilst Bolker (2008, chapter 6, p. 200), who fitted a Holling-type FR suggested that a Michaelis-Menten formulation of the FR may improve identifiability compared to the Holling formulation (our original guess as well), we did not find this here: correlation seemed slightly lower for the pair. However, some correlation between the two parameters always seemed to remain. One take home message from our work is that correlations between parameters are more or less problematic depending on which model components they affect; a key question to ask when noticing correlations between estimated parameters is therefore whether they occur between functional forms of interest - which is main hassle - or within functional forms, which does not really pose a problem so long as there is still a well-defined maximum in the likelihood.

Kao and Eisenberg (2018) found, in an epidemiological model fitted on short incidence time series, that in spite of the unidentifiability of individual transmission parameters (that were very strongly correlated), they could estimate the basic reproduction number successfully. Our setup here is less extreme as even with the correlations, the MLE is well-defined (for at least), thus it is doable to estimate the parameters individually with kill rate data. Still, their findings resonate with ours: while the values for the and parameters should be interpreted individually with great caution due to the correlations in the estimates, the functional response as a whole can always be estimated with precision when counts are complemented with kill rate data. In other words, and form a practically identifiable parameter combination, sensu Eisenberg and Hayashi (2014). Whether the functional forms (like here) or aggregate properties (as in Kao and Eisenberg, 2018) can be well-estimated could therefore be given primacy over estimation of individual parameters, that are sometimes of more limited ecological significance.

Identifiability without kill rate data

The absence of kill rate data substantially compromises the practical identifiability of the functional response (and its component parameters) whenever there is a fixed point equilibrium, but theoretical identifiability (sensu FIM) is maintained. This is true for simulations with a perturbed fixed point, showing overdamped convergence to a point equilibrium in absence of noise, which results in fluctuations around the fixed point (that are non-cyclic here). In the main text, we report simulations with length but additional simulations using show the same results, hence this is not just a question of time series length, but of the likelihood being too flat.

These results are true for the model and parameters that we consider in the main text, but also other model types such as a discrete-time version of the Rosenzweig-MacArthur (RMA) model (Supplement S1), as well as models with more temporally variable functional reponses (Supplement S2). Turchin and Ellner (2000), using nonlinear forecasting to fit a similar stochastic predator-prey model to ours (whose likelihood had no closed form, because of continuous time), also reported a lesser identification of parameters for stable fixed points. We are therefore confident that our main result has some degree of generality. However, we concede that quasi-cycles (Nisbet and Gurney, 1982) with more discernible cyclic fluctuations (e.g., closer to the discrete Hopf bifurcation, Wiesenfeld, 1985; Neiman et al., 1997) should logically be easier to identify without kill rate data. In fact, quasi-cycles might, in many cases, be difficult to distinguish from limit cycles (but see Louca and Doebeli, 2014). Surprisingly, our investigations of the RMA model in Supplement S1 do not quite appear to support that theory, since the functional response was poorly estimated without kill rate data for the QC parameter set. But this could be due to the specifics of the simulation or discrete-time model, an investigation with more emphasis on quasi-cycles would be needed for a definite answer.

We also tried models with a stochastic noise on the functional response but without kill rate data, and these did not perform better (typically worse) than the models with a deterministic functional response. This is likely because the model has then no information available to dissociate functional response signal vs noise, and the latent state is in this case poorly estimated - a problem that is reminiscent of state-space models without knowledge of observation error (Auger-Méthé et al., 2016).

We used prior-posterior overlap to quantify practical identifiability, together with inspection of likelihood surfaces for pairs of parameters. Another useful tool is the profile likelihood (Bolker, 2008; Raue et al., 2009; Eisenberg and Hayashi, 2014), which would likely yield similar results. However, as explained in Eisenberg and Hayashi (2014), likelihood profiles can sometimes be ambiguous when combinations of multiple parameters are at play (Raue et al., 2009). They are also very computationally intensive, since for all values of one parameter, optimization has to be performed for all other parameters. We have therefore preferred to evaluate the practical identifiability with the prior-posterior overlap, which is available as a by-product of any Bayesian estimation.

Relative performances of the data combinations

We found that adding kill rate data on just one fourth of the time series greatly increased the performances with respect to no kill rate data, in the sense of decreasing the bias and increasing the precision of the estimators. In some cases (noisy limit cycle), this kill rate data availability for only one fourth of the time series was almost as good as having it throughout the whole time series. This is remarkable, as this requires the algorithm to estimate the stochastic kill rate as a latent state for 75% of the time series, which is usually a hard task.

Time series in ecology are typically short, from a dozen to often at best a hundred of time points, especially for vertebrates (some fast-living species can have much longer time series), and functional response data are rarely available through the whole time series (e.g., Gilg et al., 2003). While very short time series of time steps seemed beyond saving for our nonlinear predator-prey models, no matter what data was available, the results for and with 25% of functional response data are quite encouraging, especially for cyclic species. For non-cyclic species, and worked well but a time series length presented many unidentifiable simulations for max intake and especially the half-saturation constant . It may be possible, in this case, to improve practical identifiability by complementing the kill rate data with other data types, typically survival and reproduction data (Péron and Koons, 2012; Barraquand and Gimenez, 2019), or more informative priors on the intrinsic growth rates.

Avenues for methodological development

The behaviour of deterministic predator-prey models is typically sensitive to the type of functional response that is used (Fussmann and Blasius, 2005; Aldebert and Stouffer, 2018). Here, we used a type II functional response of the Holling type. Different functional responses, e.g., an equivalent Ivlev-type may yield slightly different results in a deterministic setting (e.g., values yielding convergence to equilibrium or limit cycle), though it is known that stochasticity smoothes over the transition between different dynamical behaviours (Wiesenfeld, 1985; Neiman et al., 1997; Barraquand et al., 2017), making the dynamical behaviour probably less sensitive to the exact functional response formulation in stochastic models.

There would, however, be several ways to extend the present results exploiting different functional response formulations:

  • Simulating and fitting the model with the Ivlev functional response, to check the robustness of our results to the functional form used.

  • Simulating with Holling and fitting with Ivlev, or vice versa, keeping in mind that exact parameter match is impossible because of slightly different shapes.

Moreover, the functional response can depend on additional variables such as predator densities (Skalski and Gilliam, 2001) or multiple species, which could be added as well. Predator-dependence can emerge from the discreteness of time (Abrams and Ginzburg, 2000), and might be worth adding to the model. That said, regarding the addition of multiple species, a simple and elegant way to model implicitly many species is to consider a stochastic, temporally variable half-saturation constant. We have done this in Supplement 2, and such model fitting works very well.

Our models used discrete time. Continuous-time equivalents using stochastic differential equations could be considered (Gilioli et al., 2008, 2012), though two difficulties arise. First, the likelihood of a nonlinear stochastic differential equations cannot usually be written in a closed form, which requires evaluation by more elaborate algorithms. These include nonlinear forecasting using attractor reconstruction methods (Turchin and Ellner, 2000), particle filtering (Ionides et al., 2006) or other likelihood-free methods including Approximate Bayesian Computation, see Fasiolo et al. (2016) for a review. In a stochastic differential equation context, Gilioli et al. (2008, 2012)

discretized their equations with an Euler’s scheme to simplify the inferential problem; in that case, it can be as simple and more transparent with regard to biological assumptions to work with a discrete-time predator-prey model from the start. Second, the model presented here can easily be extended to incorporate environmental variation that is correlated in time. That is typically more difficult in stochastic differential equations, whose noise terms are Brownian (Wiener) processes. At the moment, discrete-time therefore seems a more convenient option in many respects, but continuous-time stochastic models also have clear advantages (e.g., no ordering of events and a good connection to most of the multi-species deterministic theory) that should make them an interesting target for future work.

Finally, we glossed over some of the difficulties in actually measuring kill rates from the number of kills, while taking into account non-replacement of prey items (Vonesh and Bolker, 2005; Fenlon and Faddy, 2006; Rosenbaum and Rall, 2018). While we believe that these observational issues deserve attention, this type of work may be best done as a follow-up on a case study. Moreover, many if not most sources of kill rate data are not numbers of observed kills, but in fact reconstructed kill rates combining observed diet with some maximum or allometrically reconstructed consumption (Christensen and Pauly, 1992; Nielsen, 1999; Gilg et al., 2006); our model framework is ideal for the latter data types and can be readily applied in this context.

In addition to the functional response, it might be worth to add to future models fitted to real data some demographic information, which amounts to formulate integrated population and community models (Besbeas et al., 2002; Abadi et al., 2010; Péron and Koons, 2012). These are also currently under development for predator-prey interactions (Barraquand and Gimenez, 2019), and it is likely that more diverse ways to combine information in dynamic models of species interactions will emerge.

Acknowledgements

FB was supported by Labex COTE (ANR-10-LABX-45) and OG by the French National Research Agency (grant ANR-16-CE02-0007).

1 Stability analysis of the deterministic model and stochastic implications

The model, as specified in log-scale by eqs. 78, can be analyzed by computing the Jacobian for the (non-trivial) fixed equilibrium point. For this, we need first the find the fixed point, which is defined by

(A1)
(A2)

The second equation gives so that . Insert this in eq. A1, and we find

(A3)

with . Unfortunately, this is a transcendental equation so there is not closed form solution. But we can say a few things about the equilibrium by the study of the function

(A4)

for . The function is zero at the fixed point. We have and , thus the function goes through zero at least once for all positive values (there is at least one fixed point). The derivative of exists and is defined by

(A5)

so for , (which follows from the definition of parameters), there is a monotonic decrease of the function from to . The function therefore passes only once through zero, which means that there is a unique strictly positive equilibrium point . This equilibrium can be found numerically by iterative convergence, we have done this using the R package ‘rootSolve’.

Derivation of the Jacobian

The equations 78 can be rewritten

(A6)
(A7)

with , and

(A8)
(A9)

The Jacobian of the log-transformed system is by definition

(A10)

We introduce modified functions for convenience

(A11)
(A12)

which leads to a Jacobian

(A13)

with and . Doing so, we obtain

(A14)

which can be rewritten using the relationships at equilibrium

(A15)

The Jacobian is then computed using the numerically found , in the case of the fixed point parameter set and in the case of the limit cycle parameter set .

Using our two parameter sets, we found the following eigenvalues for the Jacobian, reported in Table 2.

Attractor
Fixed point 0.44 0.16 0.47
Limit cycle 0.85 0.53 1.004
Table 2: Eigenvalues of the Jacobian for both parameter sets (Fixed point, ; limit cycle, . The values in columns are the real part, the imaginary part, and the modulus of the eigenvalues. A modulus below 1 indicates convergence of the fixed point in absence of stochasticity, a modulus above 1 indicates here an escape of the trajectory towards the invariant loop.

2 Complementary results on the limit cycle parameter set

2.1 Fisher Information Matrix

Figure 12: Distribution of eigenvalues of the FIM, noisy limit cycle case. First row, eigenvalues of the FIM for the model with kill rate data, second row, eigenvalues for the model without kill rate data.

2.2 Posterior correlations

(a) With kill rate data
(b) Without kill rate data
Figure 13: Plots of the MCMC samples per pair of parameters, for the noisy limit cycle parameter set.

2.3 Estimation of functional forms

The functional response curve is plotted in Fig. 14 with vs without the correlations between parameters C and D. In the next Fig. 15, results are shown for the prey density-dependence curve. These are similar to the main text results, except that the functional response can now be estimated without kill rate data for .

Figure 14: Average functional response with vs without correlation between parameters, with (top row) and without (bottom row) kill rate data. The red line is the true, simulated average functional response.
Figure 15: Prey growth rate-density curve with vs without correlation between parameters, with (top row) and without (bottom row) kill rate data. The true, simulated curve is drawn in red.

We present here plots for the prey density-dependence curve but similar results can be obtained for the predator density-dependence curve.

2.4 Data availability scenarios

The plots are similar to those of the main text, examining the posterior probability densities of and , except that here we do this for the noisy limit cycle (LC) parameter set.

Figure 16: Posterior probability densities for the parameter , for 100 simulations (each grey line is one posterior probability density).The vertical red line materializes the true parameter value.
Figure 17: Posterior probability densities for the parameter , for 100 simulations (each grey line is one posterior probability density). The vertical red line materializes the true parameter value.

3 Complementary results on the reparameterized model

Here we report the results of the Bayesian estimation of the reparameterized model of eqs. 1011 of the main text. Fig. 18 illustrates that and suffer from the same identification issues as and without kill rate data, for the fixed point parameter set (identical parameters and T = 100 simulation as used in the main text).

Figure 18: Perturbed fixed point dataset. Prior-posterior overlap for the reparameterized predator-prey model.

Although are slightly less correlated than , the average functional response curve as a whole is not better estimated (Fig. 19).

Figure 19: Average functional response with vs without correlation between parameters, with (top row) and without (bottom row) kill rate data. The red line is the true, simulated average functional response.

We can see that the prey growth rate – density curve is, on the other hand, estimated without notable correlations: independent permutations of the temporal order of parameters on chains does not affect the width of this curve anymore (Fig. 20). This is due to the use of a carrying capacity rather than a density-dependence coefficient , which removes posterior correlations (see also main text for the likelihood-based results). However, as we note in the main text, the parameterization does not yield an intrinsically more precise prey growth rate – density curve than the parameterization, when correlations between parameters in the posteriors are allowed.

Figure 20: Prey growth rate-density curve with vs without correlation between parameters, with (top row) and without (bottom row) kill rate data. The true, simulated curve is drawn in red.

References

  • Abadi et al. (2010) Abadi, F., O. Gimenez, R. Arlettaz, and M. Schaub. 2010. An assessment of integrated population models: bias, accuracy, and violation of the assumption of independence. Ecology 91:7–14.
  • Abrams (1982) Abrams, P. 1982. Functional responses of optimal foragers. American Naturalist pages 382–390.
  • Abrams (2010) ———. 2010. Quantitative descriptions of resource choice in ecological models. Population Ecology 52:47–58.
  • Abrams and Ginzburg (2000) Abrams, P., and L. Ginzburg. 2000. The nature of predation: prey dependent, ratio dependent or neither? Trends in Ecology & Evolution 15:337–341.
  • Aldebert and Stouffer (2018) Aldebert, C., and D. B. Stouffer. 2018. Community dynamics and sensitivity to model structure: towards a probabilistic view of process-based model predictions. Journal of the Royal Society Interface 15:20180741.
  • Auger-Méthé et al. (2016) Auger-Méthé, M., C. Field, C. M. Albertsen, A. E. Derocher, M. A. Lewis, I. D. Jonsen, and J. Mills Flemming. 2016. State-space models’ dirty little secrets: even simple linear Gaussian models can have estimation problems. Scientific Reports 6:26677. doi:10.1038/srep26677.
  • Barraquand and Gimenez (2019) Barraquand, F., and O. Gimenez. 2019. Integrating multiple data sources to fit matrix population models for interacting species. bioRxiv doi:10.1101/594986.
  • Barraquand and Nielsen (2018) Barraquand, F., and O. K. Nielsen. 2018. Predator-prey feedback in a gyrfalcon-ptarmigan system? Ecology and Evolution 8:12425–12434. doi:10.1002/ece3.4563.
  • Barraquand et al. (2017) Barraquand, F., S. Louca, K. C. Abbott, C. A. Cobbold, F. Cordoleani, D. L. DeAngelis, B. D. Elderd, J. W. Fox, P. Greenwood, F. M. Hilker, et al. 2017. Moving forward in circles: challenges and opportunities in modelling population cycles. Ecology Letters 20:1074–1092.
  • Besbeas et al. (2002) Besbeas, P., S. N. Freeman, B. J. Morgan, and E. A. Catchpole. 2002. Integrating mark–recapture–recovery and census data to estimate animal abundance and demographic parameters. Biometrics 58:540–547.
  • Bolker (2008) Bolker, B. 2008. Ecological models and data in R. Princeton University Press.
  • Caswell (2001) Caswell, H. 2001. Matrix population models. Sinauer Associates.
  • Catchpole and Morgan (1997) Catchpole, E. A., and B. J. Morgan. 1997. Detecting parameter redundancy. Biometrika 84:187–196.
  • Christensen and Pauly (1992) Christensen, V., and D. Pauly. 1992. ECOPATH II - a software for balancing steady-state ecosystem models and calculating network characteristics. Ecological modelling 61:169–185.
  • Cole et al. (2010) Cole, D. J., B. J. Morgan, and D. Titterington. 2010. Determining the parametric structure of models. Mathematical biosciences 228:16–30.
  • Cook et al. (2015) Cook, R. M., S. J. Holmes, and R. J. Fryer. 2015. Grey seal predation impairs recovery of an over-exploited fish stock. Journal of Applied Ecology 52:969–979.
  • Curtsdotter et al. (2019) Curtsdotter, A., H. T. Banks, J. E. Banks, M. Jonsson, T. Jonsson, A. N. Laubmeier, M. Traugott, and R. Bommarco. 2019. Ecosystem function in predator-prey food webs – confronting dynamic models with empirical data. Journal of Animal Ecology 88:196–210. doi:10.1111/1365-2656.12892.
  • de Roos et al. (1991) de Roos, A., E. McCauley, and W. Wilson. 1991. Mobility versus density-limited predator–prey dynamics on different spatial scales. Proceedings: Biological Sciences 246:117–122.
  • DeLong et al. (2018) DeLong, J. P., T. C. Hanley, J. P. Gibert, L. M. Puth, and D. M. Post. 2018. Life history traits and functional processes generate multiple pathways to ecological stability. Ecology 99:5–12.
  • Eisenberg and Hayashi (2014) Eisenberg, M. C., and M. A. Hayashi. 2014. Determining identifiable parameter combinations using subset profiling. Mathematical biosciences 256:116–126.
  • Fasiolo et al. (2016) Fasiolo, M., N. Pya, S. N. Wood, et al. 2016. A comparison of inferential methods for highly nonlinear state space models in ecology and epidemiology. Statistical Science 31:96–118.
  • Fenlon and Faddy (2006) Fenlon, J. S., and M. J. Faddy. 2006. Modelling predation in functional response. Ecological Modelling 198:154–162.
  • Ferguson et al. (2018) Ferguson, J. M., J. B. Hopkins III, and B. H. Witteveen. 2018. Integrating abundance and diet data to improve inferences of food web dynamics. Methods in Ecology and Evolution 9:1581–1591.
  • Fussmann and Blasius (2005) Fussmann, G., and B. Blasius. 2005. Community response to enrichment is highly sensitive to model structure. Biology Letters 1:9–12.
  • Gilg et al. (2003) Gilg, O., I. Hanski, and B. Sittler. 2003. Cyclic dynamics in a simple vertebrate predator-prey community. Science 301:866–868.
  • Gilg et al. (2006) Gilg, O., B. Sittler, B. Sabard, A. Hurstel, R. Sané, P. Delattre, and I. Hanski. 2006. Functional and numerical responses of four lemming predators in high arctic Greenland. Oikos 113:193–216.
  • Gilioli et al. (2008) Gilioli, G., S. Pasquali, and F. Ruggeri. 2008. Bayesian inference for functional response in a stochastic predator–prey system. Bulletin of Mathematical Biology 70:358–381.
  • Gilioli et al. (2012) ———. 2012. Nonlinear functional response parameter estimation in a stochastic predator-prey model. Math. Biosci. Eng 9:75–96.
  • Gimenez et al. (2004) Gimenez, O., A. Viallefont, E. A. Catchpole, R. Choquet, and B. J. Morgan. 2004. Methods for investigating parameter redundancy. Animal Biodiversity and Conservation 27:561–572.
  • Gimenez et al. (2009) Gimenez, O., B. J. Morgan, and S. P. Brooks. 2009. Weak identifiability in models for mark-recapture-recovery data. In Modeling demographic processes in marked populations, pages 1055–1067. Springer.
  • Greenman and Benton (2005) Greenman, J., and T. Benton. 2005. The impact of environmental fluctuations on structured discrete time population models: resonance, synchrony and threshold behaviour. Theoretical population biology 68:217–235.
  • Holling (1959) Holling, C. 1959. The components of predation as revealed by a study of small mammal predation of the European pine sawfly. Canadian Entomologist 91:293–320.
  • Hosseini (2003) Hosseini, P. 2003. How localized consumption stabilizes predator-prey systems with finite frequency of mixing. American Naturalist 161:567–585.
  • Ionides et al. (2006) Ionides, E., C. Bretó, and A. King. 2006. Inference for nonlinear dynamical systems. Proceedings of the National Academy of Sciences 103:18438–18443.
  • Ives (1995) Ives, A. R. 1995. Predicting the response of populations to environmental change. Ecology 76:926–941.
  • Ives et al. (2008) Ives, A. R., Á. Einarsson, V. A. Jansen, and A. Gardarsson. 2008. High-amplitude fluctuations and alternative dynamical states of midges in Lake Myvatn. Nature 452:84–87.
  • Jost and Arditi (2001) Jost, C., and R. Arditi. 2001. From pattern to process: identifying predator–prey models from time-series data. Population Ecology 43:229–243.
  • Jost and Ellner (2000) Jost, C., and S. P. Ellner. 2000. Testing for predator dependence in predator-prey dynamics: a non-parametric approach. Proceedings of the Royal Society of London B: Biological Sciences 267:1611–1620.
  • Kalinkat et al. (2013) Kalinkat, G., F. D. Schneider, C. Digel, C. Guill, B. C. Rall, and U. Brose. 2013. Body masses, functional responses and predator–prey stability. Ecology letters 16:1126–1134.
  • Kao and Eisenberg (2018) Kao, Y.-H., and M. C. Eisenberg. 2018. Practical unidentifiability of a simple vector-borne disease model: Implications for parameter estimation and intervention assessment. Epidemics 25:89 – 100. doi:https://doi.org/10.1016/j.epidem.2018.05.010.
  • Karban and De Valpine (2010) Karban, R., and P. De Valpine. 2010. Population dynamics of an arctiid caterpillar–tachinid parasitoid system using state-space models. Journal of Animal Ecology 79:650–661.
  • Lande et al. (2003) Lande, R., S. Engen, and B. Sæther. 2003. Stochastic population dynamics in ecology and conservation. Oxford University Press, USA.
  • Leslie (1948) Leslie, P. 1948. Some further notes on the use of matrices in population mathematics. Biometrika 35:213–245.
  • Leslie and Gower (1960) Leslie, P., and J. Gower. 1960. The properties of a stochastic model for the predator-prey type of interaction between two species. Biometrika 47:219–234.
  • Little et al. (2010) Little, M. P., W. F. Heidenreich, and G. Li. 2010. Parameter identifiability and redundancy: theoretical considerations. PloS one 5:e8915.
  • Louca and Doebeli (2014) Louca, S., and M. Doebeli. 2014. Distinguishing intrinsic limit cycles from forced oscillations in ecological time series. Theoretical ecology 7:381–390.
  • Mac Nally et al. (2010) Mac Nally, R., J. R. Thomson, W. J. Kimmerer, F. Feyrer, K. B. Newman, A. Sih, W. A. Bennett, L. Brown, E. Fleishman, S. D. Culberson, et al. 2010.

    Analysis of pelagic species decline in the upper San Francisco Estuary using multivariate autoregressive modeling (MAR).

    Ecological Applications 20:1417–1430.
  • Maunder (2004) Maunder, M. N. 2004. Population viability analysis based on combining bayesian, integrated, and hierarchical analyses. Acta Oecologica 26:85–94.
  • May (1973) May, R. 1973. Stability and complexity in model ecosystems. Princeton University Press, Princeton, USA.
  • Murrell (2005) Murrell, D. 2005. Local spatial structure and predator-prey dynamics: Counterintuitive effects of prey enrichment. American Naturalist 166:354–367.
  • Mutshinda et al. (2009) Mutshinda, C. M., R. B. O’Hara, and I. P. Woiwod. 2009. What drives community dynamics? Proceedings of the Royal Society B: Biological Sciences 276:2923–2929.
  • Neiman et al. (1997) Neiman, A., P. I. Saparin, and L. Stone. 1997. Coherence resonance at noisy precursors of bifurcations in nonlinear dynamical systems. Physical Review E 56:270.
  • Nielsen (1999) Nielsen, Ó. 1999. Gyrfalcon predation on ptarmigan: numerical and functional responses. Journal of Animal Ecology 68:1034–1050.
  • Nisbet and Gurney (1982) Nisbet, R., and W. Gurney. 1982. Modelling fluctuating populations. John Wiley & Sons.
  • Niu et al. (2014) Niu, S., Y. Luo, M. C. Dietze, T. F. Keenan, Z. Shi, J. Li, and F. S. C. III. 2014. The role of data assimilation in predictive ecology. Ecosphere 5:art65.
  • Péron and Koons (2012) Péron, G., and D. N. Koons. 2012. Integrated modeling of communities: parasitism, competition, and demographic synchrony in sympatric ducks. Ecology 93:2456–2464.
  • Planque et al. (2014) Planque, B., U. Lindstrøm, and S. Subbey. 2014. Non-deterministic modelling of food-web dynamics. PloS one 9:e108243.
  • Plummer et al. (2003) Plummer, M., et al. 2003. Jags: A program for analysis of bayesian graphical models using gibbs sampling. In Proceedings of the 3rd international workshop on distributed statistical computing, volume 124. Vienna, Austria.
  • Pritchard et al. (2017) Pritchard, D. W., R. A. Paterson, H. C. Bovy, and D. Barrios-O’Neill. 2017. Frair: an r package for fitting and comparing consumer functional responses. Methods in Ecology and Evolution 8:1528–1534.
  • Rall et al. (2012) Rall, B. C., U. Brose, M. Hartvig, G. Kalinkat, F. Schwarzmüller, O. Vucic-Pestic, and O. L. Petchey. 2012. Universal temperature and body-mass scaling of feeding rates. Philosophical Transactions of the Royal Society B: Biological Sciences 367:2923–2934.
  • Raue et al. (2009) Raue, A., C. Kreutz, T. Maiwald, J. Bachmann, M. Schilling, U. Klingmüller, and J. Timmer. 2009. Structural and practical identifiability analysis of partially observed dynamical models by exploiting the profile likelihood. Bioinformatics 25:1923–1929.
  • Rosenbaum and Rall (2018) Rosenbaum, B., and B. C. Rall. 2018. Fitting functional responses: Direct parameter estimation by simulating differential equations. Methods in Ecology and Evolution 9:2076–2090.
  • Rosenbaum et al. (2019) Rosenbaum, B., M. Raatz, G. Weithoff, G. F. Fussmann, and U. Gaedke. 2019. Estimating parameters from multiple time series of population dynamics using bayesian inference. Frontiers in Ecology and Evolution 6:234. doi:10.3389/fevo.2018.00234.
  • Rothenberg et al. (1971) Rothenberg, T. J., et al. 1971.

    Identification in parametric models.

    Econometrica 39:577–591.
  • Schaub and Abadi (2011) Schaub, M., and F. Abadi. 2011. Integrated population models: a novel analysis framework for deeper insights into population dynamics. Journal of Ornithology 152:227–237.
  • Skalski and Gilliam (2001) Skalski, G., and J. Gilliam. 2001. Functional responses with predator interference: viable alternatives to the Holling type II model. Ecology 82:3083–3092.
  • Smout et al. (2010) Smout, S., C. Asseburg, J. Matthiopoulos, C. Fernández, S. Redpath, S. Thirgood, and J. Harwood. 2010. The functional response of a generalist predator. PloS one 5:e10761.
  • Solomon (1949) Solomon, M. 1949. The natural control of animal populations. Journal of Animal Ecology 18:1–35.
  • Stouffer (2019) Stouffer, D. B. 2019. All ecological models are wrong, but some are useful. Journal of Animal Ecology 88:192–195.
  • Tanner (1975) Tanner, J. 1975. The stability and the intrinsic growth rates of prey and predator populations. Ecology 56:855–867.
  • Turchin and Ellner (2000) Turchin, P., and S. Ellner. 2000. Living on the edge of chaos: population dynamics of Fennoscandian voles. Ecology 81:3099–3116.
  • Turchin and Hanski (1997) Turchin, P., and I. Hanski. 1997. An empirically based model for latitudinal gradient in vole population dynamics. American Naturalist 149:842–874.
  • Viallefont et al. (1998) Viallefont, A., J.-D. Lebreton, A.-M. Reboulet, and G. Gory. 1998. Parameter identifiability and model selection in capture-recapture models: A numerical approach. Biometrical Journal 40:313–325.
  • Vonesh and Bolker (2005) Vonesh, J. R., and B. M. Bolker. 2005. Compensatory larval responses shift trade-offs associated with predator-induced hatching plasticity. Ecology 86:1580–1591.
  • Weide et al. (2018) Weide, V., M. C. Varriale, and F. M. Hilker. 2018. Hydra effect and paradox of enrichment in discrete-time predator-prey models. Mathematical Biosciences doi:https://doi.org/10.1016/j.mbs.2018.12.010.
  • Wiesenfeld (1985) Wiesenfeld, K. 1985. Noisy precursors of nonlinear instabilities. Journal of Statistical Physics 38:1071–1097.
  • Wilson (1998) Wilson, W. 1998. Resolving discrepancies between deterministic population models and individual-based simulations. American Naturalist 151:116–134.

1 Estimation and identifiability of a Rosenzweig-MacArthur model

Due to the general congruence between Bayesian estimation through MCMC and ML estimation demonstrated in the main text, we present here only a Bayesian analysis of the discrete-time RMA model.

The stochastic model is defined as

(S1)
(S2)

which translates into

(S3)
(S4)

We have used the following parameter notations and parameter values

Name Meaning FP QC LC
Prey growth rate 2 2 1.8
DD coefficient 1 1 1
Predator mortality 0.2 0.2 0.7
Conversion efficiency 0.1 0.1 0.1
Max kill rate 2.5 2.5 10
Half-saturation constant 1 0.5 0.6
Prey growth rate noise 0.05 0.05 0.05
Pred. growth rate noise 0.05 0.05 0.05
Functional response noise 0.05 0.05 0.05
Table S1: Notation, meaning and values of the parameters for three parameter sets: FP, fixed point (in the deterministic model); QC, quasi-cycles, fixed point with strongly oscillatory dynamics in the presence of stochasticity and LC, noisy limit cycle.

We added a parameter set leading to quasi-cycles (i.e., an excited stable focus, Nisbet and Gurney, 1982, with strongly oscillatory decay to equilibrium in the absence of noise) that is an interesting intermediate case between a noisy limit cycle and a perturbed fixed point (the latter having near-monotonic decay to equilibrium in the absence of noise).

Figure S1: System dynamics for 100 timesteps (1 simulation), after transients (100 timesteps). In panel (A), we show densities of prey in blue and predator in red , for the perturbed fixed point (FP) case, (B) corresponding trajectories in phase plane and (C) functional response: kill rate of individual predator as a function of prey density. In the second and third rows, identical panels (D-E-F) for the quasi-cycles (QC) parameter set and (G-H-I) for the noisy limit cycle (LC).

Due to the presence of important transients in this model, we have removed the first 100 timesteps and done the estimation on the following 100. The correlations between pairs of parameters belonging to the same function was similar to that shown in the main text (omitted).

Figure S2: Prior-posterior overlap with kill rate (KR) data (left) and without kill rate data (right). RMA model, FP parameter set. Top row, parameter ; bottom row, parameter . True parameter values are red vertical lines.
Figure S3: Prior-posterior overlap with kill rate data (left) and without kill rate data (right). RMA model, QC parameter set. Top row, parameter ; bottom row, parameter . True parameter values are red vertical lines.
Figure S4: Prior-posterior overlap with kill rate data (left) and without kill rate data (right). RMA model, LC parameter set. Top row, parameter ; bottom row, parameter . True parameter values are red vertical lines.

The posterior-prior overlap curves are reported in the followings Figs. S2S4. For the perturbed fixed point parameter set, the estimation of parameters and is barely feasible without kill rate data while it is very good with kill rate data (Fig. S2). The findings are similar for the QC parameter set (Fig. S3), even though surprisingly - since we are closer to cyclic dynamics - the identifiability seems even worse in this case (but we caution this is only one simulation). The limit cycle parameter set leads to a less good identification of parameters and in absence of kill rate data than for the model of the main text. Indeed, the right order of magnitude is found for but is a bit misleading (Fig. S4), and in absence of kill rate data, the posterior densities are a bit spiky even though the chains converge in the sense of .

Therefore, our results with the RMA model confirm that a much more precise and less biased estimation of parameters is possible once the data on kill rates is added to the model.

2 Stochastic functional responses with larger noise

In this appendix, we model cases where the functional response exhibit large, seemingly stochastic temporal variation. This is done by introducting temporal variation in the parameters and of the functional response . To data simulated with temporally variable and , we fit:

  • A model with a Gaussian noise on the kill rate (main texts models)

  • A model without kill rate data

  • The simulated model

Therefore, this supplement also serves as an evaluation of the robustness of the results presented in the main text, since we fit the Gaussian FR model to stochastic data provided by a more complex and noisy simulation model.

2.1 Temporal variation in

To make for a realistic point cloud in the functional response, we assume that there is a ‘hard limit’ (say, ) to the kill rate with some variation below. In this formulation of the model, parameters are equal to that of the fixed point (FP) parameter set considered previously, except that the functional response is now

(S5)

The Beta distribution

allows to make for lifelike-looking functional response point clouds (Fig. S5). We have used the FP parameter set of the main text, with added parameters and . To have commensurate noise on all components (and therefore test the general robustness of our results to increased levels of randomness), we chose a noise on growth rates with variance .

Figure S5: Functional response with temporally variable . The max FR value, , is materialized as a horizontal black line. The blue line is a least-square-fit of the functional response assuming a Gaussian noise on the kill rate (instead of a Beta noise on ); the red line the same fit in JAGS, with mildly informative priors, which performs typically a bit better.
2.1.1 Model with a Gaussian functional response, with kill rate data

JAGS estimates () of a Gaussian functional response based on the simulated FR data with temporally variable yield . Given this, the estimates of the functional response in Table S2 are reasonable, even though has some difficulties to converge (which is to be expected, since the fitted model is not fully appropriate).

Parameter mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff
1.50 0.08 1.37 1.44 1.50 1.56 1.65 1.01 230.00
0.16 0.14 0.00 0.00 0.16 0.26 0.45 1.17 41.00
10.70 3.07 5.65 8.53 10.37 12.51 17.77 1.00 4200.00
2.06 1.22 0.49 1.05 1.72 2.92 4.71 1.01 360.00
0.44 0.11 0.24 0.37 0.44 0.51 0.67 1.00 2100.00
2.50 0.56 1.47 2.06 2.49 2.97 3.43 1.01 350.00
0.41 0.06 0.31 0.37 0.41 0.45 0.54 1.00 4900.00
0.59 0.09 0.44 0.52 0.58 0.64 0.78 1.00 6000.00
6.34 0.93 4.65 5.68 6.29 6.96 8.25 1.00 3000.00
Table S2: Estimated parameters for the model with Gaussian FR.

We therefore confirm here that:

  • it is possible to provide a reasonable first approximation of complex functional response with a simple Gaussian noise

  • the main text results are robust to higher noise levels

2.1.2 Model with a deterministic functional response, without additional kill rate data

In this case, we confirm logically the main text results, the functional response is not identifiable and the chains do not converge for parameters and which are informed only by the priors.

Parameter mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff
0.00 0.00 0.00 0.00 0.00 0.00 0.00 1.28 16.00
0.01 0.06 0.00 0.00 0.00 0.00 0.00 2.92 4.00
10.76 3.12 5.68 8.52 10.46 12.63 17.87 1.00 6000.00
1.24 0.94 0.30 0.60 0.91 1.57 3.93 1.00 680.00
0.45 0.11 0.23 0.37 0.44 0.52 0.67 1.00 3100.00
1.94 0.56 1.06 1.52 1.86 2.31 3.16 1.00 750.00
0.41 0.06 0.31 0.37 0.41 0.45 0.55 1.00 5700.00
0.55 0.08 0.41 0.49 0.55 0.60 0.74 1.00 2100.00
Table S3: Estimated parameters for the model with deterministic FR, without kill rate data.
2.1.3 Model with a Beta-distributed max kill rate , with kill rate data

Here, we fit the same model that we simulated, with extra temporal variation on (eq. S5). The model is now quite well estimated, although it should be noted that some chains can have difficulties to converge.

Parameter mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff
2.37 0.11 2.17 2.28 2.36 2.44 2.61 1.00 1100.00
0.48 0.09 0.31 0.42 0.47 0.53 0.66 1.00 1200.00
10.71 3.07 5.55 8.53 10.44 12.52 17.54 1.00 6000.00
1.70 0.42 1.01 1.40 1.65 1.96 2.65 1.00 3600.00
3.70 0.66 2.54 3.23 3.66 4.12 5.11 1.00 6000.00
2.22 1.28 0.49 1.13 1.92 3.22 4.77 1.01 300.00
0.44 0.11 0.23 0.37 0.44 0.52 0.67 1.00 6000.00
2.56 0.57 1.46 2.12 2.58 3.05 3.45 1.01 310.00
0.41 0.06 0.31 0.37 0.40 0.45 0.55 1.00 6000.00
0.59 0.09 0.44 0.53 0.58 0.64 0.78 1.00 2600.00
0.09 0.01 0.07 0.09 0.09 0.10 0.10 1.00 1500.00
Table S4: Estimated parameters for the model with stochastic FR and a temporally Beta-distributed C, with kill rate data. Here is a small SD term that is needed to estimate the model in JAGS with a proper mixing of the chains.

2.2 Temporal variation in the half-saturation constant

Temporal variation in

is added using a Gamma distribution for the half-saturation constant, with mean

and standard deviation

that are then converted into shape (sh) and rate (ra) parameters. We chose and , all other parameters are equal to those of the fixed point simulations, save for the noise variances on the growth rates. To have commensurate noise on all components, we chose .

(S6)

with and .

Figure S6: Functional response with temporally variable . The blue line is a least-square-fit of the functional response assuming a Gaussian noise on the kill rate.
2.2.1 Model with a Gaussian functional response, with kill rate data

Estimates of a least square fit of the Gaussian functional response based on the simulated FR data with temporally variable yield . Given this, the estimates of the functional response in Table S5 are fairly good.

Parameter mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff
2.47 0.12 2.25 2.39 2.47 2.55 2.72 1.00 5800.00
1.05 0.21 0.69 0.90 1.04 1.18 1.51 1.00 3400.00
10.53 3.03 5.55 8.39 10.23 12.31 17.35 1.00 6000.00
1.96 1.18 0.48 1.01 1.67 2.72 4.67 1.01 300.00
0.44 0.11 0.22 0.36 0.44 0.51 0.67 1.00 6000.00
2.45 0.54 1.46 2.03 2.45 2.89 3.40 1.01 280.00
0.41 0.06 0.31 0.37 0.41 0.45 0.55 1.00 6000.00
0.59 0.09 0.45 0.53 0.58 0.64 0.78 1.00 6000.00
2.55 0.41 1.84 2.27 2.53 2.81 3.42 1.00 6000.00
Table S5: Estimated parameters for the model with Gaussian FR.
2.2.2 Model with a deterministic functional response, without additional kill rate data

Here again, we confirm the main text results, the functional response parameters are not identifiable and only informed by the priors.

Parameter mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff
0.00 0.00 0.00 0.00 0.00 0.00 0.00 1.32 15.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 1.67 6.00
10.51 2.99 5.42 8.40 10.24 12.29 17.13 1.00 2300.00
0.63 0.46 0.21 0.37 0.51 0.72 1.87 1.00 1200.00
0.44 0.11 0.23 0.36 0.43 0.51 0.65 1.00 2600.00
1.42 0.41 0.79 1.15 1.36 1.63 2.42 1.00 1400.00
0.41 0.06 0.31 0.37 0.41 0.45 0.55 1.00 6000.00
0.55 0.08 0.41 0.50 0.55 0.60 0.74 1.00 6000.00
Table S6: Estimated parameters for the model with deterministic FR and no kill rate data.
2.2.3 Model with a Gamma-distributed half-saturation , with kill rate data

Here, we fit the same model that we simulated, with extra temporal variation on (eq. S6).

Parameter mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff
2.50 0.00 2.50 2.50 2.50 2.50 2.50 1.00 1.00
10.65 3.04 5.53 8.48 10.36 12.50 17.59 1.00 2600.00
1.91 1.18 0.47 0.99 1.56 2.58 4.70 1.00 520.00
2.07 0.30 1.57 1.86 2.04 2.25 2.74 1.00 1400.00
4.56 0.70 3.39 4.07 4.48 4.97 6.17 1.00 1400.00
0.44 0.11 0.24 0.37 0.44 0.51 0.66 1.00 4600.00
2.42 0.55 1.44 2.00 2.40 2.84 3.42 1.00 530.00
0.41 0.06 0.31 0.37 0.41 0.45 0.55 1.00 6000.00
0.59 0.09 0.45 0.53 0.58 0.64 0.79 1.00 3400.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 1.08 55.00
Table S7: Estimated parameters for the model with stochastic FR and a temporally Gamma-distributed D, with kill rate data. Here is a small SD term that is needed to estimate the model in JAGS with a proper mixing of the chains.

This model is quite well estimated, with excellent convergence; for the parameters chosen here the Gaussian FR and Gamma-distributed half-saturation yield fairly similar parameter values. We therefore believe that in many cases for which the functional response has a deterministic ‘hard’ limit (e.g., the animal cannot possibly kill than items in a certain amount of time), so that is truly constant, a model with

temporally variable will be an interesting model to avoid the assumption of a normally distributed kill rate.