On the mechanistic understanding of predator feeding behavior using the functional response concept

. Functional response models describe the relationship between prey density and per capita prey consumption rate by a predator. Type II functional responses, in which density-dependent predation occurs via a decelerating feeding rate, seem to prevail in nature and are commonly described by Holling ’ s disk equation. In the derivation of the disk equation, Holling did not include digestion time. Although some authors have later extended the interpretation of handling time by also including digestion time, this violates the key assumption of the disk equation that the processes of searching for and handling prey are mutually exclusive. The steady-state satiation (SSS) equation is a functional response model that discrimi-nates between handling and digestion time. The application of the SSS equation is underutilized so far in the ecological literature, probably due to its complexity. In this study, we ﬁ rst tested the viability of the SSS equation. Second, we investigated the mechanistic basis of the SSS equation, comparing the model ’ s predictions with directly observed data. For this purpose, we used predator – prey systems of different taxa, that is, the ladybird beetle Hippodamia variegata preying on the aphid Aphis fabae , the lacewing Chrysoperla agilis preying on the aphid Myzus persicae , and the predatory mite Iphiseius degenerans preying on the thrips Frankliniella occidentalis . Our results show that the SSS equation is viable and can realistically describe type II functional response. In all predator – prey systems we tested, the model ﬁ tted the data reasonably well and provided realistic estimation of its parameters.


INTRODUCTION
Functional response models are central to predator-prey models, describing the relationship between prey density and per capita prey consumption rate by a predator (Solomon 1949), and are typically classified according to their shape (Holling 1959a). If the consumption rate of a predator rises linearly with prey density (as is assumed in many population dynamics models), we speak of a linear functional response. Such a response can theoretically be seen if both handling and digestion time are negligible. In practice, however, linear functional responses are rare, as digestion time is hardly ever negligible (Jeschke et al. 2004). If consumption rate first rises linearly with prey density, but remains constant once a threshold prey density has been reached, we speak of a type I functional response. In type II functional responses, density dependence occurs via a continuously decelerating feeding rate, and a type III functional response is sigmoid. The plateaus of type II and III responses are typically determined by handling or digestion time. Type II functional responses seem to prevail in nature (Jeschke et al. 2004). Holling's (1959b) seminal modeling approach for type II functional responses, commonly known as the disk equation, has been the base upon which much of modern functional response theory has been developed, becoming the null functional response model (Jeschke et al. 2002, Englund et al. 2011. While examining the foraging cycle of individuals, Holling used two parameters in order to develop an improved explanation of their feeding behavior: the attack rate, that is, the predation ability at low prey densities, and the handling time, that is, the time a predator spends pursuing, subduing, and eating a prey item. The disk equation is of the form: where N denotes the prey density; P the predator density; a the predator attack rate, and T h the handling time. Thereafter, several modeling approaches for type II functional responses have been pursued with respect to the attack rate (see Jeschke et al. 2002 for a review). The discrimination of these models is typically based upon their assumption on predation limitation. The limitation regarding the number of prey attacked at higher densities is often attributed to satiation or handling time constraints. However, although digestion and handling prey are discrete phenomena, digestion must be considered as a background process, acting in parallel to handling (attacking and eating) prey, which influences the predator's hunger level and consequently its probability of searching for prey (Jeschke et al. 2002). In the derivation of the disk equation, Holling did not include digestion time. Although some authors have later extended the interpretation of handling time by also including digestion time, this violates the key assumption of the disk equation that the processes searching for and handling prey are mutually exclusive, as a predator can search for prey while digesting its last meal (i.e., the processes searching for and digesting prey are not mutually exclusive). This leads to further examination of functional response models and an argument about their mechanistic basis.
The mechanistic basis of the disk equation has been challenged by Jeschke et al. (2002), as they presented a detailed examination of the fundamental process of predator feeding behavior, discriminating handling and digestion time. In this task, the authors developed a functional response model that realistically incorporates attacking and eating prey, the steady-state satiation (SSS) equation, which for digestion-limited predators, is of the form: where b denotes the corrected handling time and c the corrected digestion time of a single predator, preying on a randomly distributed single type of prey. In this model, the corrected handling time is equal to t att /e + t eat , where t att denotes the attack time per prey, e the proportion of successful attacks (i.e., attack efficiency), and t eat the eating time per prey. The corrected digestion time is equal to st dig , where s denotes the satiation per prey item, defined as the reciprocal of gut capacity, and t dig the digestion time per prey item, defined as the food transit time in the gut (Sentis et al. 2013). The key assumption of this model is a separation between handling time and digestion time, based on the observation that only handling but not digesting a prey item is incompatible with searching for additional prey. Furthermore, a predator's searching effort is proportional to its hunger level (i.e., the emptiness of a relevant part of the gut), which, by a separation-of-timescale argument, is assumed to take on a steady-state value dependent on prey density (Jeschke et al. 2002). The application of the SSS equation is underutilized so far in the ecological literature, probably due to its complexity relative to the simpler disk equation. Thus, there is only little available information on its practical viability and mechanistic basis (e.g., Tollrian 2005, Jeschke andHohberg 2008). Given that functional response models form the basis of population dynamics and foraging theory (Jeschke and Tollrian 2005), in this study we first test the viability of the SSS equation adopting a Bayesian framework, in order to quantify its estimated parameters in a coherent probabilistic manner. Second, we investigate the mechanistic basis of the SSS equation, comparing the model's predictions with directly observed data, as the parameters in a mechanistic model must have biological interpretation in order to be able to be measured independently of the functional response data. For this purpose, we used predator-prey systems of different taxa, that is, the ladybird beetle Hippodamia variegata Goeze (Coleoptera: Coccinellidae) preying on the black bean aphid Aphis fabae Scopoli (Hemiptera: Aphididae); the lacewing Chrysoperla agilis Henry et al. (Neuroptera: Chrysopidae) preying on the green peach aphid Myzus persicae Sulzer (Hemiptera: Aphididae); and the predatory mite Iphiseius degenerans Berlese (Acari: Phytoseiidae) preying on the western flower thrips Frankliniella occidentalis Pergande (Thysanoptera: Thripidae).
Chrysoperla agilis-Myzus persicae. -A M. persicae stock colony was established with individuals collected from tobacco fields in the area of Komotini (Northern Greece) and maintained on pepper (Capsicum annuum L.) plants at 25°AE 1°C, 60-70% RH, and 16:8 LD at the Department of Agricultural Development (Democritus University of Thrace). A colony of the lacewing C. agilis was established with individuals sampled in Crete (Southern Greece) and was maintained at 25°AE 1°C, 60-70% RH, and 16:8 LD at Democritus University of Thrace. Adults were reared in cylindrical cages on a liquid diet consisting of water, honey, sugar, and yeast hydrolysate, and juveniles were reared on Ephestia kuehniella eggs provided ad libitum in the rearing cages, as described by Pappas et al. (2007). Lacewing species identification was performed by Prof. Charles Henry (University of Connecticut).
Iphiseius degenerans-Frankliniella occidentalis.-Thrips (F. occidentalis) were originally collected from Chrysanthemum sp. plants in Orestiada (Northeastern Greece). The colony was maintained at the Department of Agricultural Development at Democritus University of Thrace on cucumber (Cucumis sativus L., cv Ginga F1, Geostore SA) leaves placed on wet cotton wool in Petri dishes in a climate room at 25°AE 2°C, 16:8 LD, and 60-70% RH as described in Pappas et al. (2018). An I. degenerans colony was established with individuals collected from citrus orchards in Arta (western Greece). The mites were maintained on detached bean leaves (Phaseolus vulgaris L.) kept in contact with water-saturated cotton wool in plastic cups at 25°C, 60-70% RH, and 16:8 LD at Democritus University of Thrace and fed with Typha latifolia (L.) pollen as described by D€ oker et al. (2015).

Functional response experiments
Long-and short-term functional response experiments were conducted, using the experimental approach of Papanikolaou et al. (2014). Concerning the experimental design, we used the daily foraging cycle of the predators in order to conduct long-term functional response experiments. We also performed short-time experiments that largely excluded digestion effects.
Hippodamia variegata-Aphis fabae.-The experiments were carried out under laboratory conditions at 22°AE 1°C, 65 AE 2% RH, and a 16:8 LD photoperiod. We used 20-to 30-d-old female adults of H. variegata. All individuals were starved for 24 h and placed individually in plastic containers (9 cm diameter 9 2 cm height) with various densities of nymphs of A. fabae (3to 4-d-old) on leaves and branches of V. faba plant hosts. After 2 (short-term experiment) or 24 h (long-term experiment), the predators were removed and the number of aphids remaining counted. Prey densities were as follows: 1, 2, 4, 8, 16, 24, and 32 aphids for the short-term experiment and 2, 4, 8, 12, 16, 32, 64, and 128 aphids for the long-term experiment. There were 8-12 replicates of each prey density.
Chrysoperla agilis-Myzus persicae.-Experiments were carried out under laboratory conditions at 25°AE 1°C, 60 AE 5% RH, and 16:8 LD. We used young second-instar lacewing larvae that had been reared on M. persicae nymphs provided daily in Petri dishes (5.5 cm in diameter). Before the experiments, the larvae were left to starve without food for 12 h. Afterward, each larva was transferred on a pepper leaf disk (3 cm in diameter) placed with its abaxial surface on wet cotton wool in a Petri dish with various densities of aphid (M. persicae) second-and third-instar nymphs. After 2 (short-term experiment) or 24 h (long-term experiment), the predators were removed and the number of live aphids counted. Prey densities were as follows: 2, 4, 16, 32, 64, 80, and 128 aphids for both the short-and long-term experiments. There were 8-12 replicates of each prey density.
Iphiseius degenerans-Frankliniella occidentalis. -Experiments were carried out under laboratory conditions at 25°AE 1°C, 60 AE 5% RH, and 16:8 LD. We used 2-to 3-d-old adult females of the predatory mite I. degenerans. All individuals were starved for 24 h. Afterward, each adult mite was transferred on a bean leaf disk (3 cm in diameter) placed with its adaxial surface on wet cotton wool in a Petri dish with various densities of first-instar thrips (F. occidentalis) larvae. After 3 (short-term experiment) or 24 h (long-term experiment), the predators were removed and the number of live thrips counted. The length of the short-term experiments differed due to the significantly lower predation efficiency of I. degenerans compared to H. variegata and C. agilis. Prey densities were as follows: 1, 2, 3, 4, and 5 thrips larvae for the short-term experiment and 2, 4, 8, 12, and 16 thrips larvae for the long-term experiment. There were 8 replicates of each prey density.

Direct observations on handling time
In order to directly measure predator handling time, we used the approach of Sentis et al. (2013). In this task, we defined the handling time as the time interval from the beginning of an attack of a predator on its prey to the moment when the predator finished eating, that is, end of chewing/ sucking and resumption of searching behavior (walking and head swinging). The experiments were conducted in the same experimental arenas as the functional response experiments. There were 10 replicates for each predator-prey system.

Statistical analysis
Model.-In this paper, we adopt the same modeling framework as the one presented in Papanikalou et al. (2016a). Denote by N e (t) the number of prey eaten by time t. Since a prey item is either dead or alive by time t (which often denotes the end of the experiment), we assume that N e (t) follows a Binomial distribution with parameters N 0 and p(t), where N 0 is the initial prey population size and p(t) the probability that a prey item has been eaten by time: where N(t) is given by the solution of the ordinary differential equation (ODE) that corresponds to the chosen functional response model (e.g., Holling's disk or the SSS equation) evaluated at time t. Note that in all but the simplest functional response models, the ODEs have to be solved numerically. Bayesian inference. -We adopt a Bayesian approach to quantify the uncertainty around the model's parameters in a coherent, probabilistic manner (e.g., Bolker 2008).
Fitting the Holling's disk equation model to the data of short-term experiments.-We first assign uninformative uniform priors to both the attack rate (a) and handling time (T h ) in the Holling's disk equation, namely a~U[0,4] and T h~U [0,2]. We then derive the likelihood of the data observed in the short-term experiments at T = 2 h or T = 3 h, respectively; the likelihood of the observed data X given the parameters h = (a, T h ) is of the form (Papanikolaou et al. 2016a): In this function, n j consists of the initial prey densities for j = 1, . . ., 5 or 6 and x j refers to the number of prey eaten by the predator. When the likelihood is combined with the prior densities for a and T h via the Bayes theorem, it gives rise to the posterior distribution of the parameters, p (ɑ, T h |Χ). We sample from the latter using a random-walk Metropolis algorithm (see Gamerman and Lopes 2006, for example).
Fitting the SSS equation to the data of long-term experiments.-We approximate the posterior distribution of the handling time T h (estimated from the short-term experiments for the disk equation) by a Gamma distribution (see Fig. 1), which we then use as a prior distribution for the parameter b of the SSS equation. For the parameters, a and c, we use independent exponential prior distributions with rate 0.1 (mean 10). Note that these priors are essentially flat (due to the high mean), but come from the same family of distributions as the prior for b (since the exponential is a special case of a Gamma distribution). We then fit the SSS equation to the experimental data up to 24 h, using an independence Metropolis-Hastings algorithm, where the proposal distribution is a Gaussian approximation to the posterior density with mean the mode of the posterior density and variance-covariance matrix the inverse of the negative Hessian matrix.
All statistical analyses were conducted in R (R Development Core Team 2015). The corresponding code to produce the results in the paper is available at https://github.com/kypraios/functional_ response.

RESULTS
At 24-h exposure time, the handling time estimated by the disk equation for H. variegata, C. agilis, and I. degenerans was 0.479, 0.313, and 2.885 h, respectively, decreasing to 0.133, 0.056, and 1.33 h, respectively, at the short-time exposure (Table 1). This indicates that the functional response for all three predators is limited by digestion. Similarly, estimated attack rates were higher during the short-term experiments for H. variegata, (0.963 h À1 ), C. agilis (0.555 h À1 ), and I. degenerans (1.065 h À1 ) than in the long-term experiments (0.180, 0.172, and 0.232 h À1 , respectively).
The SSS equation, as well as the disk equation, fitted the observed long-term data reasonably well (in fact, the fitted curves from the different models are nearly identical as shown in Fig. 2). The estimated parameters for each of the tested predators, as well as the corresponding credible intervals, are presented in Table 2. Based on the 95% credible intervals, there were no significant differences between the estimated attack rates of the disk equation and the SSS equation. In addition, there were no significant differences between the direct observations on handling time and handling times estimated by the SSS equation for all three predators.
The digestion time estimated by the SSS equation for H. variegata, C. agilis, and I. degenerans was 0.460, 0.303, and 2.525 h, respectively. For each predator, these values did not differ significantly from the estimated handling times by the disk equation, based on the 95% credible intervals (Tables 1 and 2).
We also fit the SSS equation and the disk equation to the long-term data assuming that the handling time is equal to the direct observations, which corresponds to 0.170 h for H. variegata, 0.047 h for C. agilis, and 0.620 h for I. degenerans. In this case, the disk equation does not explain the functional response of the predators, in contrast to the SSS equation (Fig. 2).

DISCUSSION
Knowledge of the trophic interactions between living organisms may contribute to a better understanding of the evolution of food webs. Given the importance of the functional response concept as a crucial component in community ecology, linking food webs, and determining predator-prey dynamics, the use of a mechanistic functional response model that accurately describes predation is imperative. The SSS equation goes in this direction, as its formal derivation was based on fundamental knowledge of the interactions between process variables was used and also the authors provided a detailed examination of the predation process. Our data on the functional response of three predatory arthropods confirm that the SSS equation is viable and can realistically describe type II functional responses. In all three predator-prey systems, the model fitted the data reasonably well and provided realistic estimation of its parameters. The estimated handling times are in accordance with direct observations on the time that H. variegata, C. agilis, and I. degenerans spent on pursuing, subduing, and eating A. fabae, M. persicae, and F. occidentalis individuals, respectively. Also, the plateau of the functional response curve, which determines the maximum number of prey consumed by a predator in a time interval (i.e., the maximal intake rate), is well defined by the reciprocal of the digestion time, as digestion limits the predation efficiency of these predators at long timescales.
It is worth noting that for I. degenerans, the directly observed handling time is about 50% of the corrected handling time estimated from the SSS equation, although credible intervals overlap. Considering the low prey consumption of this predator, it seems that digestion breaks may affect its predation ability immediately after the first prey item consumed, in contrast to H. variegata and C. agilis. This fact is also reflected in the short-term handling time estimated from the disk equation, which we used as a prior distribution for the corrected handling time. While the SSS equation is based on the assumption that predator satiation (or hunger level) is at a steady state, it seems that it is flexible enough to adapt a violation of its assumptions, as in our experiments predators are initially starved, and a potential steady-state satiation sets in only later. It is therefore suggested that short-term experiments must be designed over a time period that largely excludes digestion effects, especially in cases where predators consume few prey items. Holling's (1959b) basic assumption when deriving the disk equation was that predator time is divided into two separate periods, that is, "searching for" and "pursuing, subduing, and eating" prey. Hungry predatory arthropods are likely to devour the first few prey items they encounter (e.g., Hodek 1996, Maselou et al. 2015. In such circumstances, digestion effects may be negligible so that the handling time limits the predation efficiency and the functional response can be explained by the disk equation. Thereafter, digestion-limited predators become less efficient due to the on-setting digestion process, which determines their predation efficiency. As a consequence, if long-term functional responses are to be measured, the estimated handling time of the disk equation has no biological meaning and can be only used to determine the plateau of the functional response curves and therefore the maximal intake rate of the digestion-limited predators. Our results are in accordance with previous studies pointing out that real handling time and estimated handling time from the disk equation frequently have little in common (Caldow and Furness 2001, Jeschke and Tollrian 2005, Jeschke and Hohberg 2008. The discrimination between the handling and digestion time using the SSS equation could also lead to a more reliable examination of the predation process, as we show in this study. The fit of the disk equation and the SSS equation to the long-term functional response data was almost identical. The disk equation is the most frequently used functional response model in the scientific literature in either its basic or integrated form (for non-constant prey densities; e.g., Rogers 1972, Piersma et al. 1995, Caldow and Furness 2001, Sentis et al. 2012. This is probably due to its simplicity and tractability compared to other models, describing predator feeding rate in a relatively simple but realistic way. Jeschke et al. (2002) extended the basic assumptions of the disk equation when developing the SSS equation. Our data support the fact that these assumptions are reliable, leading to a biological interpretation of the estimated parameters, so that the SSS equation is a viable alternative to the disk equation.
The application of the SSS equation model is underutilized so far in the ecological literature, partly due to its complexity. Fitting such models to experiments can be challenging due to issues of practical identifiability, which is when the experimental data do not provide enough information about the parameters of interest (Papanikolaou et al. 2016b). In this paper, we overcame such issues by taking advantage of the Bayesian framework and using informative priors for the parameter that is common across models.
Ecologists often utilize functional responses to illustrate coevolutionary associations and improve functional predictive power in predator-prey dynamics (e.g., Gutierrez et al. 2008, Garay et al. 2015, Sentis and Boukal 2018. Dividing and estimating the foraging activities of predators could lead to a more efficient description and prediction of their functionality. In this study, we show that the discrimination between handling and digestion time is possible. However, several factors may affect predator feeding rate, such as alternative prey (Schenk and Bacher 2002), interference competition (Papanikolaou et al. 2016b), predator confusion (Jeschke and Tollrian 2005), and temperature (Sentis et al. 2012). For example, a further examination of how functional response parameters vary with temperature (Englund et al. 2011) could lead to a more efficient description of predator feeding rate in several biotic and abiotic conditions.
In conclusion, we provide evidence on the viability and the mechanistic basis of the SSS equation. This is critically important for utilizing functional response outcomes to accurately infer the components of predation at a population and community level, and also use the functional response as a key component to dynamically couple prey and predator populations in models of predator-prey systems.