Rejuvenating functional responses with renewal theory

TL;DR

This paper introduces a mechanistic, stochastic framework based on Renewal Theory to derive and analyze ecological functional responses from individual behaviors, improving understanding and empirical inference.

Contribution

It presents a bottom-up approach linking individual traits to population-level functional responses using Renewal Theory, with applications for modeling and data analysis.

Findings

Explicit formulas for mean and variance of functional responses.

Application to specific ecological contexts like nuptial-feeding species.

Method for inferring parameters from interaction variability.

Abstract

Functional responses are widely used to describe interactions and resources exchange between individuals in ecology. The form given to functional responses dramatically affects the dynamics and stability of populations and communities. Despite their importance, functional responses are generally considered with a phenomenological approach, without clear mechanistic justifications from individual traits and behaviors. Here, we develop a bottom-up stochastic framework grounded in Renewal Theory showing how functional responses emerge from the level of the individuals through the decomposition of interactions into different activities. Our framework has many applications for conceptual, theoretical and empirical purposes. First, we show how the mean and variance of classical functional responses are obtained with explicit ecological assumptions, for instance regarding foraging behaviors. Second, we give examples in specific ecological contexts, such as in nuptial-feeding species or size dependent handling times. Finally, we demonstrate how to analyze data with our framework, especially highlighting that observed variability in the number of interactions can be used to infer parameters and compare functional response models.

Figures

Tables

Table 1. Table 1: Examples of assumptions for recovering some classical functional responses in a 1D space with direct movement to the nearest individual (from Eq. 8 ). For the sake of simplicity, formulas are given under the assumption that y >> 1 much-greater-than 𝑦 1 y>>1 , α ​ z >> 1 much-greater-than 𝛼 𝑧 1 \alpha z>>1 and β ​ z >> 1 much-greater-than 𝛽 𝑧 1 \beta z>>1 , and denoting l = λ / 2 𝑙 𝜆 2 l=\lambda/2 . ( × Δ − 1 absent superscript Δ 1 \times\Delta^{-1} ) means that the variance is scaled with the parameter Δ Δ \Delta (see Eq. 2 )

Functional	Handling	Searching	Interaction	Mean
Response	time	probability	with $e_z$	& Variance ($\times {\Delta }^{-1}$)
$R(.)$	$c_y$	$p_s$
No handling	$0$	$1$	No: $z=0$	$y{l}^{-1}$
2 species
(Holling Type I)				$y{\left(3l\right)}^{-1}$
With Handling	$>0$	$1$	No: $z=0$	$y{\left[l+y{c}_y\right]}^{-1}$
2 species
(Holling Type II)				$y{l}^2{\left[3\left(l+y{c}_y\right)\right]}^{-3}$
With Handling	$>0$	$1$	Yes: $z>0$	$y^2{\left[\alpha z\left(c_z-ly\right)+y^2{c}_y\right]}^{-1}$
3 species			$c_z\to c_z/y$
(Holling Type III)			$\beta =0$	$y^{5}z\alpha l^2{\left[3\left(\alpha z\left(c_z-ly\right)+y^2{c}_y\right)\right]}^{-3}$
With Handling	$>0$	$1$	Yes : $z=x$	$y{\left[2l+c_{y}y+c_x\alpha x\right]}^{-1}$
2 species			$\alpha =\beta$
Predators Interference				$y{l}^2\left(y+\alpha x\right){\left[3\alpha x\left(2l+c_{y}y+c_x\alpha x\right)\right]}^{-3}$
(Beddington-DeAngelis)
With Handling	$>0$	$x^{-1}$	No: $z=0$	$y/x{\left[l+c_{y}y/x\right]}^{-1}$
2 species
Predators Competition				$y/x{l}^2{\left[l+c_{y}y/x\right]}^{-3}$
(Ratio Dependence)

Equations

$$N{}_{\Delta }=k\text{if}{T}_1+\cdots +T_k\le \Delta <T_1+\cdots +T_k+T_{k+1}.$$

$$R(x,y,z)=\frac{N{}_{\Delta }(x,y,z)}{\Delta }.$$

$$R(x,y,z)\stackrel{\mathrm{law}}{\approx }\frac{1}{\text{mathds}E\left(T(x,y,z)\right)}+\mathcal{N}\left(0,\frac{1}{\Delta }\frac{\mathrm{Var}\left(T(x,y,z)\right)}{\text{mathds}E{\left(T(x,y,z)\right)}^3}\right),$$

$$T=\sum _{a=1}^{|\mathcal{A}|}{\tau }_a^{\Sigma }$$

$${\tau }_a^{\Sigma }=\sum _{i=1}^{G_a}{\tau }_a^{(i)},$$

$$\text{mathds}E(T)=\sum _{a=1}^{|\mathcal{A}|}\frac{\text{mathds}E({\tau }_a)}{p_a}$$

$$\mathrm{Var}\left(T\right)=\sum _{a=1}^{|\mathcal{A}|}\left(\frac{\mathrm{Var}\left({\tau }_a\right)}{p_a}+\frac{1-p_a}{p_a^2}\text{mathds}E{\left({\tau }_a\right)}^2\right).$$

$$\text{mathds}E(T)=\frac{1}{p_s}\text{mathds}E({\tau }_s)+\text{mathds}E({\tau }_h)$$

$$\mathrm{Var}\left(T\right)=\frac{\mathrm{Var}\left({\tau }_s\right)}{p_s}+\frac{1-p_s}{p_s^2}\text{mathds}E{\left({\tau }_s\right)}^2+\mathrm{Var}\left({\tau }_h\right)$$

$$\text{mathds}E\left(D(y)\right)$$

$$\text{mathds}E\left(D(y{)}^2\right)$$

$$\text{mathds}E\left(R\right)$$

$$\mathrm{Var}\left(R\right)$$

$$\text{mathds}E\left(T\right)$$

$$\mathrm{Var}\left(T\right)$$

$$\text{mathds}E\left(T\right)$$

$$\mathrm{Var}\left(T\right)$$

$$\text{mathds}E\left(T\right)=\frac{1}{p_R}\left(\frac{\text{mathds}E\left({\tau }_{S_G}\right)}{p_{S_G}}+\text{mathds}E\left({\tau }_{H_G}\right)+\frac{\text{mathds}E\left({\tau }_{S_F}\right)}{p_{S_F}}+\frac{\text{mathds}E\left({\tau }_C\right)}{p_C}\right)$$

$$\mathrm{Var}\left(T\right)=\frac{1-p_R}{p_R^2}{\left(\frac{\text{mathds}E\left({\tau }_{S_G}\right)}{p_{S_G}}+\text{mathds}E\left({\tau }_{H_G}\right)+\frac{\text{mathds}E\left({\tau }_{S_F}\right)}{p_{S_F}}+\frac{\text{mathds}E\left({\tau }_C\right)}{p_C}\right)}^2$$

$$+\frac{1}{p_R}(\frac{1-p_C}{p_C^2}\text{mathds}E{\left({\tau }_C\right)}^2+\frac{1-p_{SG}}{p_{SG}^2}\text{mathds}E{\left({\tau }_{SG}\right)}^2+\frac{1-p_{SF}}{p_{SF}^2}\text{mathds}E{\left({\tau }_{SF}\right)}^2$$

$$+\frac{\mathrm{Var}\left({\tau }_{S_G}\right)}{p_{S_G}}+\frac{\mathrm{Var}\left({\tau }_{S_F}\right)}{p_{S_F}}).$$

$$\text{mathds}E({\tau }_h)$$

$$\mathrm{Var}({\tau }_h)$$

$$N_t=\sum _{n=1}^{+\infty }\text{mathds}{1}_{\{S_n\le t\}}=\sup \{n\ge 1:S_n\le t\}.$$

$$\frac{N_t}{t}\underset{t\to +\infty }{\to }\frac{1}{\text{mathds}E(T)}\text{with probability one}.$$

$$\frac{S_{N_t}}{N_t}\le \frac{t}{N_t}<\frac{S_{N_t+1}}{N_t}.$$

$$\frac{S_{N_t}}{N_t}\underset{t\to +\infty }{\to }\text{mathds}E(T)\text{and}\frac{S_{N_t+1}}{N_t+1}\frac{N_t+1}{N_t}\underset{t\to +\infty }{\to }\text{mathds}E(T).$$

$$\frac{N_t-\frac{t}{\text{mathds}E(T)}}{\sigma \sqrt{\frac{t}{(\text{mathds}E(T){)}^3}}}\underset{t\to +\infty }{\overset{\mathrm{law}}{\to }}\mathcal{N}(0,1),$$

$$\text{mathds}P\left(\frac{N_t-\frac{t}{\text{mathds}E(T)}}{\sigma \sqrt{\frac{t}{(\text{mathds}E(T){)}^3}}}<y\right)$$

$$=\text{mathds}P(S_{N_t}>t)\text{(because}\{N_t<n\}=\{S_n>t\}\text{)}$$

$$=\text{mathds}P\left(\frac{S_{n_t}-n_t\text{mathds}E(T)}{\sigma \sqrt{n_t}}>\frac{t-n_t\text{mathds}E(T)}{\sigma \sqrt{n_t}}\right).$$

$$\underset{t\to \infty }{\lim }\frac{t-n_t\text{mathds}E(T)}{\sigma \sqrt{n_t}}=\underset{t\to \infty }{\lim }\frac{t-r_t\text{mathds}E(T)}{\sigma \sqrt{r_t}}=\underset{t\to \infty }{\lim }\frac{-y\text{mathds}E(T)\sqrt{\frac{t}{(\text{mathds}E(T){)}^3}}}{\sqrt{\frac{t}{\text{mathds}E(T)}+y\sigma \sqrt{\frac{t}{(\text{mathds}E(T){)}^3}}}}=-y.$$

$$\underset{t\to +\infty }{\lim }\text{mathds}P\left(\frac{N_t-\frac{t}{\text{mathds}E(T)}}{\sigma \sqrt{\frac{t}{(\text{mathds}E(T){)}^3}}}<y\right)=\text{mathds}P(Z>-y)=\text{mathds}P(Z<y).$$

$$\frac{N{}_{\Delta }(x,y,z)}{\Delta }-\frac{1}{\text{mathds}E\left(T(x,y,z)\right)}$$

$$\frac{1}{\Delta }\frac{\mathrm{Var}\left(T(x,y,z)\right)}{(\text{mathds}E\left(T(x,y,z)\right){)}^3},$$

Full text

Running title. Stochastic derivation of functional responses

Keywords. stochastic model, ecology, predation, cooperation, inference, emerging properties

Corresponding author: Sylvain Billiard: sylvain.billiard@univ-lille1.fr, Univ. Lille, CNRS, UMR 8198 - Evo-Eco-Paleo, F-59000 Lille, France.

1 Introduction

Interactions between individuals affect all ecological processes, and how fast they occur determine resources exchange rates in an ecosystem. Interactions are generally considered at the population (or macroscopic) level and supposed to vary along with species densities following an interaction function, generally called a functional response. The seminal papers introducing functional responses in population ecology followed a reductionistic approach and aimed at giving the underlying mechanisms [1, 2, 3]. Yet, population ecologists generally follow a phenomenological approach to justify functional responses [4], and functional responses are rarely derived from the individuals’ traits and behaviors. The form given to functional responses is crucial since it can dramatically affect the dynamics and stability of populations and communities. For instance, the stability of predator and prey populations strongly depends on whether predator consumption rates increases linearly (Holling type I functional response) or following a saturating function (Holling Type II and III functional responses) with prey densities [5]. It is thus critical to modelize within and between species interactions as accurately as possible in order to have the best predictions and understanding of population and community dynamics, and eventually support wildlife management decisions [6].

The form given to functional responses is largely debated for decades, sometimes fiercely (e.g. the long-standing controversy about whether it is best to assume density-dependence or ratio-dependence in predator-prey models [7, 8]). Hundreds of functional responses have been proposed in the literature regarding all types of interactions: cooperation [9], plant-pollinators [10, 11], predation (reviewed in [4, 8]), competition [12]. Strikingly, despite a large variety of possible functional responses, Holling type II or related functional responses are most often used [4, 13], either for predator-prey [14] or mutualistic interactions [10]. There is however a general agreement that there is much room for improvement. It is for instance difficult to determine which one of alternative functional response fits the best to empirical data because of the poor statistical power of fitting different functions to data [14, 15]. Holling Type II functional response is certainly preferred not because it adequately modelizes ecological interactions but rather because it is a saturating function with a single parameter. Many authors argue that bridging the gap between interactions at the level of the individuals (the microscopic scale) and functional responses (the macroscopic level) would be a critical milestone in the field [6, 16], first because it is necessary to mechanistically justify functional responses, and second because evidence of individual traits variation in functional responses accumulate [e.g. 6, 17].

Few studies aimed at making the link between microscopic and macroscopic scales in the context of functional responses. Some authors used a mean-field approximation to derive functional response in consumer-resources relationships [e.g. 1, 18]. In particular, inspired by his famous “Disc experiment”, Holling showed that considering mechanisms such as searching times proportional to prey density and constant handling times of prey by predators, gives the well-known Holling’s Type II functional response [1]. Holling later showed that Holling Type III functional response can be derived by introducing the mechanism of predators learning [2]. Other authors used deterministic approaches derived from chemical reactions equations to show which assumptions must be verified for a specific functional response to be valid, e.g. Holling functional responses [19], Beddington-DeAngelis functional response [20], plant-marking pollinators interactions [11]. However, these approaches have strong limitations. First, they can be used only to derive simple functional responses. Indeed, approximations of deterministic functional responses can be obtained from a system of ordinary differential equations under assumptions of slow/fast processes. Approximations might generally be difficult to obtain because the number of equations and simplifying assumptions to make increase with the complexity of the system. It makes this approach unlikely to be general enough to embrace the large variety of possible ecological contexts. Second, since only the mean rates of interaction at the population level are considered, variability of traits and behaviors between individuals can not be considered. Third, since these approaches are deterministic, they are only valid in very large populations and the importance of stochasticity can not be taken into account.

Deriving stochastic models of functional responses have two purposes. First, it makes explicit the assumptions underlying the average interaction rates, i.e. the mean of functional responses. It has been performed in few studies in specific situations. [21] analyzed a Markov-chain model of a predator-prey system adapted from chemical reactions, and showed how Holling functional responses could emerge at the population level. [22] and [23] used a similar approach to derive a functional response in the specific cases of kleptoparasitism and a Beddington-DeAngelis functional response, respectively. Second, the development of stochastic models are needed for inferring functional responses from empirical data in order to clearly identify the processes and mechanisms underlying the variability of interaction rates, which appears to be large in experiments [e.g. 17]. Three sources of variability are possible: i) environmental or exogenous variability, e.g. temperature or observation errors; ii) inter-individual variation of behavior or traits, e.g. size, color or run speed; iii) endogenous variability due to stochastic fluctuations of the interaction themselves, hereafter called interaction stochasticity.

Modelling interactions processes at the microscopic level is necessary to evaluate the part of variance due to interaction stochasticity. [12] is the only study to our knowledge which proposed an explicit expression for an approximation of variance due to interactions, in the specific case of competition for resources. They assumed a site-based models where individuals are randomly distributed in patches every generation. They showed how different functional responses for competition can be derived depending on the assumptions about how resources are shared between competitors. A general stochastic framework for the derivation of functional responses is however still lacking. This would help to better justify the choice of the form given to functional responses, in interpreting data and making statistical inference, and evaluating the importance of interaction stochasticity for the dynamics and stability of populations and communities. In addition, since the variance due to interaction stochasticity comes from the interactions process itself, it is a valuable source of informations for inference, and would help for parameters estimation and models selection.

In this paper, we propose a general stochastic framework which allows accounting for interaction mechanisms at the level of individuals and the derivation of functional responses at the population or community level. It is based on the modelling of the distribution of the times separating two interactions and on the use of the so-called Renewal Theory, a well-known mathematical stochastic theory (for the reader’s convenience, we provide a brief account of this theory in Supp. Mat. A.1). It is classically used in foraging theory but has never been used, to the best of our knowledge, to derive functional responses and bridge the gap between behavioral ecology and population ecology. We first show how functional responses and their stochastic fluctuations can be approximated under a wide range of possible individual behaviors and interaction types. Second, we show how our model can be used to derive stochastic versions of classical functional responses, such as Holling functional responses. We then show how to derive stochastic functional responses in many different ecological contexts, through two examples: the rate of successful copulations for males in a nuptial-feeding species, and the feeding rate of predators when handling of the prey depends on its size. Finally, we apply our framework for inference through a model comparison framework by reanalyzing a dataset on grey partridges [24]. We show in particular that information can be extracted from observed fluctuations in order to estimate model parameters and improve inference, given that variability due to interactions themselves are large enough relatively to other sources of variations.

2 Functional responses from Renewal Theory: A general framework

Consider a community with three species denoted by $e_{x}$, $e_{y}$ and $e_{z}$. We take three species for concreteness, one can of course consider more species, and one can consider genotypes, phenotypes, substrates or resources instead of species. Our goal is to determine the number of times $N\!_{\Delta}(x,y,z)$ a focal individual of species $e_{x}$ successfully interacts with other individuals of a given species $e_{y}$ (possibly the same species), during a time span $\Delta$, with $x,y,z$ the size or density of each species in the environment.

We assume $1\ll N\!_{\Delta}(x,y,z)\ll x,y,z$, i.e. i) the number of interactions is low compared to the size of the populations during $\Delta$, and ii) the number of interactions during $\Delta$ is large. In other words, it is assumed that interactions between individuals have a negligible effect on the population sizes and the time span $\Delta$ is large enough for many interactions to occur. Under these assumptions, the time between the $k-1^{th}$ and $k^{th}$ interactions, denoted by $T_{k}(x,y,z)$, is a random variable whose distribution is independent of $k$, but generally depends on $x,y,z$. The number of interactions $N\!_{\Delta}$ during $\Delta$ is thus a random variable defined as follows (Figure 1):

[TABLE]

Approximation of the mean and variance of a functional response.

Functional responses, i.e. the expected number of successful interactions per unit of time (in the more specific case of predator-prey interactions the number of prey killed by a predator per unit of time, e.g. [1, 7, 8]), can be defined as

[TABLE]

Under the assumptions that the number of interactions is large (i.e. $\Delta$ is large) but the environment does not change during this time (i.e. variation in densities $x$, $y$ and $z$ is negligible), the mean and variance of the functional response $R(x,y,z)$ can be approximated in law according to Renewal Theory by

[TABLE]

where $\mathscr{N}(0,\sigma^{2})$ denotes a Normal distribution with mean [math] and variance $\sigma^{2}$ (see Supp. Mat. A.1 for mathematical details). Since all the random variables $T_{k}(x,y,z)$ have the same distribution, we dropped the index $k$ and simply wrote $T(x,y,z)$ inside the mean and the variance. This formula shows that knowing the expected time between two successful interactions $\mathds{E}\left[T(\cdot)\right]$ and its variance $\mathrm{Var}\left[T(\cdot)\right]$ provides the mean of a functional response and a confidence interval. It is important to notice that the stochastic fluctuations of the functional response around its mean are only due to the interactions per se and not to external sources: since we consider interactions between individuals as a stochastic process, it implies stochastic fluctuations. We will call this intrinsic source of variability interaction stochasticity in order to differentiate it with possible others, such as environmental stochasticity or between individuals variability.

In order to give an explicit form of a functional response, the next step is to decompose the time $T(x,y,z)$ according to the different events taking place between two interactions.

Decomposing the time between two interactions.

We suppose that the time between two interactions, namely $T(x,y,z)$, can be decomposed into a sequence of a given number of independent steps, in a definite order. We start with a simple example, and give the general version afterwards. Consider a predator $e_{x}$ and a prey $e_{y}$ (the third species $e_{z}$ is ignored for simplicity). Assume that a predator has two activities: searching a prey and handling it if its capture is successful (see Fig. 1). Suppose it takes a time $\tau_{h}$ to handle a prey and $\tau_{s}$ to find it. The first attempt of the predator to find a prey takes a random time $\tau_{s}^{(1)}$. The prey is then either caught (‘success’) or not (‘failure’). If the capture is successful, then the time of an interaction is $T(x,y)=\tau_{s}^{(1)}(x,y)+\tau_{h}(x,y)$. Otherwise, the predator has to make another attempt, which takes a new random time $\tau^{(2)}_{s}$. If there is a success, then $T(x,y)=\tau_{s}^{(1)}(x,y)+\tau_{s}^{(2)}(x,y)+\tau_{h}(x,y)$, and so on and so forth. The reader can guess that the number of steps is given by the number of trials needed to get the first success. It is well-known that this number is a random variable which follows a geometric law whose parameter is the probability of a success.

The decomposition of the time between two interactions can be generalized and applied to many different ecological contexts, for any number of activities. Denote by $\mathcal{A}$ the (finite) set of possible activities and by $|\mathcal{A}|$ its cardinality, i.e. the number of necessary activities in order to perform a new successful interaction. In the previous example, we had $\mathcal{A}=\{``\text{handling}",``\text{searching}"\}$ and $|\mathcal{A}|=2$. Let $\tau_{a}(x,y,z)$ be the random time needed to carry out a given activity $a\in\mathcal{A}$, e.g., $a=``\text{searching}"$. Let $p_{a}(x,y,z)$ be the probability that this activity ends up successfully, and $G_{a}(x,y,z)$ the random number of attempts the focal individual must perform before succeeding activity $a$. $G_{a}(x,y,z)$ is a random variable following a geometric distribution with parameter $p_{a}(\cdot)$. Then we can write (the notation $(x,y,z)$ is dropped hereafter not to overburden the notations, but the reader should keep in mind that all random variables generally depend on $x,y,z$):

[TABLE]

where $\tau^{\Sigma}_{a}$ is the total time spent for activity $a$ to end up successfully, given by

[TABLE]

where, the random variables $\tau^{(i)}_{a}$, $i=1,2,3,\ldots$, are independent and have the same distribution as $\tau_{a}$ (note that $G_{a}$ appears as a superscript because it is the random number of attempts before succeeding activity $a$).

We can compute the mean and variance of the time span $T$ between two interactions, and we find (see Supp. Mat. A.2)

[TABLE]

and

[TABLE]

We finally obtain an approximation of the functional response by plugging Eqs. (3) and (4) into Eq. (2). To achieve this, one needs to know $\tau_{a}(x,y,z)$ and $p_{a}(x,y,z)$ for all $a\in\mathcal{A}$, that is, the time taken by every activity and their probability of success, both generally depending on the species density in the community. These quantities can be determined empirically, or explicitly specified in various ecological contexts. We give examples and applications in the following.

3 Examples of functional responses with searching and handling times

In this section, we show how stochastic versions of classical functional responses can be derived assuming simple processes at the level of the individuals for searching and handling activities.

A general form of the functional response with searching and handling times.

Our aim is to determine the expectation and variance of the number of interactions between a focal individual $e_{x}$ and individuals $e_{y}$ during a time $\Delta$. We want to take into account possible perturbations due to interactions between the focal individual $e_{x}$ with individuals of a third species $e_{z}$. We suppose that the time between two interactions $T$ in a given environment $(x,y,z)$ can be decomposed into a searching time $\tau_{s}$ and a handling time $\tau_{h}$, that handling always succeeds ($p_{h}=1$) and catching an individual $e_{y}$ has a probability $p_{s}$ to succeed. If searching fails, the focal individual starts a new searching phase or stop searching if the total time $\Delta$ is reached (Figure 1). From Eqs. (3) and (4), we can compute the expectation and variance of the time separating two interactions between the focal $e_{x}$ individual and $e_{y}$ individuals:

[TABLE]

and

[TABLE]

which give an approximation of the stochastic functional response $R$ (Eq. (2)).

In order to give explicit expression of the mean and variance of the time between two interactions $\mathds{E}(T)$ and $\mathrm{Var}(T)$, Eqs. 5 and 6 show that we now need to specify the mean and variance of the time taken to find an individual $e_{y}$ ($\mathds{E}(\tau_{s})$ and $\mathrm{Var}(\tau_{s})$), the probability to effectively interact with an individual once found ($p_{s}$), and the time taken for handling the interaction ($\mathds{E}(\tau_{h})$ and $\mathrm{Var}(\tau_{h})$). In other words, we need to explicitly specify how individuals move into space and how interactions take place. In the following, we derive functional responses under explicit assumptions regarding foraging and handling times. We first detail a simple case with assumptions which are classically made to obtain a Holling type II functional response, i.e. constant handling time, a single type of prey, and foraging time as a function of prey densities in the environment. Second, we give a generalized model, in a $d$-dimensions space which can include many different classical functional responses.

A simple case: foraging in a 2D space with handling.

We assume a focal individual $e_{x}$ foraging in a 2D space of size $L^{2}$, where $y$ individuals $e_{y}$ are uniformly distributed on a square lattice, such as the distance between the two nearest $e_{y}$ individuals is $L/(\sqrt{y}-1)$. The location of the individual $e_{x}$ is randomly chosen in the 2D space and we want to calculate the expectation and variance of the distance $D(y)$ to the closest $e_{y}$ individual. The horizontal $u_{1}$ and vertical distances $u_{2}$ between $e_{x}$ and $e_{y}$ follow uniform distributions in the range $\left[0,\delta_{2}\right]$ with $\delta_{2}=L/2(\sqrt{y}-1)$, which gives

[TABLE]

where $C_{2}^{E}$ and $C_{2}^{V}$ are two constants (explicit values are given in Supp. Mat. A.3 ). It is assumed that the focal individual $e_{x}$ has a perfect knowledge of the spatial distribution of individuals $e_{y}$ and goes to the nearest $e_{y}$ individual following a straight line at speed $v$. Hence, the expectation and variance of the searching time are $\mathds{E}\left(\tau_{s}\right)=\mathds{E}\left(D(y)\right)/v$ and $\mathrm{Var}\left(\tau_{s}\right)=\mathrm{Var}\left(D(y)\right)/v^{2}$. Assuming that searching always succeeds ($p_{s}=1$) and that handling time is a constant $c_{y}$ ($\mathrm{Var}\left(\tau_{h}\right)=0$), we finally get from Eqs. 5 and 6 the expectation and variance of the functional response $R$

[TABLE]

with $\lambda=L/2v$, the scaled size of the environment (the size of the environment $L/2$ relative to foraging speed $v$). Predictions provided by Eqs. 7 are in very good agreement with individual-based simulations of the present model (Figure 2).

The functional response given by Eq. 7 is a saturating function of the density $y$ which looks like a Holling II functional response. Classically, the form of the Holling II functional response is justified by two mechanisms: a searching time depending on the density of the prey, and a constant handling time [1]. This is thus not surprising that we recover a functional response close to Holling II here. However, the exact form differs: we obtain a function of $\sqrt{y}$ instead of $y$. This illustrates that being explicit about how individuals forage into the environment, with specific justifications about the mechanisms underlying interactions between individuals, can give rise to alternative functional responses (other alternatives are given in the following section). This also shows that adopting a bottom-up approach allows to estimate how variable is the number of interactions in a given ecological context due to intrinsic interaction stochasticity. In the present case, the variance of the functional response decreases with $y$.

Foraging and handling one or two species in a $d$-dimensions space.

We give here a generalization of the framework with handling and foraging in $d$-dimensions. The individual $e_{x}$ now forages for two possible species $e_{y}$ or $e_{z}$, with constant handling times $c_{y}$ and $c_{z}$. We make similar assumptions than in the previous section, but we introduce the parameters $\alpha$, denoting a possible preference of the species $e_{x}$ for species $e_{z}$, and $\beta$ denoting a different availability or vulnerability of species $e_{z}$ (e.g. if $\beta>1$, $e_{z}$ is easier to be detected relatively to $e_{y}$, see Supp. Mat. A.3 for details). We will moreover compare two different movements followed by the focal $e_{x}$ individual: a straight line to the nearest individuals (as in the previous section) or a Brownian motion. We want to calculate the time between two interactions between $e_{x}$ and $e_{y}$ which depends on: i) the number of occurrence where $e_{x}$ interacts with $e_{z}$ instead of $e_{y}$, which follows a geometric distribution with probability of success $y/(y+\alpha z)$; ii) the two first moments of the time taken to reach a given species in a $d$-dimensions space $\theta^{E}_{d}(w)$ and $\theta^{V}_{d}(w)$, respectively (see Supp. Mat. A.3 for details); iii) the handling times $c_{y}$ and $c_{z}$. Finally, the time between two interactions between a focal individual $e_{x}$ and individuals $e_{y}$ are (Supp. Mat. A.3):

[TABLE]

in the case of a straight movement to the closest patch and

[TABLE]

in the case of a Brownian motion.

With appropriate assumptions regarding underlying mechanisms, Eqs. 2, 8 and 9 are general enough to recover some classical functional responses and their variance. As shown in Table 1, assuming a 1D space and direct foraging to the nearest individual yields the equations for Holling Type I, II and III, Beddington-DeAngelis or the Ratio-Dependence functional responses, given further assumptions about handling time and possible interference between species. On a side note, we were not able to find appropriate assumptions to recover the exact equation of the Holling Type III functional responses as defined in the literature: We found a function of the form $y^{2}/(1+y+y^{2})$ instead of the classical form $y^{2}/(1+y^{2})$. It can be due either to a lack of generality of Eq. 8, or because there is no realistic biological assumptions which allows to recover the same exact form under our framework, or because Holling Type III is only correct when $y$ is large, i.e. $y$ is negligible relative to $y^{2}$. Even if both functions have similar sigmoidal shapes, our results illustrate that deriving functional responses with a bottom-up approach highlights implicit and often hidden assumptions regarding the underlying mechanisms.

Interestingly, alternative underlying mechanisms can give rise to similar functional responses, in two ways. First, the exact same equation can emerge for different hypotheses. There are for instance two ways to introduce competition between predators to obtain the ratio-dependence functional response, assuming either that the probability of searching success is $p_{s}=1/x$ or that the searching time is proportional to the density of individuals $e_{x}$, $\mathds{E}(\tau_{s}(x,y,z))=\lambda\ x/y$ and the probability of success is $p_{s}=1$. Second, we can obtain similar form of the functional responses with alternative foraging strategies in a given environment. Fig. 3 shows that decelerating functional responses can be obtained without handling time ($c_{y}=0$) but with foraging in a 2D or 3D space. Sigmoidal functional responses (Holling Type III-like forms) can emerge if foraging follows a Brownian motion in a 1D space. Note that not only forms of functional responses vary with assumptions, but also their range (interactions rates have different scales even though parameter values are identical, see the $y$-axes on Fig. 3). This illustrates the relevancy of the bottom-up approach developed here. Since the same functional responses can be obtained under different hypotheses, deducing underlying mechanisms from an observed form at the macroscopic level is limited; for instance, a sigmoidal functional response does not necessarily mean that a predator is able to learn. Instead, we propose to modelize interactions at the microscopic level (the level of the individuals) in order to get corresponding functional responses.

Another merit to adopt a bottom-up approach is that the variance of the functional responses can be estimated (Tab. 1 and Fig.3). The intensity of stochastic fluctuations of the functional responses, depicted by the 95% confidence interval on Fig. 3, depends on the supposed underlying mechanisms (e.g. the large fluctuations in the Ratio-Dependence functional response are due to the probability of searching success of the form $1/x$), and on the density of species (e.g. fluctuations increase linearly with $y$ when there is no handling). Also note that even if, as shown before, we can recover the same mean functional response with different underlying hypotheses, variance can be different: the ratio-dependence functional response shows different variance whether $p_{s}=1/x$ and $\mathds{E}(\tau_{s}(x,y,z))=\lambda/y$ or $p_{s}=1$ and $\mathds{E}(\tau_{s}(x,y,z))=\lambda\ x/y$. This is noteworthy because it shows that, in order to make inference from data, both the mean and the variance of the number of interactions per unit of time can be used to discriminate between concurrent functional responses models.

4 Build the $<$insert your name$>$ functional response

In addition to recover classical functional responses, we argue that our theoretical framework is general enough to be applied to many different ecological contexts. In this section, we give two detailed examples inspired from the empirical literature. This illustrates how anyone with a specific question can derive functional responses and its stochastic fluctuations from basic knowledge about individual traits and behaviors, following a bottom-up approach.

Nuptial feeding: how many successful matings for a male?

In many species such as insects, spiders or birds, males bring food to females in order to increase their probability to mate [25, 26, 27]. Let us imagine that five successive steps must be fulfilled by a male to mate with one female: find a free a gift (e.g. a prey), handle it, find a free female, court it and copulate. Let denote $\tau_{S_{G}}$, $p_{S_{G}}$ and $\tau_{H_{G}}$ respectively the time taken and probability to successfully search for a free gift (i.e. not already handled by another male or female), and to handle it; $\tau_{S_{F}}$, $p_{S_{F}}$ and $\tau_{H_{C}}$, $p_{H_{C}}$ the time taken and the probability to successfully search for a female, and to court it; finally let $p_{R}$ be the probability to successfully copulate with it. Applying Eqs. (3) and (4) gives the expectation and the variance of the time taken for a male to successfully mate with a female:

[TABLE]

Assuming male-male and male-female competition for finding prey, the probability to successfully find a free gift by a focal male can be supposed as $p_{S_{G}}=1/(m+f)$, with $m$ and $f$ the male and female densities, respectively. The form given to $p_{S_{G}}$ is arbitrary and, assuming that all females and males have equal chance to find prey, it only reflects that the higher the density of competitor males and females, the more difficult for the focal male to find a gift. Note that the time spent by the focal male to find a free gift also depends on prey density through $\tau_{S_{G}}$. Similarly, we assumed that males compete for finding a free female, with equal chance, which gives the probability to successfully find a female to court is $p_{S_{F}}=1/m$. Handling and court times, and the probability to copulate are assumed constant (respectively equal to $C_{G}$, $C_{C}$ and $p_{R}=p$). An approximation of the expectation and variance of the number of successful matings by a male can then be obtained using Eq. 2. Figure 4 shows that the variance of the number of successful copulations decreases when the number of males and females decreases, and that the expected number of successful copulations has a non-monotonous variation with the density of females. Such an approach can be useful for studying the evolution of interaction behavior in the context of nuptial-feeding behaviors [26].

Trait-dependent interaction times: prey size and functional responses.

While most functional responses make the implicit assumption that interaction rates only depend on the number or density of individuals, many investigations suggest that interaction rates might also depend on individual traits and their distribution in a population, especially body mass or size [e.g. 28, 29]. For instance, [30] showed that functional responses generally depend on the ratio between predator and prey masses, and [15] showed that predators prefer prey with a particular body mass. We show here that our theoretical framework can be used to derive functional responses dependent on quantitative traits.

As an illustration, we will focus on the derivation of handling times accounting for body mass. We neglect searching times for the sake of simplicity, but it might be relevant to include the effect of size on the time and probability to find and catch a prey: larger prey can be more easily detected or caught. The effect of prey size on searching time can for instance be taken into account into the spatial scale $L$ introduced before. Let us note $s$ the size of the prey and $m(s)$ its distribution, supposed following a truncated normal distribution defined on $[0,+\infty]$, with mean $\mu_{s}$ and standard deviation $\sigma_{s}$. $\pi(s)$ is defined as the density probability that a predator catches a prey of size $s$. Finally, we suppose that the expected handling time conditional on the prey size $s$ is $\mathds{E}(\tau_{h}|s)=t(s)$. The expectation and variance of handling times are thus

[TABLE]

We can now be more specific in order to give explicit expression of the functional response. Let $P(s)$ be the preference of the predator for prey with size $s$, assumed to follow a truncated Gaussian function defined on $[0,+\infty]$, with $\mu_{P}$ the most preferred size and $\sigma_{P}$ its tolerance. The probability to catch a prey of size $s$ is therefore $\pi(s)=m(s)P(s)ds/\int m(u)P(u)du$. We can reasonably assume that handling time is a bounded increasing function of prey size $s$, such as $\mathds{E}(\tau_{h}|s)=\tau_{\text{max}}(1-e^{-s})$. Finally, using Eq. (2), an approximation of the functional response and its variance can be obtained (Fig. 5). Figure 5 shows how the functional response and its variance change as a function of the variance of the body size of the prey $\mu_{s}$. This illustrates that functional responses not only depend on the average body size, as shown several times [e.g. 31, 32], but also on the variance of prey body size, possibly giving non-monotonous relationships.

5 Functional responses inference: A model comparison framework

Many works compare different functional responses models with likelihood ratio test or the Akaike’s Information Criterion in order to determine which one fits the best [e.g. 14, 24, 15, 33]. However, it is generally difficult to discriminate between alternative functional responses [14, 15]. A limit of these approaches is that the variance of the functional responses is directly estimated from the data, while variance depends on interactions themselves, as shown by our framework (Eq. 2). Here, we reanalyze a dataset from [24] to estimate the feeding rate of grey partridges on seeds in a controlled experiment. [24] especially aimed at testing whether vigilant behaviors significantly affect the feeding rate. They compared two models using the AIC method: Holling type II with and without vigilant behaviors. They found no statistical difference between models with or without vigilance, and concluded that vigilant behaviors do not affect feeding rate. Here, we first derive a functional response with vigilant behaviors using our bottom-up approach. Second, we use a likelihood approach to estimate the parameters and test whether vigilance significantly affects grey partridges feeding rates (see Supp. Mat. A.4 for details).

We supposed that the time between two successful interactions, i.e. between two eaten seeds, can be decomposed into vigilance bouts with constant duration $c_{v}$, occurring with probability $p_{v}$, searching bouts $\tau_{s}$, successful with probability $p_{s}$, and handling bouts with constant duration $c_{h}$. We fitted and compared models with ($p_{v}\neq 0$) or without vigilance ($p_{v}=0$) using a maximum likelihood approach. Note that parameters of the model were estimated both from the mean and variance (since mean and variance of the functional responses are explicitly expressed as a function of the ecological parameters). The best model was estimated to be in a 2D space, with a direct movement to the nearest seed, with vigilance (Likelihood-ratio of the models with and without vigilance = 641.3). We estimated the probability of entering a vigilant bout $p_{v}\simeq 0.12$, which roughly corresponds to the estimation from direct observations in [24]. A comparison of the resulting functional responses is shown on Fig. 6.

To illustrate the importance of the information included in the variance, we also fitted the models assuming a fixed variance independent of the parameters. In this case, the model without vigilance is the best (using an AIC comparison) because it has a lower number of parameters to estimate. In practice, only using the expectation of the functional response makes difficult to distinguish between different functional responses. This is however made possible using the variance. In our case, the large variability in the functional responses are better explained by the model with vigilance because the probability to enter into a vigilance bout mechanically increases variance in our model (Eq. 2).

These results must yet be taken with caution. The large variation in the functional response shown in the data (Fig. 6) can be due to other sources of variability, observation errors or variability between individuals, not taken into account into the present analyses. Indeed, as shown in [17] for a predator fish, most of the variation in the number of prey consumed by unit of time can be due to variation in the predator individual traits: searching and handling times both decrease with predator body size. This suggests, in the case of the grey partridges treated here, that the large observed variance is not necessarily due to the probability to enter into a vigilance bout, but to between-individuals variability. In order to properly analyze such type of data and improve our ability to infer functional responses parameters, an adapted statistical framework must be developed taking into account both between- and within-individuals variability. Within-individuals variability is taken into account by our model (Eq. 2) but a mixed-effect model is needed to take into account between-individuals variability. Ideally, taking account of between prey variability would also be needed, for instance to take into account the effect of the variability in prey body size. Classically, mixed effects models include between-individuals variability only in the mean part of the model, and not in the variance. In our model, since between-individuals variability would affect both the mean and the variance of the model (Eq. 2), a specific mixed-effect statistical framework is needed (development in progress).

6 Discussion

Dealing with identifiability. The rate at which individuals interact, within or between species, is central in most ecological processes: it affects individuals’ growth, birth and death, populations’ and communities’ dynamics, and how energy and matter flow through ecosystems. Predicting ecological dynamics and supporting management decisions depend on our understanding of how individuals interact. It is thus crucial to use functional responses models that fit the best with data and with mechanisms underlying interactions [6]. However, functional response models are hardly identifiable: often, different functional responses fit well to data and inferences show no or few statistical power [14, 15]. This identifiability problems of functional responses sometimes even resulted in violent debates [7, 8]. Adopting a bottom-up approach, we show in the present paper that such an identifiability is not surprising since, on the one hand, various underlying mechanisms can give rise to similar or identical functional responses (Fig. 3). On the other hand, our results show that choosing a functional response with a phenomenological approach does not give insight about the underlying mechanisms.

Ignoring the various sources of variability of the interactions rates is certainly responsible in part of the difficulty to infer functional responses from data. It is now recognized that both interspecific and intraspecific variations are important in ecology [34, 35], especially for functional responses [6]. For instance, [17] and [16] respectively showed that functional responses depend on predators and prey body mass. Here we show that, in addition to within and between species variability, a third source of variability should also be considered: interaction stochasticity (Fig.7). Since ecological processes at the level of the individuals are stochastic, the number of interactions per unit of time is random and make functional responses randomly fluctuate. We argue that the variation in the number of interactions per unit of time can have three internal sources: within species variations for both types of interacting individuals (e.g. the predator and the prey), and interaction stochasticity (Fig. 7). Inferring functional responses from data would need to take those three sources of variability into account. Our model would allow the development of a statistical framework adapted to this problematic since it follows a bottom-up approach, from the individual to the macroscopic level. Since our model allows to relate the variance of functional responses to ecological parameters and individual traits (Eq. 2), variance in data can be used as a source of information to infer parameters (as we showed in reanalyzing [24]’s dataset). In many cases, one can expect interaction stochasticity to be the weakest source of variability in real data. As it is illustrated in our reanalysis of the grey partridges dataset, relying on it alone is certainly not relevant for model identification. Developing a statistical framework allowing to disentangle all sources of variability and to extract information from variance is a new challenge.

However, one of the main message of our paper is to question the relevancy of making inferences by using phenomenological functional responses: is it really important to know whether an hyperbolic Holling Type II functional response fits better the data than a sigmoidal Ginzburg-Arditi functional response [36]? We think it is more relevant to make the link between behavioral ecology, foraging theory, population ecology and community ecology. We argue that our model can help in doing this by i) stopping searching which adhoc functional response fits best to data, and ii) constructing one’s own functional response from one’s own specific case, using basic knowledge about species and their ecological contexts. Our model is general enough to adapt to many different situations, given it is possible to decompose interactions into a sequence of different activities.

Assumptions and time scales. Many situations in natural populations or experiments are obviously expected to depart from our model. We made some choices for the sake of simplicity, without really being necessary. For instance, we supposed that individuals move in straight line or following a Brownian motion, which can be modified accordingly to a specific ecological context. However, two assumptions are inherent to our theoretical framework, and for which departures can be generally found in nature. First, we assume that the interactions can be decomposed into a sequence of activities whose duration are random, independent and identically distributed, which is a fundamental hypothesis of Renewal Theory used here. Yet, individuals are in general able to learn or sets mark in the environment in order to improve their ability to search, or change their preference for a prey or another, which violates the independence in time [1, 37]. Individuals can also be clustered in space, which would for instance affect the time taken by a predator to find two succeeding prey [36]. Population size of each interacting species can also significantly vary during the considered time frame, because of prey depletion or blooming. Second, we assume that the number of interactions and the considered time frame are large. Yet in general, experiments and observations generally involve a few dozens or hundreds of interacting individuals during a limited time frame.

What is the impact of violating these assumptions on the resulting functional response is an open question. Performing stochastic simulations might obviously help in addressing this question. However, we think it can be addressed more generally from a theoretical point of view. Indeed, we think that it mostly relies on how do the different time scales of the ecological processes involved in the functional response relate with each other. For instance, if individual learning occurs at a slower time scale than it consumes prey, then we can expect that learning little affects the resulting functional response. Similarly, if the time taken by a predator to move from one patch of prey to another is large enough, then it should not affect much the resulting functional response in a short time frame. If prey reproduction is fast enough, then depletion should have a low impact. Some authors proposed for example to take into account prey depletion into account [38], but the approach is still phenomenological, i.e. different time scales are generally not considered. The relative importance of the time scales is central in ecology. It deserves particular attention in the case of functional responses in particular [39], and it is a natural perspective of the current work.

The last perspective of this work is about how individuals interactions translates into birth and death, in other words how functional responses are related to numerical responses [6]. Since most population and community dynamics models are systems of differential equations, functional responses and numerical responses are generally assumed to play at similar time scales (except when slow/fast processes are explicitly assumed). Some authors showed in specific cases that considering different time scales for functional responses and numerical responses [e.g. 40, 41] can dramatically affect population dynamics and stability. Finally, functional response should be considered as the result of stochastic ecological processes. We highlight in the present paper the phenomenon we called interaction stochasticity, in analogy with demographic stochasticity. It is well-known that demographic stochasticity can affect population stability [42]. How interaction stochasticity can affect population and community stability, through its effect in numerical responses, is an open question.

Author Contributions

All authors designed the theoretical framework. SB wrote the first draft of the manuscript, and all authors contributed substantially to revisions.

Acknowledgements

We are very grateful to David Baker for sharing its dataset on grey partridges. We want to thank Clotilde Lepers, Elisabeta Vergu and Arnaud Senlis for helpful discussions and Geoffroy Berthelot for providing simulations. SB benefited from the support by the “H.+C.+Y.+M.+V.+D.†+M. Moens-Sacchettini organization for the advancement of knowledge”.

Data Accessibility

Data are deposited on a public server (Dryad : doi:10.5061/dryad.c73bm2q)

Funding Statement

This work was supported by the Chair “Modélisation Mathématique et Biodiversité” of Veolia Environnement-Ecole Polytechnique-Museum National d’Histoire Naturelle-Fondation X and by the ANR ABIM (ANR-16-CE40-0001).

Figures

Appendix A Supplementary Materials of the paper “Rejuvenating Functional Responses with Renewal Theory”

by S. Billiard, V. Bansaye and J.-R. Chazottes.

A.1 Some basics in Renewal Theory

For the convenience of the reader, we state and prove the two results we are using in this paper. We refer to e.g. [43] or [44] for more details. Let $T_{1},T_{2},\ldots$ be a sequence of independent non negative random variables. We assume that the $T_{i}$’s are all distributed as a non negative random variable $T$ with distribution $F$. We assume that $F(0)\neq 1$ (i.e., the random variables are not identically zero) and $\mathds{E}(T)<+\infty$. Finally, let $S_{n}=\sum_{i=1}^{n}T_{i}$ for $n\geq 1$. The renewal process $(N_{t})_{t\geq 0}$ associated to $T_{1},T_{2},\ldots$ is defined by

[TABLE]

Since $\mathds{P}(T\geq 0)=1$ and $\mathds{P}(T=0)<1$, we have $\mathds{E}(T)>0$. By the strong law of large numbers, we know that $S_{n}/n$ converges to $\mathds{E}(T)$ with probability one. Therefore we deduce that $S_{n}\to+\infty$ with probability one. In particular, with probability one, $S_{n}<t$ for only a finite number of $n$’s, showing that $\mathds{P}(N_{t}<+\infty)=1$. It follows that, for all practical purposes, we can put a ‘max’ instead of a ‘sup’ in the definition of $N_{t}$.

Strong law of large numbers for renewal processes. We have

[TABLE]

Proof. By the very definition of $N_{t}$ we have $S_{N_{t}}\leq t<S_{N_{t}+1}$, implying that

[TABLE]

Now $N_{t}\to+\infty$ with probability one, hence

[TABLE]

To state the second result we use, we have to assume in addition that the variance of $T$, which we denote by $\sigma^{2}$, is finite.

The central limit theorem for renewal processes. We have

[TABLE]

where $\mathscr{N}(0,1)$ is the standard Gaussian distribution (centred at [math] and with variance equal to $1$).

Proof. Let $r_{t}:=\frac{t}{\mathds{E}(T)}+y\sigma\sqrt{\frac{t}{(\mathds{E}(T))^{3}}}$. If it is an integer, let $n_{t}=r_{t}$, if it is not an integer, let $n_{t}=\lceil r_{t}\rceil+1$. Then

[TABLE]

We now use the central limit theorem for the sequence $T_{1},T_{2},\ldots$, and the fact that

[TABLE]

It follows that, if $Z$ is a random variable distributed according to $\mathscr{N}(0,1)$, we have

[TABLE]

This ends the proof. Loosely speaking, $N_{t}$ is approximately Gaussian with mean $\frac{t}{\mathds{E}(T)}$ and variance $\frac{t\sigma^{2}}{(\mathds{E}(T))^{3}}$.

Application to determine the stochastic fluctuations of the functional responses.

To apply the two previous theorems in our context, set $t=\Delta$ and apply them for each given triplet $(x,y,z)$. More precisely, the functional response, defined as $R(x,y,z)=N\!_{\Delta}(x,y,z)/\Delta$, is a random variable and we want to determine its fluctuations around its expected value $1/\mathds{E}\left(T(x,y,z)\right)$. This can also be achieved by applying Renewal Theory. Colloquially, we get that the distribution of

[TABLE]

is approximately a Gaussian distribution centred at [math] with a variance equal to

[TABLE]

where \mathrm{Var}\left(T(x,y,z)\right)=\mathds{E}\left(T^{2}(x,y,z)\right)-\big{(}\mathds{E}\left(T(x,y,z)\right)\big{)}^{2} stands for the variance of $T(x,y,z)$. In particular, one can say that standard deviation of (A.1) is approximately

[TABLE]

The larger the $\Delta$, the better these approximations. Finally, we can condense what precedes by writing

[TABLE]

where $\mathscr{N}(0,\sigma^{2})$ is a random variable having a Gaussian distribution centred at 0 and with variance $\sigma^{2}$.

A.2 Some probability facts

Sum of i.i.d random variables $X_{i}$ a random number of times $Y$. Let $Y$ be a random variable assuming positive integer values such that $E(Y)<+\infty$. Let $X_{1},X_{2},\ldots$ be a sequence of independent identically distributed random variables which are also independent of $Y$, and such that $\mathds{E}(X_{1})<+\infty$. Then

[TABLE]

Furthermore, assume that $\mathrm{Var}(Y)<+\infty$ and $\mathrm{Var}(X_{1})<+\infty$. Then

[TABLE]

Expectation and variance of the time between two interactions $T$ (Eqs. 3 and 4). Assuming that the probability of success of activity $a$ is $p_{a}$, then $Y_{a}$ the number of times activity $a$ is performed before success follows a geometric distribution with expectation and variance respectively $\mathds{E}\left[Y_{a}\right]=1/p_{a}$ and $\mathrm{Var}\left[(1-p_{a})/p_{a}^{2}\right]$. Assuming that activity $a$ takes a time $X_{a}=\tau_{a}$, with expectation $\mathds{E}\left[\tau_{a}\right]$ and variance $\mathrm{Var}\left[\tau_{a}\right]$, which finally gives the expectation and variance of the time between two interactions $T$ in Eqs. 3 and 4.

Expectation and variance of the time taken to search and find $e_{y}$ (Eqs. 8 and 9). We want to calculate $T$ the time taken by a focal individual $e_{x}$ to interact with an individual $e_{y}$ given there are individuals $e_{z}$ in the environment. $Y_{s}$ is the the number of times the individual $e_{x}$ forages the environment before encountering $e_{y}$, hence $e_{x}$ interacts $Y_{s}-1$ times with $e_{z}$ before interacting with $e_{y}$. Assuming that $e_{x}$ has a probability $y/(y+\alpha z)$ to choose to searching for individual $e_{y}$, the expected number of time $e_{x}$ interacts with $e_{z}$ before interacting with $e_{y}$ is $\mathds{E}\left[Y_{s}-1\right]=(y+\alpha y)/y-1$. Further,in the case of a straight movement to the nearest prey, assuming that the expected time taken to search for individuals $e_{z}$ is $\theta^{E}_{d}(\beta z)$, that handling time of $e_{z}$ by $e_{x}$ is a constant $c_{z}$, that the expected time to search for individual $e_{y}$ is $\theta^{E}_{d}(y)$, and that handling time of $e_{y}$ by $e_{x}$ is a constant $c_{y}$ then the total expected time taken by interaction between $e_{x}$ and individuals $e_{z}$ is given by

[TABLE]

which corresponds to Eq. 8 (similar calculations gives Eq. 9 in the case of a Brownian motion).

In order to calculate the variance of the time taken by a focal individual $e_{x}$ to interact with an individual $e_{y}$, we need to calculate the second moment of $T$. Assuming that the second moment of the time taken to search for individuals $e_{z}$ is $\theta^{V}_{d}(\beta z)$ and to search for individuals $e_{y}$ it is $\theta^{V}_{d}(y)$, the second moment of $T$ is given by

[TABLE]

which after some calculations gives the variance of $T$ in Eq. 8.

A.3 Foraging in a $d$-dimension space

The individual $e_{x}$ forages an environment where two species $e_{y}$ or $e_{z}$ are present with densities $y$ and $z$, respectively. The environment is a $d$-dimensions space with size $L^{d}$. The handling times of the individual $e_{x}$ with species $e_{y}$ and $e_{z}$ are assumed constant, respectively denoted by $c_{y}$ and $c_{z}$. We assume that the focal individual $e_{x}$ arrives randomly at a given position in space, and that it has a perfect knowledge of the environment. The parameter $\alpha$ denotes a possible preference of the species $e_{x}$ (if $\alpha>1$) or rejection (if $\alpha<$) for resource $e_{z}$ relatively to species $e_{y}$. Because of this possible preference, the individual $e_{x}$ chooses to interact next with an individual of species $e_{y}$ with probability $y/(y+\alpha z)$ or with an individual of species $e_{z}$ with probability $\alpha z/(y+\alpha z)$. Given the species with which individual $e_{x}$ will interact next, say $e_{z}$, we have to compute the time taken for $e_{x}$ to reach an individual $e_{z}$. The parameter $\beta$ denotes an availability or vulnerability of species $e_{z}$ relatively to $e_{y}$ during the foraging process (for instance species $e_{z}$ is easier to detect if $\beta>1$). The expected time to reach an individual $e_{y}$ is denoted by $\theta^{E}_{d}(y)$, and the time taken to reach an individual $e_{z}$ is denoted by $\theta^{E}_{d}(\beta z)$. Two different movements are considered: following a straight line to the nearest patch or following a Brownian motion at speed $v$. Once the interaction has taken place, foraging starts again from a random position in the environment. There is thus no memory in the foraging process, which is a basic assumption of Renewal Theory used in our framework.

A.3.1 Regular repartition in a $d$-dimensions space, movement to the nearest individuals

The focal $e_{x}$ individual has a random location in the environment and moves to the nearest individual. We denote $D_{d}(y)$ the expected distance between the $e_{x}$ individual and and individual of species $e_{y}$ given there is a number $y$ individuals in the environment (similarly we denote $D_{d}(\beta z)$ the expected distance between $e_{x}$ and $e_{z}$). Assuming that $e_{y}$ and $e_{z}$ individuals have a regular distribution in a $d$-dimensions space of size $L^{d}$, the distance between two individuals $e_{y}$ is ${\delta_{d}(y)}=\frac{L}{2}\frac{1}{y^{1/d}-1}$, and the distance between two individuals $e_{z}$ is ${\delta_{d}(\beta z)}=\frac{L}{2}\frac{1}{(\beta z)^{1/d}-1}$ (accounting for $\beta$, the vulnerability or availability of $e_{z}$). The two first moments of the distance $D_{d}(w)$ are given by

[TABLE]

with $C^{E}_{d}$ and $C^{V}_{d}$ two constants depending on the dimension, given by

[TABLE]

Given that the $e_{x}$ individual moves at speed $v$, the two first moments of the time taken by the $e_{x}$ individual to encounter an $e_{w}$ individual in a $d$-dimension space are given by

[TABLE]

Finally, given the focal individual $e_{x}$ has a preference $\alpha$ for $e_{z}$ individuals, and given that the handling time of $e_{z}$ individuals by $e_{x}$ is $c_{z}$, then the total expected searching time $\tau_{s,d}(x,y,z)$ of an $e_{y}$ individual by an $e_{x}$ individual in a $d$-dimensions space is given by

[TABLE]

A.3.2 Regular repartition in a 1-dimension space, Brownian motion

We assume as before that $e_{z}$ have weights $\alpha$ and $\beta$, and that $e_{y}$ and $e_{z}$ are regularly distributed in space, with a distance $\delta_{1}(y)=L/(2(y-1)$ between individuals $y$ and $\delta_{1}(\beta z)=L/(2(\beta z-1)$ between individuals $e_{z}$. Note that we focus here on a $1$-dimension space for simplicity. Let us denote $H(l,a,b)$ the random variable representing the time taken by the $e_{x}$ individual moving following a Brownian motion to first hit the boundary of a segment $[a,b]$, given that its initial location $l$ is random in $l\in[a,b]$. The Laplace transformation of $H(l,a,b)$ is

[TABLE]

The two first moments of $H(l,a,b)$ are respectively given by $\frac{\partial\mathcal{L}(u)}{\partial u}|_{u=0}$ and $\frac{\partial^{2}\mathcal{L}(u)}{{\partial u}^{2}}|_{u=0}$ [45]. Hence, the two first moments of the time taken by an $e_{x}$ individual to reach either an individual $e_{w}$ are respectively (with $a$=0 and $b=\delta_{1}(w)$): $\theta_{\mathrm{{\scriptscriptstyle BM}}}^{E}(w)=C^{E}_{\mathrm{{\scriptscriptstyle BM}}}\,{\delta_{1}(w)}^{2}/v$ and $\theta_{\mathrm{{\scriptscriptstyle BM}}}^{V}(w)=C^{V}_{\mathrm{{\scriptscriptstyle BM}}}\,{\delta_{1}(w)}^{4}/v^{2}$ with $C^{E}_{\mathrm{{\scriptscriptstyle BM}}}=\frac{1}{6}$, $C^{V}_{\mathrm{{\scriptscriptstyle BM}}}=\frac{1}{15}$. Finally, given a weight $\alpha$ of $e_{z}$ individuals an handling time of individuals $e_{z}$ by $e_{x}$ is $c_{z}$, the expected searching time of an $e_{y}$ individual by an $e_{x}$ individual in a $1$D space with Brownian motion is given by

[TABLE]

and its variance is

[TABLE]

A.4 Analysis of Baker et al. (2010)’s data

In the [24] experiments, the behavior of grey partridges was observed and recorded in controlled conditions, with variable seeds densities. The feeding rate was calculated, i.e. the number of seeds eaten per second per individual. Other behaviors were also recorded such as the duration of vigilance bouts or handling times. The duration of observations bouts was at most of 240 seconds. [24] derived a version of Holling Type II functional response with vigilance in a phenomenological manner and fitted it with their data. They showed no statistical difference between models with or without vigilance and concluding that vigilant behaviors do not affect the feeding rates of grey partridges. Here, we derive a functional response with vigilance from the level of the individuals (following our bottom-up approach). We first give the expectation and the variance of the feeding rate, and second we fit our models on [24]’s data following a maximum likelihood approach, using both the mean and the variance of the functional response.

A.4.1 Derivation of a functional response with vigilance

We derive the functional response with vigilance by decomposing the time $T$ between two successful interactions, i.e., two eaten seeds, in three steps. First, we suppose that the time between two eaten seeds is decomposed into a time taken to approach a seed and a time to handle it, respectively denoted by $\tau_{a}$ and $\tau_{h}$. Supposing that the handling time is fixed $\mathds{E}(\tau_{h})=c_{h}$, that handling and approaching are always successful, we have

[TABLE]

which gives an approximation of the stochastic functional response as given by Eq. (2). Second, we suppose that the time taken to approach a seed is decomposed into times to successfully foraging for a seed, denoted by $\tau_{f}$ which succeeds with probability $p_{f}$. We thus have

[TABLE]

Third, we suppose that the time taken by successful foraging can be decomposed into a succession of vigilance bouts, with duration $\tau_{v}$ and searching time denoted by $\tau_{s}$, such as $\tau_{f}=\tau_{v}+\tau_{v}+\cdots+\tau_{s}$. We denote $p_{v}$ the probability that the animal enters in a vigilance bout, hence the expected number of vigilance bouts before encountering a seed is $1/(1-p_{v})-1$. Finally, assuming that the vigilance times are fixed, $\mathds{E}(\tau_{v})=c_{v}$ and $\mathrm{Var}(\tau_{v})=0$, we have

[TABLE]

The last step consists in determining mean and variance of the searching times, which depends on the movement of the individual and the dimension of space and the density of seeds $y$, as shown in the main text (see Eqs. (8) and (9)).

A.4.2 Likelihood ratio test

From the model given before, we can derive an approximation of the expectation and the variance of the functional response as a function of seeds density $y$ and the set of parameters $\Theta$: the time taken by a vigilance bout $c_{v}$ and the probability to enter in a vigilance bout $p_{v}$, the time taken by handling seeds $c_{y}$, the probability that foraging succeeds $p_{f}$. Our goal is to test whether vigilant behaviors significantly affect the functional response given a data set $\mathcal{(}D)$ from the experiments in [24]. To answer this question, a commonly used method is to estimate the parameters set $\Theta$ of two different models, without ($p_{v}=0$) or with ($p_{v}\neq 0$) vigilance bouts, and to compare their likelihood by a likelihood ratio test since both models are nested.

The likelihood of the functional response $R(\Theta,y)$ given the dataset $\mathcal{D}$ is defined as

[TABLE]

where $\mathds{P}(E)$ is the probability of $E$. An approximation of $\mathds{P}\left(\mathcal{D}|R(\Theta,y)\right)$ is given by Eq. (2) and its derivation in App. A.4.1. The dataset $\mathcal{D}$ contains the feeding rate as calculated for a given individual in eight treatments with different densities from 5 to 400 seeds.$m^{-2}$. We looked for the maximum of $\log\ell$ for different models (1D, 2D, or 3D space with movements to the nearest and Brownian motion in a 1D space), with or without vigilance, using the newton method in the “FindMaximum” procedure of Mathematica 10.1 [46]. Parameters estimation was constrained to adequate with their definition (for instance probabilities must lie between 0 and 1) and with our information about the protocol of the experiments (the time duration $\Delta$ was estimated in our model and supposed to be lower than 240 seconds; the scaled size of the environment $\lambda$ was supposed between 1 and 20; handling times and vigilance bouts were measured in seconds). The best model was the 2D space with movement to the nearest with vigilance (see main text).

To illustrate the importance of the fact that the parameters are estimated not only by fitting the expectation of the functional response but also its variance, we also fitted the models with and without vigilance using nonlinear model fitting (“NonLinearModelFit” in Mathematica 10.1 [46]), assuming a fixed variance independent of the parameters. Following this procedure, the model without vigilance had the lowest AIC and thus considered as the best model, in agreement with [24]’s results and conclusions. However, this is due only to the fact that it has one less parameter to estimate, hence the AIC chooses the simplest model.

As shown in the main text, the AIC of the model fitted with this latter procedure is lower than the model fitted with our framework (using Eq. (2)). The estimation of the parameters are more constrained when both variance and expectation depends on the ecological parameters: the fitting procedure must find a parameters set satisfying the best trade-off between both the expectation and the variance of the functional response. When using the information from the variance of the data to estimate parameters, the different models show largely different log-likelihoods, showing that the difference of the AIC both procedures is not essentially due to the difference in the number of parameters to be estimated. Our results suggest a different conclusion than the one by [24]: since the best model is with vigilance, we conclude that it significantly affects the feeding rate of grey partridges in the context of this experiment. Our results and conclusions illustrate the importance to use the information from the variance of the feeding rate, and not only its mean. Especially, our framework makes possible to express this variance as a function of the ecological parameters, because of the bottom-up approach: describing the behaviors and properties at the level of the individuals to make emerge functional responses at the macroscopic level.