Predators' Functional Response: Statistical Inference, Experimental Design, and Biological Interpretation of the Handling Time

Citation: Papanikolaou NE, Kypraios T, Moffat H, Fantinou A, Perdikis DP and Drovandi C (2021) Predators’ Functional Response: Statistical Inference, Experimental Design, and Biological Interpretation of the Handling Time. Front. Ecol. Evol. 9:740848. doi: 10.3389/fevo.2021.740848 Predators’ Functional Response: Statistical Inference, Experimental Design, and Biological Interpretation of the Handling Time


INTRODUCTION
Functional responses describe the predator feeding rate with increasing prey density (Solomon, 1949) and are central to ecology, quantifying the energy transfer across trophic levels. Holling's approach (Holling, 1959a,b) has been the base upon which many of the critical aspects of predator-prey interactions can be detected (e.g., Abrams, 1980Abrams, , 1989. The most frequently observed and widely used functional responses in describing predator-prey relationships are that of the type II and III (Jeschke et al., 2002(Jeschke et al., , 2004, characterized by a curvilinear and a sigmoidal increase in feeding rate with prey abundance, respectively. Accurate and robust approaches for quantifying functional responses are critical to the investigation of predator-prey coexistence (e.g., Aldebert and Stouffer, 2018;Uszko et al., 2020;Barraquand and Gimenez, 2021;Coblentz and DeLong, 2021). Therefore, the estimation, as well as a mechanistic understanding of the parameters that determine predator feeding behavior is of importance. In this paper, we summarize advances related to experimental design, statistical analysis, and the mechanistic interpretation of the predation process that are central to the robust quantification of functional responses and hence should be adopted broadly.

STATISTICAL INFERENCE
In many cases, several functional response models are fitted to experimental data using methods such as non-linear least squares optimization (e.g., Juliano and Williams, 1987;Pervez and Omkar, 2005). However, such an approach provides no information about the uncertainty around the estimates and it may well be the case that there are other plausible parameter values that offer an equally good fit. Furthermore, a frequentist approach to uncertainty quantification (most often using maximum likelihood estimation-MLE) assesses the performance of a statistical estimation procedure on the basis of the expected long-run performance given a hypothetical series of datasets collected under identical conditions. Furthermore, the methods by which the uncertainty is quantified are typically constructed under parametric assumptions of the MLEs and with increasing accuracy observed as the size of a dataset increases. However, in many cases an experimentalist will have a dataset of a fixed size and may not always be sufficiently large for these asymptotic results to hold. Novak and Stouffer (2021) have recently highlighted and demonstrated using a large compilation of public datasets that there is systematic bias in the statistical comparison of functional response models and the estimation of their parameters which are rooted in a lack of sufficient replication, or in other words, small sample sizes. Furthermore, although it is important to account and quantify the uncertainty around the model parameters, one should not ignore or forget that this is done under the assumption of a particular functional response model (e.g., Holling's Type II). However, there is also uncertainty around the structure of functional response model (e.g., Type II vs. Type III and beyond). We discuss below how one can jointly perform parameter inference and model selection within a coherent probabilistic framework.
We advocate the use of a Bayesian framework. Such an approach treats all the unknown model parameters as random variables and first assigns them a prior distribution that represents our beliefs about the unknown parameters before any experimental data are collected. This prior information is subsequently updated in light of experimental data using Bayes theorem, leading to the posterior distribution that contains all information regarding the model parameters given both the experimental data and prior knowledge.
Bayesian statistical inference is not limited to parameter estimation. We are often interested in assessing a particular scientific hypothesis related to the functional response. For example, discriminating between type II and III functional responses is a very common question, as type II functional responses are known to destabilize predator-prey dynamics, in contrast to type III (Oaten and Murdoch, 1975). A Bayesian approach to model selection treats the model itself (as well as its parameters) as unknown and hence in addition to quantifying uncertainty about the model parameters, uncertainty about the model too is also taken into account. Given a series of plausible models, representing for instance, different forms of functional response, we specify a prior distribution for each model and prior distributions for the model parameters for each model and in light of experimental data, we can then obtain posterior model probabilities which represent our beliefs. That is, having obtained experimental data, what is the chance that a particular model out of the pool of models we are considering is the true one. This can be formalized using the notion of Bayes Factors, which is a summary of the evidence provided by the data in favor of one hypothesis represented by a statistical model as opposed to another. When there are more than two different model/hypotheses considered, e.g., different types of functional response such as type I, type II or type III, then it is best to consider the posterior model probabilities to identify to what extent each model is supported by the data.
A detailed description of the proposed methodology above and its application to several functional response data can be found in Bolker (2008) and Papanikolaou et al. (2016aPapanikolaou et al. ( ,b, 2020Papanikolaou et al. ( , 2021.

EXPERIMENTAL DESIGN
The previous section provides an overview of inferring model parameters and performing model selection after a dataset has been collected. However, it is often of interest to determine how best to conduct an experiment so that the resulting dataset that is collected will be most informative about the goals we are trying to achieve in an experiment, for example selecting the most appropriate model and estimating its parameters as precisely as possible. Such methods are referred to as optimal experimental design (OED, see for example Pukelsheim, 2006;Ryan et al., 2016). Formally, we define a utility function that encodes the goal of the study, and we plan an experiment that maximizes the expected utility function, i.e., the utility function averaged over all datasets that we might see for the planned experiment. By planning to collect the data in this statistically principled manner, we can reduce the amount of experimentation needed to achieve a particular statistical analysis goal, and hence reduce costs and required resources.
In OED there must be some variables that we can control in an experiment. In the context of functional response experiments, we are able to specify, for example, the initial prey density and time interval to use for an experiment. Therefore, OED tries to solve the problem of what initial prey densities and time intervals should be used such that the data collected will be as informative as possible. Once the initial prey densities and time intervals are determined, data can be collected, and the statistical inference methods in section 2 can be applied.
There are two main types of OED; static design and sequential design. In a static design, optimal input variables for all planned experiments are determined at the outset. In a sequential design, the optimal input variables for subsequent experiments are determined from the data obtained in previous experiments, by updating model structure and parameters. The sequential approach is generally more statistically efficient, since we use information from data collected sequentially to update our decisions, whereas the static design can only use information available prior to all experiments being conducted. One drawback of the sequential approach is that the overall data collection process may take longer since the experiments only take place one at a time or in small batches.
The first OED approach for functional response experiments is given in Zhang et al. (2018). The authors develop an approach to determine an optimal static design for estimating parameters of a functional response model. Here the utility function is based on the Fisher information matrix, which is one way to quantify how much information about model parameters we expect to gain from an experiment. This type of utility function assumes that maximum likelihood estimation will be performed on the data once the experiments at the optimal design are completed. Moffat et al. (2020) develop a Bayesian sequential design method for functional response experiments. This approach allows for multiple competing functional response models, not just a single model as in Zhang et al. (2018). Here, the utility function is based on a quantity called the total entropy, which computes the expected change in posterior distributions for both model probabilities and parameters. We prefer initial prey densities that lead to larger changes in posterior distributions, as this allows us to learn more about the preferred model and its parameters with less experimentation. Through extensive simulation studies, Moffat et al. (2020) show that the Bayesian sequential design approach leads to substantially more informative data compared to a Bayesian static design and a random design where initial prey densities are generated randomly. The improvement of the sequential over the static design results from updating information about the predator-prey system to make better decisions about future experiments.  (c = 0). Therefore, when a fixed handling time that is determined through separate short-term experiments thus excluding digestion effects is used, the disc equation no longer performs well for modeling consumption rates observed over a longer time period because the model does not take satiation into account.

BIOLOGICAL INTERPRETATION OF THE HANDLING TIME
Holling defined the handling time as the time a predator spends in pursuing, subduing and eating a prey item. This definition has been extended in many studies to include digestion time. However, the time a predator spends on pursuing a prey is followed by the time of subduing and, subsequently, consuming it. Digestion is a process that is likely to occur in parallel with these activities, in the sense that a predator can digest its prey while handling it. While the processes of searching for and handling prey are mutually exclusive, digestion is not (Jeschke et al., 2002). It should be noted that some predators can search for prey while handling others. In addition, pursuing a prey often does not result in capture, as digestion prevents successful hunting in many species. This variation in components does not produce Holling's disk equation formula (Abrams, 1982;Anholt et al., 1987).
Literature suggests that the majority of functional response experiments are conducted using a time interval that digestion effects are likely to be included in the estimated handling times, i.e., the predator daily foraging cycle (e.g., Cabral et al., 2009;Jalali et al., 2010;Fathipour et al., 2018;Islam et al., 2021). Figure 1 depicts the functional response curve of a predatory ladybird beetle on its aphid prey over a 24-h time interval (data from Papanikolaou et al., 2014). We fitted the data to the disc equation, as well as the steady-state satiation (SSS) equation model presented by Jeschke et al. (2002) that incorporates Holling's original definition of handling time while explicitly accounting for satiation and digestion time through a separate parameter. Both models showed almost identical fit to the data, realistically explaining the functional response of the predator. We also fitted the two models assuming handling was known, i.e., the handling time obtained from short-term functional response experiments, so that digestion effects are largely excluded. In this case, the disc equation does not appear to explain the functional response of the predator well. Therefore, in experimental procedures where digestion acts in parallel with prey handling (long term experiments), the interpretation of the "handling time" (T h ) estimated by the disc equation must be limited to the calculation of maximum attack rate (T/T h ), i.e., the maximum number of prey that can be consumed by a predator during the time interval (T) considered.

CONCLUSIONS
Over half a century after its conceptualization, predator functional responses remain a core feature in ecology. Recently, Novak and Stouffer (2021) revealed that bias in model comparison, as well as in parameter estimation are common in functional response studies, mainly attributed to the relatively small sample sizes used. OED has great potential for conducting more efficient functional response experiments. In addition, Bayesian inference enables the quantification of model and parameter uncertainty in a coherent, probabilistic manner, through the use of probability distributions.
Our understanding of predation could be further improved by elucidating the components of handling time. For example, Sentis et al. (2013) revealed that the processes of handling and digesting prey have different thermal responses. However, most functional response models incorporate handling time in a way that does not permit a biological interpretation, combining handling and digestion time (Jeschke et al., 2002). Although several modeling approaches has been presented (e.g., Mills, 1982;Abrams, 1990), using the SSS equation model would permit a mechanistic interpretation of these components of the predation process .
In conclusion, we advocate for the adoption of OED, a Bayesian framework, and the use of the SSS equation to efficiently and robustly infer predator functional responses moving forward. We anticipate that OED in conjunction with Bayesian inference will improve the predictive power of functional response experiments and reduce the logistical burden. Furthermore, as the comprehension of predator feeding behavior can be improved discriminating different predation processes such as handling and digestive prey, we call for the application of the SSS equation in functional response studies, which can lead to a better understanding of predator-prey interactions.

AUTHOR CONTRIBUTIONS
NP, AF, and DP provided ecological expertise. CD, TK, and HM provided statistical expertise. All authors contributed to the article and approved the submitted version.