Fitting functional responses: Direct parameter estimation by simulating differential equations

The feeding functional response is one of the most widespread mathematical frameworks in ecology, marine biology, freshwater biology, microbiology, and related scientific fields describing the resource‐dependent uptake of a consumer. Since the exact knowledge of its parameters is crucial to predict, for example, the efficiency of biocontrol agents, population dynamics, food web structure, and subsequently biodiversity, a trustworthy parameter estimation method is highly important for scientists using this framework. Classical approaches for estimating functional response parameters lack flexibility and often only provide approximations of the correct parameters. Here, we combined ordinary differential equation (ODE) models that were numerically solved using computer simulations with an iterative maximum likelihood fitting approach. We compared our method to classical approaches of fitting functional responses using data both with and without additional resource growth and mortality. We found that for classical functional response models, such as the frequently used type II and type III functional responses, the established fitting methods are reliable. However, by using more complex and flexible functional responses, our new method outperforms the traditional methods. Additionally, our method allows the incorporation of side effects such as resource growth and background mortality. Our method will enable researchers from different scientific fields who are measuring functional responses to calculate more accurate parameter estimates. These estimates will enable community ecologists to parameterize their models more precisely, thus allowing a deeper understanding of complex ecological systems, and will increase the quality of ecological prediction models.


| INTRODUC TI ON
Understanding the interactions between organisms is crucial to understand the patterns of stability and biodiversity of whole communities (McCann, 2000). Theory suggests that antagonistic interactions, such as feeding interactions, may have negative effects on stability (May, 1972) and, subsequently, on biodiversity (e.g., Rall, Guill, & Brose, 2008). Moreover, a good knowledge of interactions may allow ecologists to create realistic predictions of empirical communities ranging from microcosms to whole lakes (Boit, Martinez, Williams, & Gaedke, 2012;Fussmann, Schwarzmuller, Brose, Jousset, & Rall, 2014;Schneider, Scheu, & Brose, 2012). It is therefore of utmost importance to have reliable methods for quantifying the strength of interactions amongst organisms, especially consumer-resource pairs such as predators and their prey.

| Functional response models
The functional response (Holling, 1959b;Solomon, 1949) is one of the oldest and most commonly used mathematical frameworks to describe and estimate the feeding interaction between a consumer and a resource (see Jeschke, Kopp, & Tollrian, 2002 for an overview).
It describes the density-dependent feeding rate of a consumer on its resource (Holling, 1959b). Despite the plethora of possibilities, the most important models used today are the type II and type III functional responses (Jeschke, Kopp, & Tollrian, 2004;Skalski & Gilliam, 2001). The type II response is a hyperbolic saturating curve (Figure 1, blue line), where the per capita feeding rate F depends on the number of prey in the environment N: Here, a is the instantaneous rate of discovery (Holling, 1959b) that is commonly known and subsequently used here as the attack rate (note that depending on the scientific field, it may be known as capture rate (e.g., Kalinkat, Rall, Vucic-Pestic, & Brose, 2011), maximum clearance rate (e.g., Hansen, Bjornsen, & Hansen, 1997), maximum per capita interaction strength (e.g., McCann, Hastings, & Huxel, 1998), or others); and h is handling time (Holling, 1959a). See Jeschke et al. (2002Jeschke et al. ( , 2004) for a comprehensive introduction and discussion of the ecological meaning of these parameters. The attack rate a mainly controls the initial increase of feeding at low densities ( Figure 1a, grey dashed line), whereas the handling time h mainly controls the feeding at high densities, where the feeding curve is saturated (Figure 1a, grey dotted line).
The inverse of the handling time is often referred to as the maximum feeding rate F max , but it is also called the maximum ingestion rate (e.g., Hansen et al., 1997), the handling rate (e.g., Englund, Öhlund, Hein, & Diehl, 2011), or similar terms.
If the attack rate is not a constant but instead depends linearly on resource density (a = bN), the functional response becomes a sigmoid curve (e.g., Juliano, 2001, see Figure 1b, red line) that can be described by where b is the attack coefficient. Similar to the type II functional response, the attack rate a controls the feeding rate at low prey densities (Figure 1b, grey dashed line), and the handling time h controls the feeding rate at high prey densities (Figure 1b, grey dotted line).
See, for example, Juliano (2001) for more possibilities to model a type III functional response and Okuyama (2010) for resource density-dependent handling times. Real (1977Real ( , 1979 incorporated the possibility to gradually shift between the above types of functional responses by using the enzyme kinetics model of Barcroft and Hill (1910). In Real's adaptation, the attack rate depends on a power law of resource density (a = bN q ), leading to a functional response model that can be written as where q is an exponent that influences the shape of the functional response from a hyperbolic type II functional response (q = 0) to a strict type III functional response (q = 1) and beyond these bounds (Figure 1c). Note that the original formulation of this functional response is F(N) = F max N 1+q ∕(N 1+q half + N 1+q ), where F max is the already mentioned maximum feeding rate (h = 1/F max ) and N half is the half saturation density-the prey density at which the predator consumes half the amount of its maximum possible feeding rate (b = F max ∕N 1+q half ; Figure 1c, grey dotted-dashed line, please see the Supporting Information 1 (Section 1) for the equivalence of both parametrizations and Real (1977)). Moreover, the attack exponent q is often not used directly but in its transferable version, the Hill exponent H (H = q + 1).

| The problem of prey depletion-and its possible solutions
The above-mentioned functional response models (Equations 1-3) describe how the feeding rates depend on the number of prey in the environment (Figure 1a-c). This means that the independent variable (displayed on the x-axis) must be a constant prey density. In experiments, however, it is difficult to keep the prey (2) (2) the independent variable-the prey density-decreases over time. This phenomenon is commonly known as prey depletion (Rogers, 1972).
Describing the change in population density over time is classically done by using ordinary differential equations (ODEs); this approach is applied to model a wide range of ecological systems, from single populations to multitrophic communities (e.g., Delmas, Brose, Gravel, Stouffer, & Poisot, 2017;Rall et al., 2008;Rosenzweig & Mac Arthur, 1963;Verhulst, 1838  This equation describes the change in prey density dN over time dt, which is proportional to the feeding rate, that is, the generalized functional response (Equation 3), multiplied by predator density P.
ODE models must be integrated over time and density to yield the density at a desired point in time-here, the prey density at the end of the experiment N end (Figure 1d-f, black solid lines). This integration process can be done analytically (i.e., using mathematical rules to derive a new function N(t) having time as an independent variable and density as a dependent variable); alternatively, if the ODE is too complex to solve analytically, it can be solved numerically by using computer algorithms that approximate the analytical solution with a high accuracy. After integration, the number of prey eaten over time N e is calculated by simply subtracting the remaining prey density at the end of the simulation (as for experimental data) N end from the starting density N 0 (N e = N 0 − N end ). As mentioned above, this generalized Equation 4 becomes a type II functional response if the attack exponent is zero (q = 0, Figure 1d) or a type III functional response if the attack exponent becomes unity (q = 1, Figure 1e). In the following sections, we give an overview of the existing methods that use the functional response to predict the number of consumed prey N e .
The simplest approximation that could be made to describe the depletion of prey over the course of the experiment is to replace N by the known number of initial prey N 0 in the feeding rate F(N) and to multiply this feeding rate F(N 0 ) by the number of predators P (more generally, consumers) and the duration of the experiment T, yielding the number of prey eaten at the end of the experiment N e : Here, any of the above-mentioned functional response models (Equations 1-3) can be used for the feeding rate F. This method is a linearization of the nonlinear process that is precisely described by the ODE (Equation 4) and yields a linear decrease of live prey in the experiment over time (Figure 1d-f, dotted lines).
Instead of linearizing the ODE (Equation 4), Royama (1971) and Rogers (1972) presented an analytical solution of the ODE incorporating a type II functional response (q = 0), commonly known as With the unknown number eaten N e appearing on both sides of this implicit equation (RRPE-II), the problem arises that a simple nonlinear fitting algorithm is not applicable. The most frequently used solution to this problem is to use an iterative Newton's method (Juliano, 2001;Juliano & Williams, 1987). However, the state-of-the-art solution is to instead use the LambertW function, which allows an explicit solution of the implicit RRPE-II (Bolker, 2008): see also Figure 1d, dashed line. This approach has already been incorporated into a r-package for easy application (Pritchard, 2016;Pritchard, Paterson, Bovy, & Barrios-O'Neill, 2017). Juliano (2001)  Regardless of whether this equation is solved implicitly or explicitly, due to the approximation a = bN ≈ bN 0 , the predicted number of eaten prey N e is an approximation of the type III ODE and is not analytically correct (Figure 1e, dashed line).
As an alternative to the type III functional response, Hassell, Lawton, and Beddington (1977) presented a sigmoid-shaped feeding curve that accounts for prey depletion. Since the classical type III response is a special case of their wider class of curves (see While the above models can be used to fit either type II or type III functional responses, the following equation can also be used to fit the generalized functional response (Equation 3) to data with prey depletion. The power law dependency of the attack rate a = bN q ≈ bN q 0 was incorporated into the RRPE-II (e.g., Barrios-O'Neill et al., 2016;Kalinkat et al., 2013;Vucic-Pestic et al., 2010), yielding where the parameters are as above: the number eaten N e , attack constant b, attack exponent q, handling time h, total time T and number of predators P. We hereafter refer to this approach as the gener- As an alternative, Uszko et al. (2015Uszko et al. ( , 2017b) generalized a method developed by Frost (1972), in which (1) an approximation of the average resource density N is calculated (N 0 ≥N ≥ N end )-taking prey depletion into account-and (2) the functional response parameters are estimated using (for further information see Supporting Information 1 Section 4):

| Possible drawbacks of the different methods
Using the simplest approach, where the feeding rate is multiplied by the predator density and duration of the experiment (Equation 5), prey depletion during the experiment is neglected. Because this feeding rate is higher than the actual feeding rate, which is based on a reduced number of prey later in the experiment (N 0 > N causes Figure 1d-f, compare the dotted lines to the solid line created by the ODE representing the natural depletion of prey), the attack rate a will be underestimated, regardless of whether a type II, type III, or generalized functional response is used. If the feeding realized over the course of an experiment (i.e., the total number of consumed prey) approaches a saturation level, the estimate of the maximum feeding rate F max (or handling time h = 1 F max ) should only be marginally affected.
The RRPE-II is the exact analytical solution of an ODE that describes the depletion of prey over time with a type II functional response. The resulting depletion of prey over time is consequently identical when using the RRPE-II or the numerical simulation of an ODE (Figure 1d, dashed vs. solid line). This method therefore delivers reliable parameter estimates if the data can be described by a type II functional response.
As the RRPE-II is the exact solution for an ODE with a type II functional response, there can be negative consequences when it is used to fit nontype II functional responses, such as the type III and generalized functional responses. The attack rate a = bN q is, in this case, itself dependent on the prey number. Approximating N by using a larger N 0 creates a larger attack rate and, therefore, a greater decrease in prey than given in a simulated ODE (Figure 1e-f, dashed vs. solid line). A greater decrease in prey should lead to an underestimation of attack rates (i.e., attack coefficient b), but the handling time should not be affected if the number of consumed prey in the experiment reaches a saturation level. An increasing attack exponent q decreases the number of prey consumed at low prey numbers, as prey depletion is overpredicted by using Rogers' Random Predator equation with a type III or generalized functional response.
We assume that the attack exponent q is overestimated to counteract this effect.
We were not able to create strong expectations for the method by Frost (1972) and Uszko et al. (2015) since this two-step process cannot be visualized as descriptively as the previous methods.
However, as the original method is based on a type I functional response, we also assume deviations from the correct solution in all cases.

| The problem of nonconsumer-mediated prey growth and mortality
Another problem that often arises in experiments is that background mortality or growth of the prey could occur. For example, the assumption of constant prey numbers in the experiment could be violated not only because of prey removal by predators but also because of prey mortality or population growth, which is often the case if using microbial or algal prey (e.g., Fussmann, Rosenbaum, Brose, & Rall, 2017;Uszko et al., 2017b). Standard ways to deal with this are to exclude any data that are biased by mortality or growth, leading to a substantial loss of data, or to ignore the problem if it is not too pronounced, accepting that the resulting bias in the parameters is likely to be small (e.g., Uszko et al., 2017b;Vucic-Pestic et al., 2010). However, it is also possible to minimize bias without excluding data by fitting the parameters using simulations of ODEs (see, for example, Rall & Latz, 2016). Additional mortality can be accounted for by adding a mortality rate, m, to the ODE (Equation 4). This is a standard way to model mortality in population and food web modelling (e.g., Fussmann et al., 2014;Kalinkat et al., 2013;Rall et al., 2008;Rosenzweig & Mac Arthur, 1963;Williams & Martinez, 2004) for arbitrary functional response models F(N).
The growth of prey is often described by simple phenomenological growth models, such as logistic growth (Verhulst, 1838) and the Gompertz growth model (Gompertz, 1825;Paine et al., 2012). The logistic growth can be added to the ODE model by simply adding the growth term as the mortality above. For the logistic growth, the ODE becomes where K is the carrying capacity, that is, the number of prey a system can sustain, and r is the intrinsic growth rate that controls prey growth at low prey densities and determines prey mortality if the prey density is above the carrying capacity. The rest of the parameters are as above in the case of a simple functional response ODE.
In cases of both prey mortality and combined prey growth and mortality, we expect biased functional response parameters if the estimation is not corrected for these additional effects. In experiments where prey mortality occurs (Equation 14), this could lead to an overestimation of the attack rate or attack coefficient, an underestimation of the attack exponent (the lower the exponent is, the higher the feeding at low prey numbers is), and an underestimation of handling time (the lower handling time is, the higher the maximum consumed number of prey is) to counteract the additional effect of mortality. If the prey is growing and dying according to the above-mentioned growth model (Equation 15), there would likely be an underestimation of the attack rate (which describes feeding at low prey numbers) or attack coefficient and, potentially, an overestimation of the attack exponent to counteract the additional growth of prey at low numbers. At high prey numbers exceeding the carrying capacity, background prey mortality occurs, and there would likely be an underestimation of the handling time (which describes feeding at high prey numbers) to counteract this effect. Using extended ODE models and fitting them numerically to the data should overcome the problem. We further expect that using additional control data (no predators present) will improve the accuracy of the parameter estimation because the effects of natural growth (or death) and predation can be disentangled.
Here, we present a new framework for fitting functional response models to data. We combine numerical simulations of ODE models with an iterative maximum likelihood estimator using r we expect that our method will outperform traditional methods.
Moreover, fitting numerical simulations of ODE models allows us to cope with side effects such as prey growth and mortality that may otherwise bias or even prevent the estimation of functional response parameters.

| MATERIAL S AND ME THODS
We simulated 1,000 datasets based on ODE models (Equation 4) for each of the type II, type III, and generalized functional responses.
For numerical simulations, we used the deSolve package (Soetaert, Petzoldt, & Setzer, 2010  Through growth, the number of prey present at the end of the experiment N end can exceed the number of initial prey N 0 . Negative numbers of eaten or dead prey N e = N 0 − N end occur and make the binomial distribution unsuitable for this model. Instead, we use N end , which is always nonnegative. For each observation (N (i) 0 , N (i) end ) (i = 1, …, n), the prediction of the number of prey present N (i) end at the end of the experiment was calculated by a numerical simulation of the ODE. The observation N (i) end was modelled as being log-normally distributed with a location parameter log (N (i) end ) and scale parameter σ. This scale parameter is also automatically estimated in the iterative maximum likelihood routine (Crawley, 2012, chap. 7.3.3). The log-normal distribution considers that residuals generally have a larger variance for large prey abundances. We also chose it because microbial systems, where prey are measured as density, usually feature prey growth; however, for pure count data, a negative binomial or Poisson distribution can be used. See Supporting Information 2 Section 7.6 for a detailed description of this likelihood function.
For data without prey growth or mortality, we fitted all corresponding feeding models and compared the original simulation parameters to the estimates. For simulated data including background prey mortality, we tested the performance of fitting the ODE model neglecting mortality (Equation 4) and the ODE model including mortality (Equation 14). We also used the mortality model to fit simulated data, which used half of the 960 observations as control data; that is, no predators are present, and prey death occurs by natural mortality only. For data including prey growth and mortality, we also tested the performance of fitting the ODE model neglecting growth (Equation 4) and the ODE model including growth (Equation 15). We furthermore used the growth model to fit simulated data, which used half of the 960 observations as control data; that is, no death by predation occurs. We evaluated the performance of all models by comparing the log 10 -ratios of estimated and true parameters (Berlow, Navarrete, Briggs, Power, & Menge, 1999;Rall & Latz, 2016 ) with and without using control data.
In a last example, we used the dataset D8 (Fussmann, 2017), which features natural prey growth and death. The example also includes control data without predators present. Similar to the simulated data, we fitted the generalized ODE neglecting growth (Equation 4) and the model including a growth term (Equation 15) with and without using control data.

| Feeding models
The results of fitting type II, type III, and generalized functional response feeding models to the simulated data are depicted in Figure 2.

| Prey mortality and prey growth
We compared three new ODE approaches of dealing with natural prey mortality ( Figure 3). When fitting a mortality-free model (Equation 4) to data that include natural mortality, the attack coefficients and maximum feeding F max = 1 h are systematically overestimated to compensate for the missing mortality term. The estimates of the attack exponents are also biased, as expected. Using a model including mortality (Equation 14) leads to an unbiased prediction of all parameters. However, the error ranges for handling time h and the mortality coefficient m are still relatively large because the effects of predation and natural mortality are correlated. However, if control data are used, the accuracy of the predicted h and m is considerably improved because predation and mortality effects can be disentangled.

F I G U R E 3
Evaluation of the ODE model with an additional prey mortality term. Plots show the error distributions of fitted parameters log 10 fit simulated in 1,000 simulated datasets Similarly, we compared three new ODE approaches of fitting data including natural prey growth and death (Figure 4). Neglecting growth in the model (Equation 4), the handling times h are unexpectedly unbiased but feature a comparatively large error range. The attack exponents q are also unexpectedly slightly underestimated. In combination with the severely underestimated attack coefficients b, we anticipate that the attack rates a = bN q are generally underestimated because they must correct for the missing growth term in the model. When including this term (Equation 15), all functional response parameters, as well as the growth term's parameters, growth rate r and carrying capacity K, are unbiased. By fitting this growth model to data and using half of the number of observations as control data, the accuracy of unbiased estimates is drastically improved.

| Experimental data
When fitting our new ODE model (Equation 4) to experimental datasets, the predicted curves of eaten prey were almost identical to the studies' original methods ( Figure 5). Both methods performed similarly in fitting observed data. However, our focus lies on the parameters that produce these curves and on their discrepancies. Therefore, we do not provide model comparisons (e.g., AIC scores) but report these differences in the estimates (Table 1) and contrast them to our results obtained from the simulated data.
The handling time h estimates are generally similar when comparing the original and our new ODE method for all datasets. When considering a type III (dataset D1) or generalized response (datasets D2, D3) with similarly fitted attack exponents q, the attack coefficients b are underestimated by the RRPE and by Frost's method compared to the new ODE approach. This confirms our findings of the comparisons using simulated data and our assumptions about biased attack rates of the RRPE. By ignoring prey depletion and approximating a = bN by bN 0 (or a = bN q by bN q 0 in the generalized case), N 0 > N causes an underestimation of the attack coefficients b.
In datasets D4-D6, we observed both an under-and overestimation of attack coefficients b by the RRPE-gen. The attack exponents are also underestimated. We emphasize that for generalized functional responses, attack coefficients and attack exponents in a = bN q can only be interpreted jointly. For an alternative parametrization using the half saturation density to describe feeding at low densities independently from q, see the Supporting Information 1 Section 1 and Real (1977).
TA B L E 1 Functional response parameters for experimental data (D1-D6). n = number of observations, b = attack coefficient, h = handling time, q = attack exponent. Standard errors are shown in parentheses. p-values: · p < 0.1, *p < 0.05, **p < 0.01, ***p < 0.001. Units are: For the dataset D7 featuring natural prey mortality, we compared our generalized ODE model neglecting mortality (Equation 4) to the ODE model including a mortality rate (Equation 14) with and without using control data ( Figure 6, Table 2, see also Supporting Information 2 Section 3). Again, the different models produce similar curves when fitting data with predators present, but the parameter estimates are quite different. As expected, the maximum feeding rate F max = 1 h compensates for natural prey death and is therefore highly overestimated when fitting a model without a mortality term; that is, the handling time h is underestimated. When fitting the mortality model without using control data, a negative mortality rate m is predicted, which actually indicates natural prey growth instead of prey mortality. However, the estimate is not significant and features a high standard deviation because the effects of feeding and natural prey mortality (or growth) cannot be disentangled. Using the same mortality model and including additional control data, a positive mortality rate m is estimated, and we observe a reduced uncertainty (smaller standard error) in the estimates for m. Control data enable the precise separation of feeding and natural death.
Dataset D8 includes natural prey growth and death, and we compared our generalized ODE model ignoring background growth F I G U R E 6 Fits of our new ODE-based models to experimental datasets D7, including background prey mortality (top), and D8 including background prey growth and mortality (bottom); feeding experiments (left) and control experiments (right)  Fitting the growth model without control data produces heavily biased estimates. A large carrying capacity K and a large growth rate r lead to a decrease in the predicted dead prey for large initial prey densities. Using control data, the uncertainty in all parameters is significantly reduced (smaller standard errors) because the effects of feeding and natural growth (N < K) (or natural loss, N > K) can be separated precisely.

| D ISCUSS I ON
We compared different popular methods and our new ODE-based approach for fitting functional response models to feeding experiments using both simulated and experimental datasets. We estimated the parameters of the type II, type III, and generalized functional response models via maximum likelihood. Additionally, we investigated models including background prey mortality or growth.
Based on our results, we give the following recommendations:  (Hassell et al., 1977), as the fit of an analytical solution is much faster than fitting simulations to data (see Supporting Information 1 Section 6).
However, when including the attack exponent q as a free parameter, only the numerical simulation of ODEs guarantees unbiased estimates in our simulated experiments. This is of utmost importance for ecological research, as even a small shift of the attack exponent, q, from 0 (strict type II functional response) to approximately 0.2 (e.g., Williams & Martinez, 2004) increases the stability of trophic ecological networks (i.e., food webs) and thereby increases the species coexistence and biodiversity (Rall et al., 2008). Studies empirically investigating the attack exponent have increased in number over the last several years (Barrios-O'Neill et al., 2014b, 2016Kalinkat et al., 2013;Uszko et al., 2015Uszko et al., , 2017bVucic-Pestic et al., 2010 needed to understand the stability and diversity patterns in complex ecological communities more accurately.
In addition to these obvious and needed advancements, our ODE-based method is highly flexible. It is not limited to the three functional response models described here-indeed, any type of functional response describing prey-dependent feeding rates can be integrated. Using simulations of ODEs that are fitted to data further allows ecological modellers to include additional terms such as prey mortality or growth, as presented in this study. We showed, however, that control data are required to disentangle the effects of natural mortality or growth and predation.
In addition to this main article, we provided an in-depth description of our method as a manual, R-source files, and data. This will provide researchers who want to apply our method to their own functional response data with an easy entry into our methodology and the possibility to practice using the method with working examples before using it to analyse their own data. Moreover, the source files attached in the supplement can be adapted to include functional response models other than those used here (see Juliano, 2001 andJeschke et al., 2002 for a variety of possible models). In addition, specific mortality and growth models can be incorporated, for example, using the Gompertz growth model instead of the logistic growth model-a growth model that is often used in microbiology Paine et al., 2012).
We also show in the Supporting Information 1 Section 7 how data availability changes the outcome of parameter fitting. Generally, we found that even sparse data can be sufficient to estimate parameters correctly, but in particular, the estimation of handling time is highly uncertain if higher densities of prey are not measured in an experiment (Supporting Information 1 Figure 2).
While maximum likelihood estimation produces trustworthy point estimates, standard errors and confidence intervals rely on a quadratic approximation of the log-likelihood function (Bolker, 2008). For more accurate confidence intervals, resampling methods such as bootstrapping can be used (e.g., Pritchard et al., 2017).
Bayesian data analysis is more complex to implement and computationally more expensive, but it provides full probability distributions of the estimates. The deBInfer package (Boersch-Supan, Ryan, & Johnson, 2017) allows fitting ODEs to data in a Bayesian frame- work. An application of fitting ODE models to feeding experiments using Bayesian methods can be found in Fussmann et al. (2017).
In conclusion, we found that our new method to estimate functional response parameters from a feeding experiment with prey depletion-the decrease of prey number over time without replacement-has advantages over traditional approaches. To our knowledge, it is the only method that allows the precise parameter estimation of the generalized functional response, which is currently one of the most frequently used descriptions of feeding interactions in ecology. It moreover allows easy adaptation and the inclusion of unwanted effects such as prey growth and mortality, which is impossible with traditional approaches. We believe that our approach to estimate functional response parameters will enable researchers to better estimate parameter values using mortality-and growth-biased data in the future and thereby deepen the understanding of how interactions drive stability and biodiversity of complex ecological communities.