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 capturerecapture 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 communitylevel 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 nonintuitive 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 predatorprey 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 continuoustime ODE models where stochasticity has been added to all parameters (Turchin and Hanski, 1997; Turchin and Ellner, 2000) and in stochastic individualbased predatorprey 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 communitylevel time series of densities and kill rate data, we consider a ground truth dynamical model for predatorprey 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 posteriorprior 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
Predatorprey 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 BevertonHolt function for densitydependence 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 continuoustime counterpart (see Weide et al., 2018 on connecting discretetime to continuoustime predatorprey models). Our model writes
(1) 
(2) 
The roots of this discretetime formulation can be traced back to Leslie (1948); Leslie and Gower (1960) who included a BevertonHolt regulation for the prey and predator to make discretetime models more similar to their continoustime 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 Leslietype parameterization over the RosenzweigMacArthur 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 NicholsonBailey 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 predatorprey model with lognormal 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 datagenerating 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 predatorprey 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 predatorprey 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 predatorprey case where logdensities for both predator and prey are gathered in a logcount 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 discretetime 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 predatorprey 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 predatorprey model
In the simplest case highlighted by our discretetime dynamical systems of the sections above, is the product of the two gaussian pdf for logdensities conditional on past densities and auxiliary information. Using the equations (1) and (2), we obtain the following difference equation on the logscale:
(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 
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 quasicycles sensu Nisbet and Gurney, 1982). This is a crucial example because not all fluctuating predatorprey 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 welldefined 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 discretetime limit cycle (see, e.g., Caswell, 2001, chapter 16). Representative time series of the model for both parameter sets are provided in Fig. 1.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 priorposterior 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 repository^{1}^{1}1https://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 welldefined likelihood function also allows to find point estimates of the parameters through quasiNewton methods (hillclimbing 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 nonsingular (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 loglikelihood (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 loglikelihood, 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 variancecovariance 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 variancecovariance 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.
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 logscale 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 blockdiagonal, which helps invertibility, we also constructed two equivalent blockdiagonal 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 illbehaved Hilbert matrices. FIMs therefore be considered to have nonzero eigenvalues with and without kill rate data. Similar results are presented in Appendix 2 for the noisy limit cycle parameter set.
VarianceCovariance Matrices
The variancecovariance 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 withinfunctional forms correlations are very strong but not necessarily too problematic. Fig. 3 shows that the withinfunctions 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 betweenfunctions correlations between the functional response and prey growth parameters () start to emerge.
Likelihood surfaces
Based on the variancecovariance matrix results, we output now the (negative) loglikelihood 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 loglikelihood in Fig. 4 for consistency and simplicity.
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 loglikelihood. 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 welldefined.
Bayesian analysis
Priorposterior overlap
For one simulation, we now plot the overlap between prior and posterior distributions in Fig. 4(a). For parameters and , the nearabsence of overlap without kill rate data and nearperfect priorposterior 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.
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 priorposterior 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
The correlation in the posteriors mirrors the frequentist results obtained on the variancecovariance 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
Examination of the shape of the growth ratedensity 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 ratedensity 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 densitydependence curve but similar results can be obtained for the predator densitydependence curve. Similar plots are also provided in Appendix 2 regarding the noisy LC parameter set.
Reparameterization of the model
In this section, we consider a reparameterized form of the model given by a prey carrying capacity as in the classical BevertonHolt model, and a transformation of to following the same logic. The functional response has for its part being changed from (classic Monod or MichaelisMenten 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).
In this case, the reparameterization in terms of carrying capacities worked well, especially for the predator, and to some degree for the Hollingtype 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. 12.
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 halfsaturation 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 covary 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.
Discussion
We have simulated data from a stochastic predatorprey 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 priorposterior 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 variancecovariance matrix in a frequentist setting.
Reparameterizing using carrying capacities substantially decreased parameter correlations within densitydependence functions. Regarding the functional response (FR), whilst Bolker (2008, chapter 6, p. 200), who fitted a Hollingtype FR suggested that a MichaelisMenten 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 welldefined 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 welldefined (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 wellestimated 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 noncyclic 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 discretetime version of the RosenzweigMacArthur (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 predatorprey 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 quasicycles (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, quasicycles 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 discretetime model, an investigation with more emphasis on quasicycles 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 statespace models without knowledge of observation error (AugerMéthé et al., 2016).
We used priorposterior 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 priorposterior overlap, which is available as a byproduct 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 fastliving 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 predatorprey 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 noncyclic species, and worked well but a time series length presented many unidentifiable simulations for max intake and especially the halfsaturation 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 predatorprey 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 Ivlevtype 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. Predatordependence 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 halfsaturation constant. We have done this in Supplement 2, and such model fitting works very well.
Our models used discrete time. Continuoustime 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 likelihoodfree 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 discretetime predatorprey 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, discretetime therefore seems a more convenient option in many respects, but continuoustime stochastic models also have clear advantages (e.g., no ordering of events and a good connection to most of the multispecies 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 nonreplacement 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 followup 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 predatorprey 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 (ANR10LABX45) and OG by the French National Research Agency (grant ANR16CE020007).
1 Stability analysis of the deterministic model and stochastic implications
The model, as specified in logscale by eqs. 7–8, can be analyzed by computing the Jacobian for the (nontrivial) 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
(A6)  
(A7) 
with , and
(A8)  
(A9) 
The Jacobian of the logtransformed 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 
2 Complementary results on the limit cycle parameter set
2.1 Fisher Information Matrix
2.2 Posterior correlations
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 densitydependence curve. These are similar to the main text results, except that the functional response can now be estimated without kill rate data for .
We present here plots for the prey densitydependence curve but similar results can be obtained for the predator densitydependence 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.
3 Complementary results on the reparameterized model
Here we report the results of the Bayesian estimation of the reparameterized model of eqs. 10–11 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).
Although are slightly less correlated than , the average functional response curve as a whole is not better estimated (Fig. 19).
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 densitydependence coefficient , which removes posterior correlations (see also main text for the likelihoodbased 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.
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 processbased model predictions. Journal of the Royal Society Interface 15:20180741.
 AugerMéthé et al. (2016) AugerMéthé, M., C. Field, C. M. Albertsen, A. E. Derocher, M. A. Lewis, I. D. Jonsen, and J. Mills Flemming. 2016. Statespace 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. Predatorprey feedback in a gyrfalconptarmigan 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 steadystate 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 overexploited 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 predatorprey food webs – confronting dynamic models with empirical data. Journal of Animal Ecology 88:196–210. doi:10.1111/13652656.12892.
 de Roos et al. (1991) de Roos, A., E. McCauley, and W. Wilson. 1991. Mobility versus densitylimited 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 predatorprey 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 predatorprey 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 markrecapturerecovery 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 predatorprey 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. Highamplitude 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 timeseries data. Population Ecology 43:229–243.
 Jost and Ellner (2000) Jost, C., and S. P. Ellner. 2000. Testing for predator dependence in predatorprey dynamics: a nonparametric 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 vectorborne 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 statespace 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 predatorprey 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 predatorprey 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. Nondeterministic modelling of foodweb 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. BarriosO’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. VucicPestic, and O. L. Petchey. 2012. Universal temperature and bodymass 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 capturerecapture 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 tradeoffs associated with predatorinduced 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 discretetime predatorprey 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 individualbased simulations. American Naturalist 151:116–134.
1 Estimation and identifiability of a RosenzweigMacArthur 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 discretetime 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  
Halfsaturation 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 
We added a parameter set leading to quasicycles (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 nearmonotonic decay to equilibrium in the absence of noise).
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).
The posteriorprior overlap curves are reported in the followings Figs. S2–S4. 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) 
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 
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 
2.1.3 Model with a Betadistributed 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 
2.2 Temporal variation in the halfsaturation constant
Temporal variation in
is added using a Gamma distribution for the halfsaturation constant, with mean
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 .
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 
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 
2.2.3 Model with a Gammadistributed halfsaturation , 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 
This model is quite well estimated, with excellent convergence; for the parameters chosen here the Gaussian FR and Gammadistributed halfsaturation 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.
Comments
There are no comments yet.