Incorporating nonlinearity with generalized functional responses to simulate multiple predator effects

Generalized functional response framework

The generalized functional response framework integrates nonlinear predator functional responses, nonlinear modifiers of risks (e.g., size- or temperature-dependent predation), and changes in somatic growth to generate null expectations for the combined effects of any number of independently acting predators. Specifically, we use a partial differential equation model that allows predation to change as a function of both prey density and a nonlinear modifier (s) of predation rates (E_i) by k different predator types (P_i); and when relevant, it can also account for changes in size-dependent and density-dependent prey somatic growth (g).

(1) ${\displaystyle \frac{\mathrm{\partial }C\left(s,t\right)}{\mathrm{\partial }t}=-C{\sum }_{i=1}^k{E}_i\left(C,s,P_i\right)+\left(-g\left(s\right){\displaystyle \frac{\mathrm{\partial }C}{\mathrm{\partial }s}}\right)+D\left({\displaystyle \frac{{\mathrm{\partial }}^{2}C}{\mathrm{\partial }{s}^2}}\right)}$

In this equation, the vector $\mathit{C}$ denotes the distribution of prey abundances across levels of s. For simplicity, we assume that predator body sizes within a given species are approximately constant over time.

We assume prey consumption (E_i) by each predator (P_i ) follows a Type II saturating functional response (Holling, 1959; Lafferty et al., 2015; McCoy, Stier & Osenberg, 2012; Rogers, 1972), which describes the density-dependent foraging dynamics of most natural enemies (Jeschke, Kopp & Tollrian, 2004; Lafferty et al., 2015). We incorporate modifier-dependent predation by modeling the attack rate, $a_i$, and handling time, $h_i$, as functions of a prey risk modifier $s$. Details about the derivation of our size distributed functional response are in Supplement 1.

The total attack rate from predator type i (predators act independently, i.e., there is no interference competition) is

(2) $E_i={\displaystyle \frac{a_i\left(s\right)C\left(s\right)}{1+A_i{H}_i},}$where A_i = ${\sum }_{s=min}^{max}{a}_i\left(s\right)C\left(s\right)$ is the total attack rate for predator type i (i.e., the expected per-predator rate of prey consumption when prey are rare or handling is not limiting) and H_i is the propensity-weighted average handling time, i.e.,

$$H_i={\displaystyle \frac{{\sum }_{s=min}^{max}{a}_i\left(s\right)C\left(s\right)h_i\left(s\right)}{A_i}}$$

These expressions take into account the fact that prey-handling by predators is limited by the entire range of prey they consume, not just the prey in a particular class (Supplement 1).

A wide array of potential processes and functional forms can be used to model modifier-dependent predation risk (e.g., McCoy et al., 2011). For instance, temperature is known to modify predation in nonlinear ways and could be modeled by making the attack rate a quadratic function of temperature and handling time a power function of temperature and predator body mass (Davidson et al., 2021; Sentis & Boukal, 2018; Uiterwaal & DeLong, 2020).

Here we focus on body size as an exemplar of a modifier that has nonlinear effects on predation risk. Size-specific risk is often a nonlinear function of predator and prey body sizes and can sometimes have stronger implications for per-capita predation risk than prey density (McCoy et al., 2011; McCoy & Bolker, 2008; Rudolf, 2008; Woodward et al., 2005). Therefore, failure to account for variation in the sizes of prey in experiments (e.g., differences in prey size structure), or changes in prey size as a result of somatic growth over the course of an experiment (e.g., a study of a single cohort over time), can bias predictions about the combined effects of multiple predators. Most species increase in size by orders of magnitude during their lifetimes (Fenchel, 1974; Werner & Gilliam, 1984); thus, failure to appreciate the dynamic effects of body size on predator-prey interactions can complicate our understanding of predator prey dynamics.

We will explore two empirically supported functions for changes in predation risk as functions of prey size, an exponential function (Alford, 1999; Aljetlawi, Sparrevik & Leonardsson, 2004; McCoy et al., 2011; McCoy & Bolker, 2008) of the form

(3) $$a=\alpha \cdot e^{\left(1-{\displaystyle \frac{s}{d}}\right)}$$and a unimodal function (Kalinkat et al., 2013; McCoy et al., 2011; McCoy & Bolker, 2008; Vonesh & Bolker, 2005; Vucic-Pestic et al., 2010; Wahlström et al., 2000) of the form

(4) $$a=\alpha {\left({\displaystyle \frac{s}{d}\cdot e^{\left(1-\left({\displaystyle \frac{s}{d}}\right)\right)}}\right)}^{\gamma }$$where $\alpha$ is maximum attack rate, $d$ scales the most vulnerable prey size s, and $\gamma$ determines the width of the window of maximum vulnerability (i.e., the level of prey size specialization for a predator) (Fig. 1). For simplicity, we will assume predator size is not variable, however predator size or predator:prey body size ratios could also be considered within this framework. We assume handling time is independent of prey body size, but this assumption could be relaxed (e.g., Eq. (2) and McCoy et al., 2011).

When predation is size-dependent we may also need to consider prey body size growth through time. We represent prey growth in our generalized functional response formulation via an advection term that moves the prey through different size classes. In order to integrate these equations over time in a numerically stable way, we implement the advection term by assuming the locations of size bins (rather than sizes of individuals) change over time. We first solve the growth equation ( $ds/dt=g\left(s\right)$) to find size as a function of time, then substitute those size values into the predation term at each time step. As with size-dependent predation, many functional forms are available for modeling growth; we use a standard Gompertz growth model. This model assumes that relative growth rate decreases exponentially with size, leading to a growth curve

(5) $$s\left(t\right)=\rho \cdot e^{-\theta \cdot e^{\left(-\delta t\right)}}$$where $\rho \mathrm{i}\mathrm{s}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{m}\mathrm{a}\mathrm{x}\mathrm{i}\mathrm{m}\mathrm{u}\mathrm{m}\mathrm{s}\mathrm{i}\mathrm{z}\mathrm{e}\mathrm{a}\mathrm{n}\mathrm{d}\theta$ is the negative log of scaled initial size (i.e., $\theta =-\mathrm{l}\mathrm{n}\left(s\left(0\right)/\rho \right)$) and $\delta$ is a constant that describes the rate of decrease in relative growth rate as size decreases. For example, if the first size bin initially includes individuals in the range [s1(0),s2(0)], at time t it will include individuals in the range [s1(t),s2(t)].

We also incorporate individual variation in growth via a diffusion term, the third term in Eq. (1). This term describes diffusion along the size axis that naturally arises from continuously occurring, independent variation around growth rate (i.e., Brownian motion Brooks, McCoy & Bolker, 2013). Because the width of the size bins in the model varies over time according to the specified growth curve, the effective diffusion rate is also size-dependent.

For the purposes of this study, we focus on two predator-one prey interactions with a few specific functional forms. However, Eq. (1) is generalizable to model assemblages consisting of many (i = 1 to k in Eq. (1)) concurrent predators and a wide array of functional forms for the predator functional response (e.g., Holling type III, Holling, 1959), risk modifiers, and somatic growth models (Gompertz, 1825; Kahm et al., 2010; Richards, 1959; West, Brown & Enquist, 2001).

We conduct a literature review to examine the prevalence of size dependent interactions in multiple predator effects studies. Next, we simulate predation by two predators using our generalized functional response model with prey size-dependent predation risk and compare the simulation results with predictions generated by the multiplicative risk model (MRM). We then explore how different functional forms of size-dependent risk influence the MRM predictions relative to simulation outputs. To contextualize our analysis we used functional relationships and estimates of size variation reported in studies of size dependent predation on red eyed treefrog larvae (Brooks, McCoy & Bolker, 2013; McCoy et al., 2011).

Reviewing size dependence

We searched web of science on May 8, 2018 using the terms “(“multiple pred*”) OR (“risk reduction” & pred*) OR (“risk enhancement” & pred*) OR (MPE & pred*) & experiment”, and retained results for the previous 20 years, a total of 492 papers. We then screened the studies to include only those that were experimental, manipulated presence of at least two predators, and had treatments of control, predator monocultures, and multiple predators foraging together. After screening, we retained 121 observations from 119 studies (Supplements 2 and 3). We then reviewed the full text of these 121 studies to determine how or if prey size was reported and if so, how much size changed over the duration of the experiment.

Comparing the MRM predictions and generalized functional response simulations

The MRM predicts the expected proportion of prey surviving in the presence of two independent predators as the product of survival (proportion alive) in the presence of each predator alone ( $S_i$ and $S_j$) corrected for survival in the absence of a predator ( $S_c$) in additive

(6) $${\hat{S}}_{ij}={\displaystyle \frac{S_i\cdot S_j}{S_c}}$$or substitutive

(7) $${\hat{S}}_{{\displaystyle \frac{ij}{2}}}={\displaystyle \frac{S_i^{0.5}\cdot S_j^{0.5}}{S_c}}$$experimental designs (Billick & Case, 1994; Griffen, 2006; McCoy, Stier & Osenberg, 2012; Sih, Englund & Wooster, 1998; Soluk & Collins, 1988; Vonesh & Osenberg, 2003). Predicted predation rates based on the MRM assume that per-capita risk is the same for all prey individuals and that it remains constant over time (McCoy, Stier & Osenberg, 2012; Sih, Englund & Wooster, 1998). However, when predators have nonlinear and saturating functional responses (Jeschke, Kopp & Tollrian, 2004; Lafferty et al., 2015) and prey are depleted, or when risk is a nonlinear function of some other modifying factor (e.g., temperature or size), per capita risk is not constant (Juliano, Scheiner & Gurevitch, 2001; McCoy, Stier & Osenberg, 2012; Rogers, 1972). Thus, applications of the MRM based on static estimates of survival can bias estimates of per-capita risk from individual predators, which carry over to bias predictions of combined predators (McCoy, Stier & Osenberg, 2012; Sentis & Boukal, 2018). By accounting for nonlinear modifiers of prey risk and depletion, we expect that the generalized functional response framework presented here could be used to generate more accurate predictions for the combined effects of multiple predators.

We simulated foraging trials for predators in monoculture and in mixed predator species assemblages by numerically solving Eq. (1) (example code available in Supplement 4), providing estimates of prey survival for each of two predators alone and in combination (i.e., S₁, S₂, and S_1,2). Estimates of S_1,2 from these simulations provided the null model for independent interactions between two predators with size-dependent predation rates. We only simulate additive combinations of predators; thus, all comparisons are based on predictions from the MRM described by Eq. (6) and the generalized functional response in Eq. (1). The degree of mismatch defined as the rate difference between generalized functional response simulations and the predictions from the MRM (Eq. (6)) provides a metric for the potential bias in predictions when prey are being depleted and predation risk is size-dependent. In other words, differences between simulation results and MRM predictions are caused by a failure of the MRM to account for the non-linear, but fundamentally predictable (and therefore not emergent), effects of density- and size-dependent predation.

For a broad range of plausible parameter values for each predator, we explored how failure to account for size-dependence can bias inferences about multiple predator effects. For attack rates that decayed exponentially with prey size, we varied the size scaling parameter $d$, which describes how risk decreases with prey size (Fig. 1A). For attack rates that were unimodal functions of size, we varied two key parameters regulating size-dependent risk, while keeping all other parameters constant: $d$, the size scaling parameter that determines the most vulnerable prey size (Fig. 1B), and $\gamma$ the parameter that determines the degree of size specialization (i.e., width of the window of maximum vulnerability) (Fig. 1C). Because each of these parameters affects size-dependent risk due to their independent and interactive effects, we explored how they affected the accuracy of predictions in three ways. First, to examine the interactive effects of the unimodal shape parameters, we set one predator to be a size generalist (i.e., a low value $\gamma$; Eq. (4)) that was most effective at killing small prey (i.e., a low value of d; Eq. (4)) and then allowed the second predator to have increasingly higher values of these two parameters (i.e., increasingly more size specialized and more effective in consuming larger prey) in each simulation. Second, to explore the independent effects of the values of the size scaling (d) and size specialization ( $\gamma$) parameters on bias, we fixed both predators to have either high or low values of each parameter (i.e., to be more efficient consumers of large or small prey, d; or to be size specialist or generalist, $\gamma$). Third, we allowed the value of the free parameter to vary between the predators.

Given the sparse data on prey size in studies of MPEs, we grounded our simulations by using prey sizes and growth rates published for one particular empirical system (McCoy et al., 2011) along with attack rates and handling times to match those used in McCoy, Stier & Osenberg (2012) (i.e., maximum attack rate $\alpha$ = 0.5, and handling time $h$ = 1). We initiated all simulations with a prey population size of 100 individuals evenly divided into 10 prey size classes ranging between 1 and 45 mm. We varied the size scaling parameter d for attack rates that decayed exponentially with size (via Eq. (3)). We varied the value of $\gamma$ between 0.001 and 3.5, and the value of $d$ from 3.5 to 24 (see table 2 in McCoy et al., 2011). For simulations of growing cohorts we set maximum size $\rho$ (Eq. (5)) at 35 mm, the maximum growth rate parameter $\theta$ to 0.05 and the growth rate decay parameter $\delta$ to a value of −0.05 (Eq. (5)), which are within the ranges of estimates of tadpole growth reported in McCoy et al. (2011) and Albecker & McCoy (2019).

Prey size refuge

When predation risk decayed exponentially with prey size, the MRM predicted either higher or lower risk than expected from the generalized functional response model depending upon the specific parameters of the risk function (Fig. 1A) and the experimental duration (Figs. 2A and 2B). For small values of d (i.e., strong size dependence for small prey) the MRM predicted higher levels of risk than expected (i.e., predicted lower survival), but this pattern was reversed for larger values of d (i.e., strong size dependence for larger prey) with the MRM predicting lower risk (higher survival) than expected (Fig. 2). This pattern emerged because the MRM overestimates predation risk for larger prey when risk is strongly size-dependent and fails to capture the higher combined size specific rates of depletion for small and intermediate sized prey when the rate of risk decay is shallower (i.e., larger values of d). These patterns were similar for growing cohorts and size structured prey, but with the magnitude of bias for growing cohorts strongly regulated by prey growth rates.

Size specialization

When size dependence was a unimodal function of prey size and small prey were most vulnerable (i.e., low values for d – Fig. 1B), the MRM predicted lower survival as the predators became more specialized relative to the generalized functional response predictions (Fig. 3A). However, this pattern reversed when either predator was a size generalist, with the MRM predicting higher survival than the generalized model (Fig. 3A). Interestingly, the MRM predicted higher survival over a much broader parameter space as the most vulnerable size class increased, with the pattern only reversing for the most size-specialized predators (Figs. 3B and 3C). This pattern is likely due to the increasing influence of depletion on per capita risk, as a large range of size classes are vulnerable to predation (McCoy, Stier & Osenberg, 2012). The strong effects of specializing on small prey also emerge in Fig. 4. When predators do not have strong size preferences (Fig. 4A; small γ) the MRM predicts higher survival for all values of d because it fails to account for the increasing risk that occurs with depletion (McCoy, Stier & Osenberg, 2012). However, this pattern reverses as the predator size specialization becomes more humped shaped and peaked, because the MRM vastly over estimates mortality of the smallest and largest prey (Figs. 4B and 4C).