# On the use of functional responses to quantify emergent multiple predator effects

Arnaud Sentis, David Boukal

<figure><img src="image_1.png" /></figure>

## To cite this version:

Arnaud Sentis, David Boukal. On the use of functional responses to quantify emergent multiple predator effects. Scientific Reports, 2018, 8 (1), 10.1038/s41598-018-30244-9. hal-03346852

HAL Id: hal-03346852

https://hal.inrae.fr/hal-03346852v1

Submitted on 6 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
SCIENTIFIC REPORTS

<figure><img src="image_3.png" /></figure>

OPEN

# On the use of functional responses to quantify emergent multiple predator effects

Arnaud Sentis<fnref n="1" />,²,³,⁴ & David S. Boukal²,³

Received: 20 October 2017

Accepted: 20 June 2018

Published online: 06 August 2018

Non-independent interactions among predators can have important consequences for the structure and dynamics of ecological communities by enhancing or reducing prey mortality rate through, e.g., predator facilitation or interference. The multiplicative risk model, traditionally used to detect these emergent multiple predator effects (MPEs), is biased because it assumes linear functional response (FR) and no prey depletion. To rectify these biases, two approaches based on FR modelling have recently been proposed: the direct FR approach and the population-dynamic approach. Here we compare the strengths, limitations and predictions of the three approaches using simulated data sets. We found that the predictions of the direct FR and the multiplicative risk models are very similar and underestimate predation rates when prey density is high or prey depletion is substantial. As a consequence, these two approaches often fail in detecting risk reduction. Finally, parameters estimated with the direct FR approach lack mechanistic interpretation, which limits the understanding of the mechanisms driving multiple predator interactions and potential extension of this approach to more complex food webs. We thus strongly recommend using the population-dynamic approach because it is robust, precise, and provides a scalable mechanistic framework to detect and quantify MPEs.

Determining the factors and mechanisms influencing the distribution and strengths of species interactions is crucial to understand and predict the dynamics and stability of ecological communities<fnref n="1" />⁻⁵. However, quantifying species interaction strengths (i.e., the effect of one species on the abundance of a second one, such as the effect of a predator on its prey) is challenging because species are embedded in complex communities and interaction strength between two species often depends on direct and indirect interactions with other species in the community. Observed community dynamics thus often differ from simple predictions based on pairwise interactions⁶⁻⁸. These differences emerge due to multiple-species interactions, often referred to (across various contexts) as higher-order interactions, trait-mediated indirect interactions, non-trophic interactions, and non-consumptive predator effects⁷⁻¹⁰.

Most previous studies of emergent multiple-species interactions have focused on multiple predator species with a shared prey⁶,⁷,⁹,¹¹. They revealed that interactions among different predators (e.g., competition or intraguild predation), as well as phenotypic and behavioral responses of prey to predation risk, often lead to emergent multiple predator effects (MPEs), where prey consumption rates by multiple predators are inconsistent with the assumption of independent effects of each predator on prey survival⁹,¹¹. MPEs are widespread in nature and can have profound effects on prey survival and transfer of energy and nutrients across trophic levels in food webs¹²⁻¹⁴. For instance, previous studies in agroecosystems showed that greater predator diversity can increase herbivore suppression and increase plant growth when predators facilitate each other¹⁵. The combined predation rate of multiple predators can sometime double the sum of their individual predation rates, showing that synergistic interaction among predators can be substantial and can improve biological control of agricultural pests¹⁶. Better understanding of these higher-order interactions is thus important for a full understanding of food web dynamics and ecosystem functioning⁶,¹⁰,¹⁷.

<footnote n="1">¹Unité Mixte de Recherche 5174 "Evolution et Diversité Biologique," Centre National de la Recherche Scientifique - Université de Toulouse III - Institut de Recherche pour le Développement, 31062, Toulouse, France. ²University of South Bohemia, Faculty of Science, Department of Ecosystem Biology, Branišovská 31a, 37005, České Budějovice, Czech Republic. ³Czech Academy of Sciences, Biology Centre, Institute of Entomology, Branišovská 31, 37005, České Budějovice, Czech Republic. ⁴Ecological Networks and Global Change Group, Experimental and Theoretical Ecology Station, UMR5321, CNRS - University Paul Sabatier, Moulis, France. Correspondence and requests for materials should be addressed to A.S. (email: Arnaud.sentis@sete.cnrs.fr)</footnote>
The most common approach used to quantify MPEs is to (1) experimentally estimate prey survival in the presence of each predator species in monoculture, (2) use this information to predict prey survival in the presence of multiple predator species using a mathematical model (hereafter; null model) assuming no emergent effects (i.e., predator species effects combine independently), and (3) compare null model predictions to empirical observations in which prey are exposed to a combination of the predator species. When predictions and observations do not differ, multiple predator species combine independently indicating the absence of emergent multiple predator effects. However, if predictions overestimate or underestimate observations, multiple predator species combine antagonistically (leading to prey risk reduction, where fewer prey are killed than predicted) or synergistically (leading to prey risk enhancement, where combined predators kill more prey than predicted by their individual effects). As the sign and magnitude of the differences among model predictions and empirical observations determine the nature and strength of MPEs, the accuracy and robustness of the null model predictions are crucial to detect and quantify MPEs.

Various null models have been proposed to predict the expected independent effects of multiple predators on prey survival<sup>9,10,18–20</sup>. A consensus emerged in the literature<sup>19,21,22</sup> that the multiplicative risk model<sup>18</sup> is the most accurate, and it has been widely used in multiple predator studies. The multiplicative risk model states that, if two predators have independent effects, the expected number of prey eaten by both predators can be predicted using the following equation:

$$N_{AB} = N_0(P_A + P_B - P_A P_B) \quad (1)$$

where $N_{AB}$ is the predicted number of prey eaten by predators A and B when foraging together, $N_0$ is the initial number of prey, and $P_i$ ($i=A, B$) is the probability of being eaten by predator $i$. However, recent studies argued against the use of the multiplicative risk model as it assumes constant prey mortality rate over the course of the experiment, which is only valid when predators have linear functional responses, when prey are continuously replenished (i.e. no prey depletion), or experimental durations are sufficiently short that depletion is inconsequential<sup>23–25</sup>. These assumptions are routinely violated in empirical studies of MPEs: functional responses of predators are nearly always a saturating function of prey density<sup>26</sup> and prey depletion is often non-negligible. For example, McCoy et al.<sup>13</sup> examined 100 multiple predator studies reviewed by Vance-Chalcraft et al.<sup>11</sup> and found that, on average, prey were depleted by 70% over the course of these experiments. A simulation study by McCoy et al.<sup>13</sup> demonstrated that these violations can severely affect the interpretation of the results of MPE studies. In particular, the above limitations of the multiplicative risk model can bias results at high prey densities or when prey depletion is substantial<sup>13,24</sup>.

To overcome the biases of the multiplicative risk model, recent MPE studies have stressed the importance of relationships based on mechanistic models and experimental designs that incorporate a gradient of the key ecological variable, such as prey density, rather than simple factorial designs that do not account for time dependence and nonlinear interaction strengths<sup>12–14</sup>. In particular, the shape of the predator functional response (i.e., the relationship between the number of available prey and the number of prey eaten by a predator) and predator-prey dynamics have been identified as crucial to (1) adequately assess the presence and strengths of MPEs, (2) identify the mechanisms driving changes in species interactions, and (3) extrapolate the results to larger temporal and spatial scales<sup>12,13</sup>.

Two approaches following these ideas have recently been proposed to investigate emergent MPEs. Both are based on functional response modelling and aim to circumvent the limitations of the multiplicative risk model by accounting for non-linear predation rates and prey depletion. In the approach proposed by Wasserman et al.<sup>20</sup> (hereafter; direct FR approach), the null model (to which predictions are compared with experimental observations involving multiple predators to estimate MPEs) combines a functional response model with the multiplicative risk model (eqn 1). In the approach proposed by McCoy et al.<sup>13</sup> (hereafter; population-dynamic approach), the null model combines a functional response model with a population dynamic model. The population-dynamic approach accounts for prey depletion and non-linear feeding rates (see Material and Methods for more details) and thus produces unbiased predictions about independent effects of multiple predators on prey survival<sup>13</sup>. It can thus be considered as the reference null model to which other approaches can be compared.

The direct FR approach is also appealing, but it remains unclear if and when its predictions differ from the multiplicative risk model and whether the underlying assumptions of the direct FR approach are robust and do not suffer from the same biases as the multiplicative risk model. We argue that comparing the different models/ approaches used to predict multiple predator independent effects on prey survival is an important issue for (1) correctly assigning the nature (risk reduction versus enhancement) and magnitude of emergent MPEs and (2) identifying modelling approaches that are best suited to quantify MPEs. The objective of this study was thus to compare the strengths and limitations of the multiplicative risk model, the direct FR model and the population dynamic model, and illustrate their use and compare their predictions using simulated data sets. We believe that this comparison will help advance the study of MPEs and improve the realism of predator-prey and food-web models.

## Results

When we compared the data from the simulated trials including two independent predators and the model predictions of each modelling approach for selected parameter values (Fig. 1), we found that the three modelling approaches captured well the non-linear relationship between prey density and the number of prey eaten by two predators. The population-dynamic model predictions also fitted well the simulated data, but the two other approaches (multiplicative risk model and direct FR model) tended to underestimate predation rates, especially for high attack rate and low handling time values. Moreover, predictions of the direct FR approach and the multiplicative risk model were nearly identical across the explored parameter space (Fig. 1).
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Simulated number of prey eaten over 24 hours by two conspecific predators under the assumption of independent predator effects on prey survival (black dots) for different values of attack rate and handling time. Lines represent predictions from the multiplicative risk model (solid black line), the direct FR model (dotted blue line) and the population-dynamic model (dashed red line). Range of values used for attack rate and handling time were chosen to reflect values found in previous studies conducting both a functional response experiment and multiple predator trials<sup>10,20</sup>.</figcaption></figure>

Further comparisons of model predictions to the simulated numbers of prey eaten across different values of the functional response parameters and experimental durations confirmed that the predictions of the direct FR approach and the multiplicative risk model tended to systematically underestimate predation rates by both conspecific (Figs 2 and 3) and heterospecific predators (Figs 4 and 5). Magnitude of this underestimation increased with experimental duration and predator attack rate *a*, decreased with handling time *h*, and either increased with or had a unimodal dependence on prey density in the studied range (Figs 2 and 3). As a result, the largest underestimations of predation rate were observed for 24-hour predation trials involving predators with the lowest values of handling time and high prey densities (Fig. 2). The lowest underestimations of predation rate were observed for 1-hour predation trials in which differences between model predictions and simulated data were negligible irrespective of the functional response parameter values and prey density (Fig. S1). Interestingly, underestimations tended to be less variable and less extreme for heterospecific than for conspecific predator assemblages. For instance, when handling time varied and attack rate remained constant we found that, at the highest prey density and 24-hour experimental duration, predation rate underestimations by the multiplicative risk model ranged from -3 to -71 for conspecific assemblages (Fig. 2) whereas they ranged from -10 to -51 prey for heterospecific assemblages (Fig. 4).

Unlike the multiplicative risk model and the direct FR model, the population-dynamic model predicted well the simulated numbers of prey eaten by two conspecific or heterospecific predators (Figs 2–5). Moreover, we found no systematic directional bias (over- or underestimation) for the predictions of the population-dynamic model across prey densities, functional response parameter values, experimental duration and predator assemblage type (conspecific or heterospecific). The population-dynamic model predictions sometimes diverged from
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Pairwise differences expressed as differences in the number of eaten prey between model predictions and simulated data for two conspecific predators for the multiplicative risk model, the direct FR model and the population-dynamic model (top, middle and bottom row, respectively) in experiments lasting 6, 12 and 24 hours (left, middle and right column, respectively) for different values of handling time *h*. Attack rate *a* = 2 in all simulations.</figcaption></figure>

the simulated data at high prey densities, but the difference was usually minor and much smaller than in the two other modelling approaches (Figs 2–5).

Discussion

Non-independent interactions among predators can have important consequences for the structure and dynam-
ics of ecological communities<sup>10,13,27</sup>. For instance, interference among predators decreases *per capita* trophic inter-
action strengths which, in turn, can increase community stability and dampen the destabilizing effect of resource
enrichment<sup>28</sup>. Quantifying emergent MPEs and determining when and how they change along ecological gra-
dients (e.g., resource density) is thus central for a full understanding of the effects of natural and anthropogenic
pressures on food web structure and dynamics<sup>6,9–11,13,17</sup>. As the quantification of MPEs relies on the comparison
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Pairwise differences expressed as differences in the number of eaten prey between model predictions and simulated data for two conspecific predators for the multiplicative risk model, the direct FR model and the population-dynamic model (top, middle and bottom row, respectively) in experiments lasting 6, 12 and 24 hours (left, middle and right column, respectively) for different values of attack rate *a*. Handling time *h* = 0.005 in all simulations.</figcaption></figure>

of empirical observations (prey survival under predation by multiple predators) with predictions from a model assuming independent predator individual effects, the accuracy and robustness of these predictions are crucial to assess the sign and magnitude of MPEs. In this study, we investigated the strengths and limitations of three modelling approaches (the popular multiplicative risk model, the direct FR model, and the population-dynamic model) and compared their predictions using simulated data sets. Our findings emphasize the importance of accounting for the non-linearity and time-dependence of species interactions to improve the detection and quantification of emergent multiple predator effects.
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Pairwise differences expressed as differences in the number of eaten prey between model predictions and simulated data for two heterospecific predators for the multiplicative risk model, the direct FR model and the population-dynamic model (top, middle and bottom row, respectively), experimental time of 6, 12, or 24 hours (left, middle and right column, respectively), and different values of handling time *h*₂ (see inside figure). Other functional response parameters (handling time *h* and attack rate *a*) fixed at *a*₁ = *a*₂ = 2 and *h*₁ = 0.005.</figcaption></figure>

**Model comparisons based on simulated data.** The multiplicative risk model proposed by Soluk¹⁸ has been used extensively in the ecological literature to quantify MPEs. However, McCoy *et al.*¹³ showed that this model can lead to false positives and detect MPEs even when predators affect prey independently. Our results based on simulated data corroborate those of McCoy *et al.*¹³ by showing that the multiplicative risk model tends to underestimate predation rates, especially when experimental duration is long, prey density is high and when the predators have short handling times and high attack rates (top rows in Figs 2–5). We also show that these qualitative conclusions apply to both conspecific and heterospecific predators (Figs 2–5).
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Pairwise differences expressed as differences in the number of eaten prey between model predictions and simulated data for two heterospecific predators for the multiplicative risk model, the direct FR model and the population-dynamic model (top, middle and bottom row, respectively), experimental time of 6, 12, or 24 hours (left, middle and right column, respectively) and different values of attack rate *a*₂ (see inside figure). Other functional response parameters (handling time *h* and attack rate *a*) fixed at *a*₁ = 2 and *h*₁ = *h*₂ = 0.005.</figcaption></figure>

Surprisingly, we have found that the predictions of the direct FR model and the multiplicative risk model are very similar for most of the explored parameter space. As a result, the systematic biases of the multiplicative risk model described above apply to the direct FR model as well: we found that this model tends to underestimate simulated predation rates by two conspecific or heterospecific predators (middle rows in Figs 2–5) except for very short experiments (Fig. S1). In the direct FR approach, a functional response model is fitted to the expected numbers of prey eaten by multiple predators. These expected numbers are generated using the multiplicative risk model, which explains the closely matching predictions of the two approaches. Based on our simulated data, we thus conclude that the direct FR approach suffers from similar biases as the multiplicative risk model and can
only be used when prey depletion is negligible and predators handle and digest prey rapidly (i.e., handling time is negligible). These restrictive conditions are likely to be violated in most experimental and natural systems as prey depletion is often substantial<sup>11,21</sup> and most predators have a nonlinear type II functional response (i.e., non-negligible handling times)<sup>26</sup>.

Moreover, the magnitude of the underestimation of the number of prey eaten in heterospecific predator assemblages produced by the multiplicative risk model and the direct FR model generally lie between the magnitudes observed in the two conspecific assemblages of these predators (Figs 2 and 3 vs Figs 4 and 5). This difference is linked to the fact that one of the two heterospecific predators has fixed handling time and attack rate values that are not at the limit of the range used in our simulations. In other words, the likelihood of having two predators with extreme values of handling time or attack rate is lower for heterospecific than for conspecific predator assemblages, which explains why model biases are less extreme in the former assemblage type.

To overcome the biases of the multiplicative risk model, McCoy *et al.*<sup>13</sup> proposed the population-dynamic model as an alternative approach based on functional response modeling that accounts for prey depletion and non-linear feeding rates. In a simulation study, they showed that the population-dynamic model produces unbiased predictions of the independent effect of predators on prey survival. In line with the results of McCoy *et al.*<sup>13</sup>, we found that the population-dynamic model predicts well the number of prey eaten by two conspecific or heterospecific predators (bottom rows in Figs 2–5). Although predictions of this model were less accurate at high prey densities, we found no systematic and substantial bias in the predictions of the population-dynamic model. This indicates that the predictions are robust against variations in functional response parameter values, experimental duration and composition of predator assemblages. Overall, we found that the population-dynamic model is more accurate and less biased than the multiplicative risk model and the direct FR model. We thus recommend using the former to predict the expected number of prey eaten by multiple predators under the assumption of predator independent effects.

**Implication for the detection and quantification of MPEs.** As mentioned in the Introduction, MPEs occur when the predictions based on single-predator feeding trials do not match observations from feeding trials involving multiple predators. The latter combine antagonistically (e.g., through predator interference) and synergistically (e.g., through predator facilitation) when the predictions overestimate and underestimate the observations, respectively. Our results indicate that the application of the multiplicative risk model and the direct FR model will likely lead to false positives (i.e., identify predators as synergistic) in the absence of MPEs because they systematically underestimate predation rates by multiple predators. This further implies that these two modelling approaches are less likely to detect antagonistic predator effects.

Even when the multiplicative risk model and the direct FR model get the MPEs qualitatively right, the two approaches will systematically overestimate synergistic predator effects and underestimate antagonist predator effects, especially at high prey densities or under more extensive prey depletion. Extensive prey depletion arises when predators have short handling time, high attack rates, or in longer experiments. Prey depletion is non-negligible in most predation studies<sup>21</sup> that last sufficiently long and do not replace eaten prey. Since 67% of the functional response experiments published so far were conducted over a 24-hour period and usually did not replenish eaten prey<sup>29</sup>, we conclude that the predictions of the direct FR model and multiplicative risk model will be unreliable in ‘typical’ MPE studies. In contrast, the population-dynamic model is more accurate and does not suffer from systematic biases. As a result, this model is much less likely to detect false positives (i.e., presence of MPEs) and the latter will not suffer from systematic biases towards synergetic or antagonistic effects.

Altogether, our results show that the detection and quantification of MPEs depends on the modelling approach used to predict the expected number of prey eaten by multiple predators. The choice of the modelling approach to detect MPEs can thus have serious consequences for our understanding of species interactions in complex food webs involving multiple predators.

**Model assumptions, biological interpretation of parameters, and possible limitations and extensions.** Functional response is a widely used mechanistic approach to model feeding rates, but its application beyond single consumers is not without problems and the need for careful and critical application has been emphasized repeatedly<sup>25,30,31</sup>. Functional response is a *per capita* rate, or in a more general sense, consumption rate per unit of consumer<sup>30</sup>. In other words, the functional response parameters (attack rate and handling time) are defined per individual predators. Their estimates from a single-predator functional response model fitted to multi-predator assemblages, as in the direct FR approach, lack mechanistic meaning and cannot be interpreted beyond providing a statistical fit to the data. In other words, these parameter estimates are aggregated attack rate and handling time values of a single predator, multiple conspecific predators or multiple heterospecific predators. They do not directly account for predator density or diversity and are thus incomparable across treatments. This prevents a meaningful interpretation of the functional response parameters and hampers the spatial or temporal extrapolation of the results, including the development of population-dynamic and food web models based on estimated functional response parameters.

To circumvent the limitations described above, the population-dynamic approach uses estimates of functional response parameters from single-predator feeding experiments to parameterise a simple population-dynamic model. These estimates are defined per predator, have a meaningful and mechanistic interpretation, and one can collect other functional parameter estimates from the literature (or from another experiment) to predict feeding rate by multiple predators using the population-dynamic approach, which is not possible with the direct FR approach or the multiplicative risk model. Moreover, the population-dynamic approach takes into account the non-linearity of feeding rates and prey depletion unlike the multiplicative risk model and predictions derived from it<sup>10,13</sup>. The population-dynamic approach can also cover situations in which the functional response shapes differ among predators (e.g., when one has a type II and another type III response), which goes beyond the direct
FR approach. For these reasons, the population-dynamic approach can be applied to more complex food-web models and longer timescales by, for instance, including prey and predator reproduction and tracking of population cycles.

The main limitation of the population-dynamic model compared to the multiplicative risk model is that the former requires estimating the functional response parameters of each predator species. This implies conducting single-predator feeding trials along a gradient of prey density, which can be time consuming and labour intensive. Using the multiplicative risk model requires single-predator feeding trials only for the experimental prey densities at which multiple-predator trials are conducted, which can make it much less labour intensive when the multiple-predator trials are conducted over only one or a few prey densities. The direct FR approach combines the limitations of both former approaches: it requires the same amount of data as the population-dynamic model but suffers from similar biases as the multiplicative risk model. We conclude that the most accurate approach is the population-dynamic model but if experimental workload is severely limiting, one may consider trading the workload against accuracy by using the multiplicative risk approach.

## Conclusions

Developing a robust and accurate modeling approach to quantify MPEs is crucial for advancing the study of MPEs and improving the realism of predator-prey and food-web models. In this study, we found that the predictions of the direct FR model and the multiplicative risk model are very similar and systematically underestimate predation rates by multiple predators at higher prey densities or when prey depletion is substantial, especially when experimental trials last longer than a few hours. Moreover, the functional response parameters estimated with the direct FR approach lack mechanistic interpretation, which limits the understanding of the mechanisms driving multiple predator interactions and prevents extensions of the direct FR approach to more complex food web structure or dynamics. Although it requires estimating the functional response of each predator species, our study indicates that the population-dynamic approach is more robust, precise, and informative than the multiplicative risk model and the direct FR approach. We thus recommend using the population-dynamic approach in future MPEs studies. In particular, studies using predators with different functional response shapes as well as experiments with substantial prey depletion should provide a critical test to the wider applicability of the population-dynamic approach to MPE studies.

## Material and Methods

**The direct FR approach.** This approach proposed by Wasserman et al.²⁰ consists of three main steps: (1) application of the multiplicative risk model (eqn 1) to data from single-predator trials to generate expected numbers of prey eaten by multiple predator, (2) fits of the expected data generated in the previous step and of data observed in multiple-predator trials with a functional response model, and (3) comparison of the fits based on the overlap of the confidence intervals (CIs) of both functional responses.

In the first step, their approach is to randomly assign a replicate number to each trial of the single-predator experiment and use these replicate numbers to generate predictions for multiple-predator treatments with the multiplicative risk model (eqn 1). For instance, replicate number 1 of a predator A feeding on a given prey density is associated with replicate number 1 of a predator B feeding on the same prey density to predict (using eqn 1) the combined predator A-predator B predation rate at the given prey density. For pairs of conspecific predators, their approach is to use the same replicate twice (note that this leads to pseudoreplication issues).

In the second step, after determining the functional response type (II or III) using standard methods²⁵, the data are fitted with the Rogers' random predator equation³² that accounts for prey depletion during the time course of the experiment:

$$N_e = N_0(1 - \exp(-a(t - hN_e))) \quad (2)$$

where $N_e$ is the number of prey eaten, $N_0$ is the initial prey density (prey.arena⁻¹), $t$ is the duration of the experiment (day), $h$ is the prey handling time (day.prey⁻¹), and $a$ is the predator attack rate (arena.day⁻¹), which is constant for type II and a function of $N_0$ for a type III functional response. Two separate instances of eqn 2 are fitted using a maximum-likelihood method; one to the predicted values for multiple-predator treatments generated at step 1 and another to the observed data from the multiple-predator trials. Non-integer predictions that often arise from the multiplicative risk model (eqn 1) need to be rounded to the nearest integer, because fitting eqn 2 using a maximum-likelihood method with a binomial distribution requires integer data³¹.

In the third step, datasets are non-parametrically bootstrapped (n = 2000 in the study by Wasserman et al.²⁰) and eqn 2 is fitted to each bootstrapped dataset using starting values of the functional response parameters ($a$ and $h$) that were obtained from the original maximum likelihood estimates (step 2). The bootstrapped fits are then used to construct 95% CIs around functional response curves using the “frair” package³³. Finally, functional response curves are visually compared based on the overlap of their 95% CIs; non-overlapping CIs are classified as significant differences between the functional response curves³³.

**The population-dynamic approach.** The approach proposed by McCoy et al.¹³ also consists of three main steps: (1) estimation of functional response parameters from single-predator trials, (2) use of these parameter estimates to predict the expected predation rate by multiple predators using a population dynamic model and (3) comparison of these predictions to the observed empirical values.

In the first step, after determining the functional response type (II or III) using standard methods²⁵, functional response parameters of each predator species are estimated by fitting eqn 2 to the single-predator experimental data using a maximum-likelihood method with a binomial distribution³¹.
In the second step, McCoy *et al.*¹³ proposed to use functional response parameter estimates from single-predator experiments (obtained in the first step) in a population-dynamic model of prey depletion by multiple predators in time:

$$ \frac{dN}{dt} = - \sum_{i=1}^{n} f_{i}(N)P_{i} \qquad (3) $$

where N is the prey population density, $P_i$ ($i=1, 2, ..., n$) are the population densities of predators i and $f_i(N)$ is the functional response of predator i. Equation 3 can accommodate any functional response shapes as, for instance, the Holling type II functional response, $f_i(N) = a_i N/(1 + a_i h_i N)$, which includes the estimates of prey handling times $h_i$ (day.prey⁻¹) and attack rates $a_i$ (arena.day⁻¹) of predators i obtained in the first step. This “null” model assumes no emergent MPEs and its predictions can be compared to multiple-predator feeding experiments to assess the sign and strength of MPEs. To generate predictions of expected prey survival in the multi-predator experiments, initial values of $N_0$ and $P_i$ ($i=1, 2, ..., n$) are set at the experimental initial prey and predator densities corresponding to the experimental treatment and the population dynamic model given by eqn 3 is integrated over the time interval of the experiment. The lower and upper 95% CIs around the predictions are estimated with a global sensitivity analysis that uses the 95% CIs of each functional response parameter estimate and their variance-covariance matrix (covariance is assumed to be zero when unknown) to generate a number of random parameter sets using a Latin hypercube sampling algorithm³⁴. For each parameter set, eqn 3 is then integrated over time and expected prey survival calculated using, e.g., the'sensRange' function in the R package 'FME'³⁴. The 2.5% and the 97.5% quantiles of the survival values obtained from these simulations are used as the lower and upper 95% CIs around the prediction.

For the third step, predictions of the population dynamic model are compared to the observed values based on the overlap of their 95% CIs for each predator treatments.

**Simulated data.** We simulated predation trials by numerically solving the population dynamic model (eqn 3) with a nonlinear type II functional response, which is the most common type of predator functional response²⁶. We used two types of simulations: deterministic to generate predictions of the three models based on single-predator trials, and stochastic to simulate 'realistic' trials including two predators to which we compared the model predictions.

In the deterministic simulations, we first generated data for single-predator trials (characterized by their species-specific handling time and attack rate values) along a gradient of prey density ranging from 1 to 800 prey.L⁻¹. For the functional response parameters used in the simulations, this prey range covers both the increasing and the asymptotic part of the functional response and covers a realistic range of prey densities such as the densities of *Daphnia* in fishless ponds in Central Europe (Sentis *et al.* unpublished data). We then (1) used these simulated data to feed the multiplicative risk model and predict cases with two predators under the assumption of pairwise independent effects, (2) fitted the functional response model to the predictions from the multiplicative risk model (i.e., the direct FR approach) and (3) used the population-dynamic model to generate predictions with two predators using the same set of parameters as the one used to generate data for the predator monoculture. Altogether, these different steps resulted in generating three predictions (one for each modelling approach) for each prey density about the expected multiple predator effects on prey survival under the assumption of predator independence, i.e., no emergent effects beyond prey depletion.

To estimate how well the different models can predict consumption rate by multiple predators under the assumption of predator independence, we used a stochastic approach to simulate predation trials for two conspecific or heterospecific predators by numerically solving eqn 3 using the same set of parameters as the one used to generate data for the predator monoculture. To mimic realistic experimental data and account for uncertainty in functional parameter estimates, we assumed that the parameters of the functional response, the handling time *h* and the attack rate *a*, are normally distributed with a standard deviation of respectively 5 and 10% of the mean parameter value reflecting typical parameter uncertainty values found in functional response experiments¹⁰,³⁵. To achieve this, we allowed values of parameters *a* and *h* to vary (following the normal distribution specified above) at each time step of the numerical integration of eqn 3, which generated variation in the simulated data. Eqn 3 was integrated 10 times for each set of parameters values to mimic empirical functional response data sets that often include up to 10 replicates per prey density¹⁰,²⁰,³⁵,³⁶. This yielded 10 replicates of simulated data (using 5 or 15 replicates did not substantially alter our conclusions, see Figs S2, S3).

We first illustrated the data for simulated trials including two predators and the predictions of the different modeling approaches for a selection of parameter values (see Fig. 1). We then systematically quantified pairwise differences between model predictions and the simulated data for two predators as the expected number of prey eaten (model prediction) minus the simulated number of prey eaten (mean of the 10 replicates). This yielded three pairwise differences, one between the multiplicative risk model and the simulated data, one between the direct FR approach and the simulated data, and one between the population-dynamic model and the simulated data. We used these differences as a measure of bias in the three modelling approaches; note that negative differences correspond to a situation in which the model predicts lower prey consumption by two predators than observed in the simulated trials.

We then determined how these differences changed with study duration *T*, handling time *h*, attack rate *a* and prey density *N*, all of which influence depletion and/or the nonlinearity of the functional response¹³,²⁴,²⁵,³². We used two predators ($P_1 = P_2 = 1$) and varied *T* from 1 to 24 hours in all simulations, which covers most of the range of experimental durations over which functional response experiment are conducted²⁹. We first assumed two conspecific predators (i.e., $a_1 = a_2 = a$ and $h_1 = h_2 = h$) and varied *h* from 0.002 to 0.02 while keeping *a* = 2. We next varied *a* from 1 to 8 while keeping *h* = 0.005. We then explored scenarios with two heterospecific predators
with a fixed handling time and attack rate values of the first predator (a₁ = 2 and h₁ = 0.005). We first varied h₂ from 0.002 to 0.02 while fixing the attack rate of the second predator (a₂ = 2) and then varied a₂ from 1 to 8 while fixing h₂ to 0.005. The range of values used for a and h were chosen to reflect values found in previous studies conducting both functional response experiment and multiple predator trials¹⁰,²⁰.

Simulations were performed with R version 3.4.1²⁸ using the package 'bbmle'³¹ to fit Roger's random predator equation and the packages 'deSolve'³⁷ and 'FME'³⁴ for the population dynamic model simulations.

**Data accessibility statement.** We confirm that the R code supporting the results will be archived in an appropriate public repository such as Dryad or Figshare.

**References**

1. McCann, K. S., Hastings, A. & Huxel, G. R. Weak trophic interactions and the balance of nature. *Nature* **395**, 794–798 (1998).

2. Berlow, E. L. Strong effects of weak interactions in ecological communities. *Nature* **398**, 330–334 (1999).

3. Tang, S., Pawar, S. & Allesina, S. Correlation between interaction strengths drives stability in large ecological networks. *Ecol. Lett.* **17**, 1094–1100 (2014).

4. Vázquez, D. P., Ramos-Jiliberto, R., Urbani, P. & Valdovinos, F. S. A conceptual framework for studying the strength of plant–animal mutualistic interactions. *Ecol. Lett.* **18**, 385–400 (2015).

5. Bascompte, J., Jordano, P. & Olesen, J. M. Asymmetric coevolutionary networks facilitate biodiversity maintenance. *Science* **312**, 431–433 (2006).

6. Schmitz, O. J. Predator diversity and trophic interactions. *Ecology* **88**, 2415–2426 (2007).

7. Werner, E. E. & Peacor, S. D. A review of trait-mediated indirect interactions in ecological communities. *Ecology* **84**, 1083–1100 (2003).

8. Bolker, B., Holyoak, M., Křivan, V., Rowe, L. & Schmitz, O. Connecting theoretical and empirical studies of trait-mediated interactions. *Ecology* **84**, 1101–1114 (2003).

9. Sih, A., Englund, G. & Wooster, D. Emergent impacts of multiple predators on prey. *Trends Ecol. Evol.* **13**, 350–355 (1998).

10. Sentis, A., Gémard, C., Jaugeon, B. & Boukal, D. S. Predator diversity and environmental change modify the strengths of trophic and non-trophic interactions. *Global Change Biol.* **23**, 2629–2640 (2017).

11. Vance-Chalcraft, H. D., Rosenheim, J. A., Vonesh, J. R., Osenberg, C. W. & Sih, A. The influence of intraguild predation on prey suppression and prey release: a meta-analysis. *Ecology* **88**, 2689–2696 (2007).

12. Okuyama, T. & Bolker, B. M. In *Trait-mediated indirect interactions: ecological and evolutionary perspectives* (eds Takayuki Ohgushi, Oswald Schmitz, & Robert D Holt) 186–204 (Cambridge University Press, 2012).

13. McCoy, M. W., Stier, A. C. & Osenberg, C. W. Emergent effects of multiple predators on prey survival: the importance of depletion and the functional response. *Ecol. Lett.* **15**, 1449–1456 (2012).

14. Denny, M. & Benedetti-Cecchi, L. Scaling up in ecology: mechanistic approaches. *Annu. Rev. Ecol., Evol. Syst.* **43**, 1–22 (2012).

15. Snyder, W. E., Snyder, G. B., Finke, D. L. & Straub, C. S. Predator biodiversity strengthens herbivore suppression. *Ecol. Lett.* **9**, 789–796 (2006).

16. Losey, J. E. & Denno, R. F. Positive predator-predator interactions: enhanced predation rates and synergistic suppression of aphid populations. *Ecology* **79**, 2143–2152 (1998).

17. Sentis, A., Hemptinne, J. L. & Brodeur, J. Towards a mechanistic understanding of temperature and enrichment effects on species interaction strength, omnivory and food-web structure. *Ecol. Lett.* **17**, 785–793 (2014).

18. Soluk, D. A. Multiple predator effects: predicting combined functional response of stream fish and invertebrate predators. *Ecology* **74**, 219–225 (1993).

19. Griffen, B. Detecting emergent effects of multiple predator species. *Oecologia* **148**, 702–709 (2006).

20. Wasserman, R. J. *et al.* Using functional responses to quantify interaction effects among predators. *Funct. Ecol.* **30**, 1988–1998 (2016).

21. Vance-Chalcraft, H. D. & Soluk, D. A. Multiple predator effects result in risk reduction for prey across multiple prey densities. *Oecologia* **144**, 472–480 (2005).

22. Vonesh, J. R. & Osenberg, C. W. Multi-predator effects across life-history stages: non-additivity of egg-and larval-stage predation in an African treefrog. *Ecol. Lett.* **6**, 503–508 (2003).

23. Collins, M. D., Ward, S. A. & Dixon, A. F. G. Handling time and the functional response of *Aphelinus thomsoni*, a predator and parasite of the aphid *Drepanosiphum planatoidis*. *J. Anim. Ecol.* **47**–487 (1981).

24. Sentis, A., Hemptinne, J. L. & Brodeur, J. How functional response and productivity modulate intraguild predation. *Ecosphere* **4**, 1–14 (2013).

25. Juliano, S. A. In *Design and analysis of ecological experiments* (eds Scheiner, S.M. & Gurevitch, J.) 178–196 (Chapman and Hall, 2001).

26. Jeschke, J. M., Kopp, M. & Tollrian, R. Consumer-food systems: why Type I functional responses are exclusive to filter feeders. *Biological Reviews* **79**, 337–349 (2004).

27. Rall, B. C., Guill, C. & Brose, U. Food web connectance and predator interference dampen the paradox of enrichment. *Oikos* **117**, 202–213 (2008).

28. R Development Core Team. *R: a language and environment for statistical computing*. (2017).

29. Li, Y., Rall, B. C. & Kalinkat, G. Experimental duration and predator satiation levels systematically affect functional response parameters. *Oikos*, https://doi.org/10.1111/oik.04479 (2018).

30. Koen-Alonso, M. In *From energetics to ecosystems: the dynamics and structure of ecological systems* (eds Neil Rooney, Kevin S. McCann & David L. G. Noakes) 1–36 (Springer, 2007).

31. Bolker, B. M. *Ecological models and data in R*. (Princeton University Press, 2008).

32. Rogers, D. Random search and insect population models. *J. Anim. Ecol.* **41**, 369–383 (1972).

33. Pritchard, D. W., Paterson, R. A., Bovy, H. C. & Barrios-O'Neill, D. frair: an R package for fitting and comparing consumer functional responses. *Methods in Ecology and Evolution* **8**, 1528–1534 (2017).

34. Soetaert, K. & Petzoldt, T. Inverse modelling, sensitivity and Monte Carlo analysis in R using package FME. *Journal of Statistical Software* **33**, 1–28 (2010).

35. Englund, G., Ohlund, G., Hein, C. L. & Diehl, S. Temperature dependence of the functional response. *Ecol. Lett.* **14**, 914–921 (2011).

36. Sentis, A., Hemptinne, J. L. & Brodeur, J. Using functional response modeling to investigate the effect of temperature on predator feeding rate and energetic efficiency. *Oecologia* **169**, 1117–1125 (2012).

37. Soetaert, K., Petzoldt, T. & Setzer, R. W. Solving differential equations in R: package deSolve. *Journal of Statistical Software* **33**, 1–25 (2010).
## Acknowledgements

We thank Viktoriia Radchuk and two anonymous reviewers for detailed comments that substantially improved the manuscript. We thank Ryan J. Wassermann for providing information on their experimental protocol and Grant Agency of the Czech Republic (project no. 14–29857 S) for financial support. AS was also founded by the French Laboratory of Excellence projects ‘LABEX TULIP’ (ANR-10-LABX-41), ANR funded Toulouse Initiative of Excellence “IDEX UNITI” (ANR11-IDEX-0002-02) and by the People Programme (Marie Curie Actions) of the European Union’s Seventh Framework Programme (FP7/2007–2013) under REA grant agreement n. PCOFUND-GA-2013-609102, through the PRESTIGE programme coordinated by Campus France.

## Author Contributions

A.S. conceived and developed the functional response models with the contribution of D.S.B., A.S. wrote the first draft of the manuscript, and D.S.B. contributed substantially to revisions.

## Additional Information

**Supplementary information** accompanies this paper at https://doi.org/10.1038/s41598-018-30244-9.

**Competing Interests:** The authors declare no competing interests.

**Publisher's note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_9.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International

License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2018