The Stochastic Nature of Functional Responses

Abstract

Functional responses are non-linear functions commonly used to describe the variation in the rate of consumption of resources by a consumer. They have been widely used in both theoretical and empirical studies, but a comprehensive understanding of their parameters at different levels of description remains elusive. Here, by depicting consumers and resources as stochastic systems of interacting particles, we present a minimal set of reactions for consumer resource dynamics. We rigorously derived the corresponding system of ODEs, from which we obtained via asymptotic expansions classical 2D consumer-resource dynamics, characterized by different functional responses. We also derived functional responses by focusing on the subset of reactions describing only the feeding process. This involves fixing the total number of consumers and resources, which we call chemostatic conditions. By comparing these two ways of deriving functional responses, we showed that classical functional response parameters in effective 2D consumer-resource dynamics differ from the same parameters obtained by measuring (or deriving) functional responses for typical feeding experiments under chemostatic conditions, which points to potential errors in interpreting empirical data. We finally discuss possible generalizations of our models to systems with multiple consumers and more complex population structures, including spatial dynamics. Our stochastic approach builds on fundamental ecological processes and has natural connections to basic ecological theory.

1. Introduction

Models of predator-prey dynamics are the archetypal description of consumer-resource interaction. Since the pioneering work of Lotka and Volterra [1,2], coupled first-order ordinary differential equations to describe the temporal variation of the abundance of two species, a consumer and a resource, are used by ecologists in all sorts of applications [3,4,5,6,7,8,9]. In this context, the term functional response [10,11,12] was introduced to capture the variation in the rate of consumption of a resource by a consumer when the density of the resource changes. Functional responses measure how per capita average rates of resource consumption respond to the density of both resources and consumers. Such functions have been derived by isolating the feeding process, that is by considering that both resource and consumer densities remain constant. These responses can be summarized by plotting the number of resource items consumed per unit time by a consumer as a function of resource density. Holling [11,12] used the term type I to describe a linear relationship between feeding rates and resource density, while the term type II describes a non-linear relationship between feeding rates and resource density, where the slope of the curve monotonically decreases with increasing resource density, saturating at a constant value of resource consumption. Such a functional form is obtained assuming consumers move from a phase where they actively search for resources to a handling phase where they are occupied with consuming/digesting and not able to search for new resource items [13]. The term type III is associated with similar non-linear functions, where the slope first slightly increases, then, after a certain threshold, steeply increases, and finally, decreases, producing a sigmoidal shape in the curve [14]. It has been shown both theoretically [15] and empirically [16] that type III functional responses have a stabilizing effect on population dynamics. Assuming that consumers can also interfere with each other, e.g., by competing for the same resource item during searching and handling, the functional response becomes a function of both resource and consumer density. Non-linear functional forms in both consumer and resource density were independently proposed by Beddington [17] and DeAngelis [18]. The Beddington–DeAngelis functional response for predator interference has been further characterized and extended in subsequent studies [19,20,21,22].

Functional responses are typically measured in feeding experiments with a controlled resource density, which isolates the feeding process and enables quantifying resource intake. Several statistical methods to infer functional response parameters from feeding experiments are available [23,24,25]. Functional response parameters can also be directly inferred from consumer resource time series obtained from empirical studies of population dynamics [16,26,27,28]. Both empirical approaches and corresponding inference methods are based on the assumption that the parameterizations of functional responses obtained by feeding experiments are analogous to the parameterizations of the 2D ODEs describing consumer resource population dynamics. The type II functional response is still the most widely used functional response, and the corresponding differential equations for predator-prey dynamics, also known as the the Rosenzweig–MacArthur equations [29], have been widely studied in the dynamical systems literature [30,31,32]. Furthermore, the dynamical properties of consumer resource equations with the Beddington–DeAngelis functional response are well studied [33,34].

A complete review of the properties of these equations is beyond the scope of this paper. Instead, here, we underline how to derive such macroscopic equations from the elementary processes affecting the dynamics of individuals. Statistical physics has developed an array of powerful tools to scale up from microscopic to macroscopic dynamics [35,36]. These methods, originally developed to describe particle systems or chemical kinetics, can be more generally applied to many-body systems, and their use in population and community ecology has become widespread [37,38,39,40,41,42,43,44,45]. Examples of such derivations were given by Dawes and Souza [46], who proposed a minimal set of reactions to derive type I, II, and III functional responses, and by Van der Meer and Smallegange [21], who used similar arguments to derive a stochastic version of the Beddington–DeAngelis functional response. We provide here a simpler derivation that can be easily extended to include a more general set of reactions and accommodate different aspects of consumer resource interactions. We made the model general by allowing immigration and resource logistic growth and by defining the formation of a consumers-resource complex and triplets composed by two consumers fighting for one resource, which describes handling and interference, respectively. Although the stochastic description of the feeding process we provide here is not exhaustive, since we are not including direct consumer interference, it provides the basis for more elaborated reaction schemes. The paper is structured as follows

• In Section 2, we give the description of a minimal set of reactions to define consumer-resource dynamics, including predators handling prey and predators interfering. We use the Van Kampen, systems size expansion to get, through the master equation of the discrete system, a mean field version of the model consisting in a set of four ODEs.

• In Section 3, we provide a detailed analysis of the dynamical properties of the 4D ODE system, including bifurcation analysis.

• In Section 4, we show how different arguments of time scales’ separation give rise to the classical 2D systems described in the ecological literature, i.e., predator-prey equations with Beddington–DeAngelis and Holling type II and type I functional responses.

• In Section 5, we derive classical functional responses (Holling type II and III and Beddington–DeAngelis) inspired by [14]. We emphasize their stochastic nature, by deriving the probability distribution of feeding events. By doing so, we clarify the assumptions that enable all these derivations. We obtained functional responses that differ in their parameterization and in their functional form (for the Beddington–DeAngelis case) from the functional responses obtained via separation of time scales in the previous section.

• In Section 6, we discuss potential limitations and possible generalizations of the model, describing future avenues of research.

2. A Minimal Set of Stochastic Processes

Here, we describe the dynamics of $S=2$ species, a resource ${\varphi }_R$ and a consumer $X_A$. We characterize the dynamics distinguishing the two species by defining a maximum number of resource individuals $N\gg 0$ ($N\in \mathbb{N}$) that can be packed in the local community. This parameter plays the role of the size of the system. This characterization introduces a limit in the total immigration rate of the resource species, ${\Lambda }_R$, between zero, when the system is full with N resource individuals, and N, when it is empty. This results in a density-dependent total immigration rate and a density-independent death rate.

(1)$\begin{array}{ccc}\hfill \varnothing & \stackrel{\hspace{1em}{\lambda }_R\hspace{1em}}{\to }& {\varphi }_R,\hfill \end{array}$

(2)$\begin{array}{ccc}\hfill {\varphi }_R& \stackrel{\hspace{1em}{\delta }_R\hspace{1em}}{\to }& \varnothing ,\hfill \end{array}$

where ${\delta }_R$ is the intrinsic death rate of the resources and ${\Lambda }_R\left(t\right)={\lambda }_R(N-n_R\left(t\right))$. Note that we can always consider the system’s size as divided into N small portions or sites, which means that ${\lambda }_R$ can be interpreted as an immigration rate per empty site. Here, $n_R\left(t\right)$ is the number of resource individuals in the local community at time t. To further characterize resource dynamics, we considered another reaction describing birth events,

(3)$\begin{array}{ccc}\hfill {\varphi }_R+\varnothing & \stackrel{\hspace{1em}\beta \hspace{1em}}{\to }& {\varphi }_R+{\varphi }_R.\hfill \end{array}$

Here, the per capita resource rate is given by $B\left(t\right)=\beta \frac{N-n_R\left(t\right)}{N}$, and $\beta$ can be regarded as a per capita birth rate when the system is almost empty.

The dynamics of the consumers is described by immigration and death reactions,

(4)$\begin{array}{ccc}\hfill \varnothing & \stackrel{\hspace{1em}{\lambda }_A\hspace{1em}}{\to }& X_A,\hfill \end{array}$

(5)$\begin{array}{ccc}\hfill X_A& \stackrel{\hspace{1em}{\delta }_A\hspace{1em}}{\to }& \varnothing ,\hfill \end{array}$

and by processes associated with interactions between consumer and resource individuals describing, e.g., the feeding process with reactions:

(6)$\begin{array}{ccc}\hfill X_A+{\varphi }_R& \stackrel{\hspace{1em}\alpha \hspace{1em}}{\to }& X_{\left[AR\right]}+\varnothing ,\hfill \end{array}$

(7)$\begin{array}{ccc}\hfill X_{\left[AR\right]}+\varnothing & \stackrel{\hspace{1em}\nu \hspace{1em}}{\to }& \hfill X_A+X_A\end{array}$

where the formation of consumer-resource pairs $X_{\left[AR\right]}$ happens at a total rate $\frac{\alpha n_R{n}_A}{N}$ and the degradation of pairs into newly born consumers happens at rate $\nu n_{\left[AR\right]}$. In Reactions (6) and (7), $\alpha$ is the encounter rate between consumers and resources and $\nu$ is a degradation rate of pairs, defining the “handling time” ${\tau }_H:=1/\nu$ of individual consumers.

We can finally consider another process related to consumer interference formalized in triplet formation, given by reactions of the kind:

(8)$\begin{array}{c}\hfill X_{\left[AR\right]}+X_A\stackrel{\hspace{1em}\chi \hspace{1em}}{\to }{X}_{\left[ARA\right]}+\varnothing ,\end{array}$

(9)$\begin{array}{c}\hfill X_{\left[ARA\right]}+\varnothing \stackrel{\hspace{1em}\eta \hspace{1em}}{\to }{X}_{\left[AR\right]}+X_A,\end{array}$

where $\chi$ is the rate of triplet formation that measures the propensity of free consumers to attack handling consumers and $\eta$ is the rate of triplet degradation, which defines the “interference time” ${\tau }_I:=1/\eta$. Note that we did not consider interference between free consumers, i.e., reactions leading to the formation of consumer pairs $X_{\left[AA\right]}$.

Derivation of Mean Field Approximations

The process defined by Reactions (1)–(9) can be formulated as a stochastic process. Fixing the number of sites to N, the system state is completely defined by a vector $\mathbf{n}=(n_R,n_A,n_{\left[AR\right]},n_{\left[ARA\right]})$, and the dynamics is given by a discrete state Markov process, where $P(\mathbf{n},t|{\mathbf{n}}_{\mathbf{0}},t_0)$ denotes the probability that the system is in state $\mathbf{n}$ at time t, given it was in state ${\mathbf{n}}_{\mathbf{0}}$ at time $t_0<t$.

The transition rates associated with elemental processes are:

(10)$\begin{array}{cc}\hfill & T_R^{+}:=T(n_R+1,n_A,n_{\left[AR\right]},n_{\left[ARA\right]}|\mathbf{n})={\lambda }_R(N-n_R)+\frac{\beta n_R(N-n_R)}{N-1},\hfill \end{array}$

(11)$\begin{array}{cc}\hfill & T_R^{-}:=T(n_R-1,n_A,n_{\left[AR\right]},n_{\left[ARA\right]}|\mathbf{n})={\delta }_R{n}_R,\hfill \end{array}$

(12)$\begin{array}{cc}\hfill & T_A^{+}:=T(n_R,n_A+1,n_{\left[AR\right]},n_{\left[ARA\right]}|\mathbf{n})={\lambda }_{A}N,\hfill \end{array}$

(13)$\begin{array}{cc}\hfill & T_A^{-}:=T(n_R,n_A-1,n_{\left[AR\right]},n_{\left[ARA\right]}|\mathbf{n})={\delta }_A{n}_A,\hfill \end{array}$

(14)$\begin{array}{cc}\hfill & T_{\left[AR\right]}^{+}:=T(n_R-1,n_A-1,n_{\left[AR\right]}+1,n_{\left[ARA\right]}|\mathbf{n})=\frac{\alpha n_A{n}_R}{N},\hfill \end{array}$

(15)$\begin{array}{cc}\hfill & T_{\left[AR\right]}^{-}:=T(n_R,n_A+2,n_{\left[AR\right]}-1,n_{\left[ARA\right]}|\mathbf{n})=\nu n_{\left[AR\right]},\hfill \end{array}$

(16)$\begin{array}{cc}\hfill & T_{\left[ARA\right]}^{+}:=T(n_R,n_A-1,n_{\left[AR\right]}-1,n_{\left[ARA\right]}+1|\mathbf{n})=\frac{\chi n_A{n}_{\left[AR\right]}}{N},\hfill \end{array}$

(17)$\begin{array}{cc}\hfill & T_{\left[ARA\right]}^{-}:=T(n_R,n_A+1,n_{\left[AR\right]}+1,n_{\left[ARA\right]}-1|\mathbf{n})=\eta n_{\left[ARA\right]}.\hfill \end{array}$

These transition rates are consistent with the set of reactions introduced above.

Let $P(\mathbf{n},t)$ denote the conditional probability $P(\mathbf{n},t|{\mathbf{n}}_{\mathbf{0}},t_0)$ obtained by imposing that the initial state, at $t=t_0$, is ${\mathbf{n}}_{\mathbf{0}}$, i.e., $P(\mathbf{n},t_0)=\delta (\mathbf{n}-{\mathbf{n}}_{\mathbf{0}})$. With this definition, the master equation reads as:

(18)$\begin{array}{cc}\hfill \frac{\partial P}{\partial t}& =({\mathbb{E}}_R^{-1}-\mathbb{I})\left(T_R^{+}P\right)+({\mathbb{E}}_R-\mathbb{I})\left(T_R^{-}P\right)\hfill \\ \hfill & +({\mathbb{E}}_A^{-1}-\mathbb{I})\left(T_A^{+}P\right)+({\mathbb{E}}_A-\mathbb{I})\left(T_A^{-}P\right)\hfill \\ \hfill & +({\mathbb{E}}_R{\mathbb{E}}_A{\mathbb{E}}_{\left[AR\right]}^{-1}-\mathbb{I})\left(T_{\left[AR\right]}^{+}P\right)+({\mathbb{E}}_A^{-2}{\mathbb{E}}_{\left[AR\right]}-\mathbb{I})\left(T_{\left[AR\right]}^{-}P\right)\hfill \\ \hfill & +({\mathbb{E}}_A{\mathbb{E}}_{\left[AR\right]}{\mathbb{E}}_{\left[ARA\right]}^{-1}-\mathbb{I})\left(T_{\left[ARA\right]}^{+}P\right)+({\mathbb{E}}_A^{-1}{\mathbb{E}}_{\left[AR\right]}^{-1}{\mathbb{E}}_{\left[ARA\right]}-\mathbb{I})\left(T_{\left[ARA\right]}^{-}P\right),\hfill \end{array}$

where we introduced the one-step operators:

(19)$\begin{array}{cc}\hfill & {\mathbb{E}}_{X}f(\dots ,n_X,\dots )=f(\dots ,n_X+1,\dots ),\hfill \\ \hfill & {\mathbb{E}}_X^{-1}f(\dots ,n_X,\dots )=f(\dots ,n_X-1,\dots ),\hfill \end{array}$

to make the notation compact. Observe that operator X acts on population abundance $n_X$, with $X\in \{R,A,[AR],[ARA\left]\right\}$.

It can be shown, via a systematic Van Kampen expansion [35] of the master equation on the system’s size (N), that the stochastic model defined above yields a deterministic approximation as the leading term of the expansion, together with a Fokker–Planck equation describing the noise in the next-to-leading order [35]. The mean field version of the process associated with the general set of reactions (1)–(9) is given by a system of coupled ODEs given by:

(20)$\begin{array}{cc}\hfill \frac{d{n}_R}{dt}& ={\lambda }_R(N-n_R)-{\delta }_R{n}_R+\frac{(\beta N-\beta n_R-\alpha n_A)n_R}{N},\hfill \end{array}$

(21)$\begin{array}{cc}\hfill \frac{d{n}_A}{dt}& ={\lambda }_{A}N-{\delta }_A{n}_A+2\nu n_{\left[AR\right]}+\eta n_{\left[ARA\right]}-\frac{(\alpha n_R+\chi n_{\left[AR\right]})n_A}{N},\hfill \end{array}$

(22)$\begin{array}{cc}\hfill \frac{d{n}_{\left[AR\right]}}{dt}& =\eta n_{\left[ARA\right]}-\nu n_{\left[AR\right]}+\frac{(\alpha n_R-\chi n_{\left[AR\right]})n_A}{N},\hfill \end{array}$

(23)$\begin{array}{cc}\hfill \frac{d{n}_{\left[ARA\right]}}{dt}& =\frac{\chi n_{\left[AR\right]}{n}_A}{N}-\eta n_{\left[ARA\right]}.\hfill \end{array}$

In Appendix A, we illustrate how the Van Kampen expansion for the master equation proceeds for the case in which no triplet formation is considered, i.e., when $\chi =0=\eta$.

3. Stability Analysis

The ODE system defined by Equations (20)–(23) has a rich behavior in terms of its stability and the qualitative analysis of the equilibrium points it exhibits. We first analyzed the equilibrium points of the full ODE system and found a transcritical bifurcation when ${\lambda }_A=0$. The system has limit cycles that emerge after a Hopf bifurcation, which we show in a particular case in Section 3.1.

At equilibrium, (23) yields $\eta n_{\left[ARA\right]}=\frac{\chi n_{\left[AR\right]}{n}_A}{N}$. Substitution into the r.h.s. of (22) yields $\nu n_{\left[AR\right]}=\frac{\alpha n_R{n}_A}{N}$. Therefore, substitution of these two relations into the r.h.s. of (21) gives, at equilibrium,

(24)${\lambda }_{A}N-{\delta }_A{n}_A+\frac{\alpha n_R{n}_A}{N}=0,$

which, together with (20), forms a non-linear system of degree three for equilibrium abundances $n_A$ and $n_R$. This system can be reduced to a cubic equation for $n_R$, which reads:

(25)$\begin{array}{cc}-f_R^3+(1-{\lambda }_R^{\prime }-{\delta }_R^{\prime }+{\delta }_A^{\prime })f_R^2+\hfill & \\ & \hfill +\left[\frac{\alpha }{\beta }{\lambda }_A^{\prime }+({\delta }_R^{\prime }-1){\delta }_A^{\prime }+({\delta }_A^{\prime }+1){\lambda }_R^{\prime }\right]f_R-{\lambda }_R^{\prime }{\delta }_A^{\prime }=0,\end{array}$

where we made the definitions ${\lambda }_A^{\prime }:={\lambda }_A/\alpha$, ${\lambda }_R^{\prime }:={\lambda }_R/\beta$, ${\delta }_A^{\prime }:={\delta }_A/\alpha$, and ${\delta }_R^{\prime }:={\delta }_R/\beta$ and used the resource abundance scaled by the system’s size $f_R:=n_R/N$. However, for the sake of simplicity, in what follows, we focus on the case of the absence of the immigration of the predator, ${\lambda }_A=0$.

In this case, (24) yields a trivial solution, $n_A=0$, which implies trivially that $n_{\left[AR\right]}=0$ and $n_{\left[ARA\right]}=0$. Therefore, the r.h.s. of (20) yields a second-order equation,

(26)${\lambda }_R^{\prime }(1-f_R)-{\delta }_R^{\prime }{f}_R+(1-f_R)f_R=0.$

This equation can be solved for $f_R$,

(27)$f_R^{\pm }=\frac{1}{2}\left(q-{\lambda }_R^{\prime }\pm \sqrt{{({\lambda }_R^{\prime }-q)}^2+4{\lambda }_R^{\prime }}\right),$

where we introduced the control parameter $q:=1-{\delta }_R^{\prime }=1-\frac{{\delta }_R}{\beta }$, which can be seen as the ratio between the intrinsic growth rate ($r=\beta -{\delta }_R$) and the birth rate ($\beta$) of resources. The solution $f_R^{-}$ is always negative, so we discarded it, and we obtained the first set of solutions $n_R=N{f}_R^{+}$ and $n_A=n_{\left[AR\right]}=n_{\left[ARA\right]}=0$.

On the other hand, for ${\lambda }_A=0$, Equation (24) yields also the constant (q-independent) solution $n_R=\frac{N{\delta }_A}{\alpha }=N{\delta }_A^{\prime }$. Again, using the r.h.s. of (20), we obtained a non-trivial equilibrium solution for $n_A$,

(28)$n_A=\frac{N\beta }{\alpha }\left(-{\delta }_A^{\prime }+\frac{{\lambda }_R^{\prime }}{{\delta }_A^{\prime }}+q-{\lambda }_R^{\prime }\right)$

which, in turn, yields the equilibrium pair abundance,

(29)$n_{\left[AR\right]}=\frac{N\beta }{\nu }\left(-{{\delta }_A^{\prime }}^2+(q-{\lambda }_R^{\prime }){\delta }_A^{\prime }+{\lambda }_R^{\prime }\right),$

and triplet abundance,

(30)$n_{\left[ARA\right]}=\frac{N{\beta }^2\chi }{\eta \nu {\delta }_A^{\prime }}{\left({{\delta }_A^{\prime }}^2-(q-{\lambda }_R^{\prime }){\delta }_A^{\prime }-{\lambda }_R^{\prime }\right)}^2.$

Now, we focus on the two equilibrium abundances for the resource, $n_R^{\left(1\right)}=N{f}_R^{+}$ and $n_R^{\left(2\right)}=N{\delta }_A^{\prime }$. These two solutions cross each other at:

(31)$q^{\star }={\delta }_A^{\prime }+{\lambda }_R^{\prime }-\frac{{\lambda }_R^{\prime }}{{\delta }_A^{\prime }}=\frac{{\delta }_A}{\alpha }+\frac{{\lambda }_R}{\beta }-\frac{\alpha {\lambda }_R}{\beta {\delta }_A},$

which is basically a restriction between model parameters,

(32)$1-\frac{{\delta }_R}{\beta }=\frac{{\delta }_A}{\alpha }+\frac{{\lambda }_R}{\beta }-\frac{\alpha {\lambda }_R}{\beta {\delta }_A}.$

Then, using the Jacobian matrix of the system (20)–(23) evaluated at the two equilibria, it is easy to show that the solution for $n_R^{\left(2\right)}=N{\delta }_A^{\prime }$ is stable for $q<q^{\star }$ and unstable for $q>q^{\star }$. Consistently, the equilibrium solution for $n_R^{\left(1\right)}=N{f}_R^{+}$ is unstable for $q<q^{\star }$ and stable for $q>q^{\star }$. Therefore, we found a transcritical bifurcation at $q=q^{\star }$ for ${\lambda }_A=0$. Figure 1 shows the transcritical bifurcation as the stability of the steady-state changes by increasing the control parameter above the threshold $q^{\star }$.

3.1. Hopf Bifurcation

We also observed the emergence of a Hopf bifurcation for several ranges of model parameters. Here, we illustrate the existence of a Hopf bifurcation for ${\lambda }_R=0$, ${\lambda }_A=0$, and ${\delta }_R=0$, in which we were able to calculate a series expansion of eigenvalues for small $\beta$. In order to derive analytical expressions in the simplest setting, we considered the model in the absence of triplets.

Let us focus on the equilibrium abundances $n_R=\frac{N{\delta }_A}{\alpha }$, $n_A=\frac{N\beta }{\alpha }\left(1-\frac{{\delta }_A}{\alpha }\right)$, and $n_{\left[AR\right]}=\frac{N\beta {\delta }_A}{\alpha \nu }(1-\frac{{\delta }_A}{\alpha })$. The bifurcation arises around this equilibrium point. Observe that these densities are meaningful if $\alpha >{\delta }_A$. Notice that the threshold given by (32) reduces to $\alpha ={\delta }_A$ in the particular case we considered (${\lambda }_R=0$, ${\lambda }_A=0$, and ${\delta }_R=0$). For $\alpha >{\delta }_A$, the solution $f_R^{\left(2\right)}={\delta }_A/\alpha$ derived above is the stable one.

At that equilibrium point, the Jacobian matrix reads:

(33)$J=\left(\begin{array}{ccc}-\frac{\beta {\delta }_A}{\alpha }& -{\delta }_A& 0\\ -\beta \left(1-\frac{{\delta }_A}{\alpha }\right)& -2{\delta }_A& 2\nu \\ \beta \left(1-\frac{{\delta }_A}{\alpha }\right)& {\delta }_A& -\nu \end{array}\right).$

For $\beta =0$, it is easy to see that the eigenvalues of J are $\lambda =0$ (double) and $\lambda =-2{\delta }_A-\nu$, so the equilibrium point is neutrally stable, and a limit cycle appears. We now expand these eigenvalues for $\beta \ll 1$. The characteristic polynomial is:

(34)${\lambda }^3+\left(2{\delta }_A+\frac{\beta {\delta }_A}{\alpha }+\nu \right){\lambda }^2+\left(\frac{3{\delta }_A+\nu }{\alpha }-1\right)\beta {\delta }_A\lambda +\left(\frac{{\delta }_A}{\alpha }-1\right)\beta {\delta }_A\nu =0.$

This cubic equation can be expanded in power series of $\beta$ (for a similar expansion, we refer the reader to [47]). First, we obtained a series expansion for the double root $\lambda =0$, obtained when $\beta =0$. The perturbation analysis starts by setting $\lambda =a{\beta }^b$ and then trying to determine the factor a and exponent b at leading order. Substitution into (34) yields:

(35)$\begin{array}{cc}-a^3{\beta }^{3b}+a{\beta }^{1+b}{\delta }_A\left(1-\frac{3{\delta }_A+\nu }{\alpha }\right)+\beta \nu {\delta }_A\left(-1+\frac{{\delta }_A}{\alpha }\right)\hfill & \\ & \hfill -\frac{{\delta }_A}{\alpha }{a}^2{\beta }^{2b+1}-(2{\delta }_A+\nu )a^2{\beta }^{2b}.\end{array}$

Therefore, the smallest powers in $\beta$ are obtained by choosing $b=1/2$. Canceling the leading terms (which turn out to be of the order of $\beta$), we can determine an expression for a as the solution of:

(36)${\delta }_A\nu \left(-1+\frac{{\delta }_A}{\alpha }\right)-(2{\delta }_A+\nu )a^2=0,$

i.e.,

(37)$a=\pm \sqrt{\frac{{\delta }_A\nu ({\delta }_A-\alpha )}{\alpha (2{\delta }_A+\nu )}}.$

Here, we see that for $\alpha >{\delta }_A$, the double eigenvalue $\lambda =0$ obtained for $\beta =0$ splits into two complex eigenvalues with the imaginary part equal to:

(38)$\mathrm{Im}\left(\lambda \right)=\pm \sqrt{\frac{\beta {\delta }_A\nu (\alpha -{\delta }_A)}{\alpha (2{\delta }_A+\nu )}}.$

We can calculate the real part (at the lowest order in the $\beta$ series expansion) by setting:

(39)$\lambda =\pm \sqrt{\frac{\beta {\delta }_A\nu (\alpha -{\delta }_A)}{\alpha (2{\delta }_A+\nu )}}i+c\beta .$

Expanding the characteristic polynomial up to order ${\beta }^{3/2}$ and setting equal to zero the coefficient of the leading term yields the following expression for c:

(40)$c=-\frac{{\delta }_A({\nu }^2+2(3{\delta }_A-\alpha )({\delta }_A+\nu ))}{2\alpha {(2{\delta }_A+\nu )}^2}.$

Therefore, for $\beta \ll 1$, we obtained the two complex conjugate eigenvalues:

(41)$\lambda =-\frac{\beta {\delta }_A({\nu }^2+2(3{\delta }_A-\alpha )({\delta }_A+\nu ))}{2\alpha {(2{\delta }_A+\nu )}^2}\pm \sqrt{\frac{\beta {\delta }_A\nu (\alpha -{\delta }_A)}{\alpha (2{\delta }_A+\nu )}}i+\mathcal{O}\left({\beta }^{3/2}\right).$

By setting ${\nu }^2+2(3{\delta }_A-\alpha )({\delta }_A+\nu )=0$, we found the threshold of a Hopf bifurcation, which can be expressed as:

(42)${\alpha }_c=3{\delta }_A+\frac{{\nu }^2}{2({\delta }_A+\nu )}.$

Observe that ${\alpha }_c>{\delta }_A$, so eigenvalues stemming from the splitting of $\lambda =0$ (for $\beta =0$) are complex. For $\alpha >{\alpha }_c$, $\mathrm{Re}\left(\lambda \right)>0$, and for ${\delta }_A<\alpha <{\alpha }_c$, $\mathrm{Re}\left(\lambda \right)<0$ in the limit $\beta \ll 1$. Therefore, the eigenvalues cross over the imaginary axis, and a Hopf bifurcation arises: for $\alpha <{\alpha }_c$, eigenvalues have a negative real part, which leads to stable oscillations around the equilibrium point. Once the threshold ${\alpha }_c$ is crossed over ($\alpha >{\alpha }_c$), damped oscillations (with $\mathrm{Re}\left(\lambda \right)<0$) no longer exist, and the equilibrium point is unstable. The dynamics converges to a limit cycle (i.e., we found a supercritical Hopf bifurcation as $\alpha$ increases).

For the sake of completeness, we provide here the series expansion for the third eigenvalue, which can be obtained in a similar way:

(43)$\lambda =-2{\delta }_A-\nu -\frac{2\beta {\delta }_A(\alpha -{\delta }_A)({\delta }_A+\nu ))}{\alpha {(2{\delta }_A+\nu )}^2}+\mathcal{O}\left({\beta }^{3/2}\right).$

Notice that this eigenvalue remains negative for small perturbations in $\beta >0$.

Figure 2 shows the stability of equilibrium solutions for arbitrary values of $\alpha$ and $\beta$. The diagram shows the transcritical bifurcation threshold (32), which separates the stable solution (1), where only the resource survives, from the coexistence equilibrium (2). This solution is first asymptotically stable, and then, damped oscillations around this stable equilibrium point arise. The rightmost threshold line in Figure 2 stands for the Hopf bifurcation, which separates stable oscillations from limit cycles.

4. Derivation of Functional Responses via the Separation of Time Scales

In this section, we discuss the asymptotic dynamics of the mean field description of the system (20)–(23). To simplify the calculations, we assumed there is no immigration in the system (${\lambda }_R=0$ and ${\lambda }_A=0$). A simple asymptotic derivation can be made for the type I and type II functional responses. In this case, there is no interference between predators ($\chi =0$, $\eta =0$), and we can write the mean field Equations (20)–(23) as:

(44)$\begin{array}{cc}\hfill \frac{d{n}_R}{dt}& =r{n}_R\left(1-\frac{n_R}{K}\right)-\frac{\alpha n_R{n}_A}{N},\hfill \end{array}$

(45)$\begin{array}{cc}\hfill \frac{d{n}_A}{dt}& =2\nu n_{\left[AR\right]}-{\delta }_A{n}_A-\frac{\alpha n_R{n}_A}{N},\hfill \end{array}$

(46)$\begin{array}{cc}\hfill \frac{d{n}_{\left[AR\right]}}{dt}& =\frac{\alpha n_R{n}_A}{N}-\nu n_{\left[AR\right]},\hfill \end{array}$

where $r=\beta -{\delta }_R$ and $K=rN/\beta$ are the growth rate and carrying capacity of the resource. Assuming that compounds disappear very fast from the system, the time scale of consumers $X_A$ is longer than the time scale of the compounds $X_{\left[AR\right]}$. In this regime, putting $\frac{d{n}_{\left[AR\right]}}{dt}=0$ in Equation (46), we obtained $\nu n_{\left[AR\right]}=\alpha n_A{n}_R/N$, which we can substitute in Equation (45) to get:

(47)$\begin{array}{cc}\hfill \frac{d{n}_R}{dt}& =r{n}_R\left(1-\frac{n_R}{K}\right)-\frac{\alpha n_R{n}_P}{N},\hfill \end{array}$

(48)$\begin{array}{cc}\hfill \frac{d{n}_P}{dt}& =\frac{\alpha n_R{n}_P}{N}-{\delta }_A{n}_P,\hfill \end{array}$

which are the classical predator-prey Lotka–Volterra equations [1,2], i.e., consumer-resource equations with a type I functional response. Note that in Equation (48), the numerical response, i.e., the per capita rate of food intake of the predator [48], is equivalent to the functional response. In other words, the conversion efficiency of resource density into predator density is one. In this case, the parameter $\alpha$ can be interpreted as the macroscopic attack rate of the predator, whose density is given by $n_P=n_A$.

In the opposite limit, when the time scale of compounds $X_{\left[AR\right]}$ is longer than the time scale of consumers $X_A$, we can assume that the degradation rate of compounds and the growth rate of resources are small compared to all the other parameters. Following [46], these assumptions can be summarized by setting $\nu =ϵ\tilde{\nu }$ and $r=ϵ\tilde{r}$ with $0<ϵ\ll 1$, while $\tilde{\nu }$ and $\tilde{r}$ remain $\mathcal{O}\left(1\right)$. For any positive parameters, Equations (44)–(46) have simpler dynamics in the leading order. Making a simple rescaling in the time of the dynamics by writing $d/dt=ϵd/d\tilde{t}$ and substituting into Equation (21), we can now derive an expression for $n_A$ as:

(49)$n_A=\frac{2ϵ\tilde{\nu }{n}_{\left[AR\right]}}{{\delta }_A+\alpha n_R/N}+\mathcal{O}\left({ϵ}^2\right),$

which can be substituted into Equations (20) and (22) to obtain an ODE system of two equations given by:

(50)$\begin{array}{cc}\hfill \frac{d{n}_R}{d\tilde{t}}& =\tilde{r}{n}_R\left(1-\frac{n_R}{K}\right)-\frac{\overline{\alpha }{n}_R{n}_P}{1+\overline{\alpha }{\overline{\tau }}_H{n}_R}+\mathcal{O}\left(ϵ\right),\hfill \end{array}$

(51)$\begin{array}{cc}\hfill \frac{d{n}_P}{d\tilde{t}}& =\frac{\overline{\alpha }{n}_R{n}_P}{1+\overline{\alpha }{\overline{\tau }}_H{n}_R}-{\delta }_P{n}_P+\mathcal{O}\left(ϵ\right),\hfill \end{array}$

which is the classical Rosenzweig–MacArthur model for predator-prey dynamics [29], where the total density of predators is slaved to the total number of compounds $n_P:=n_A+n_{\left[AR\right]}=n_{\left[AR\right]}+\mathcal{O}\left(ϵ\right)$. Equations (50) and (51) describe predator-prey dynamics characterized by logistic growth for the prey and a typical type II functional response with attack rate $\overline{\alpha }=2\alpha \tilde{\nu }/{\delta }_A$, handling time ${\overline{\tau }}_H=1/2\tilde{\nu }$, and death rate ${\delta }_P=\tilde{\nu }$. In this case as well, the numerical response in Equation (51) is given by the functional response of Equation (48). Note also that in this case, the parameters of the macroscopic equations can be expressed as a function of the parameters of the microscopic process (2)–(7), but are not intuitively related to the same parameters obtained considering typical heuristic arguments to derive functional responses [25].

Similar slaving arguments can be made for the general system (20)–(23), where also the degradation rate of triplets is small compared to the other rates (i.e, we additionally set $\eta ={ϵ}^2\tilde{\eta }$). The derivation goes as follows.

Eliminating $n_A$ from Equation (21), we obtained, up to first order in $ϵ$,

(52)$n_A=\frac{2ϵ\tilde{\nu }{n}_{\left[AR\right]}}{{\delta }_A+\alpha n_R/N+\chi n_{\left[AR\right]}/N}+O\left({ϵ}^2\right),$

which can be substituted into (20) to get an equation for the resource,

(53)$\frac{d{n}_R}{d\tilde{t}}=\tilde{r}{n}_R\left(1-\frac{n_R}{K}\right)-\frac{2\alpha \tilde{\nu }{n}_R{n}_{\left[AR\right]}}{N{\delta }_A+\alpha n_R+\chi n_{\left[AR\right]}}+O\left(ϵ\right).$

On the other hand, summing Equations (22) and (23) yields:

(54)$\frac{d}{d\tilde{t}}(n_{\left[AR\right]}+n_{\left[ARA\right]})=\frac{2\alpha \tilde{\nu }{n}_R{n}_{\left[AR\right]}}{N{\delta }_A+\alpha n_R+\chi n_{\left[AR\right]}}-\tilde{\nu }{n}_{\left[AR\right]}+O\left(ϵ\right).$

If we assumed that $n_{\left[ARA\right]}$ is proportional to $n_{\left[AR\right]}$ along the dynamics, say $n_{\left[ARA\right]}=\kappa n_{\left[AR\right]}$, and we defined the total density of consumers as $n_P:=n_{\left[AR\right]}+n_{\left[ARA\right]}=(1+\kappa )n_{\left[AR\right]}$, we finally obtained the typical functional form of the Beddington–DeAngelis functional response [19], including predator handling and interference,

(55)$\begin{array}{cc}\hfill \frac{d{n}_R}{d\tilde{t}}& =\tilde{r}{n}_R\left(1-\frac{n_R}{K}\right)-\frac{\overline{\alpha }{n}_R{n}_P}{N{\delta }_A+\alpha n_R+\overline{\chi }{n}_P}+O\left(ϵ\right),\hfill \end{array}$

(56)$\begin{array}{cc}\hfill \frac{d{n}_P}{d\tilde{t}}& =\frac{\overline{\alpha }{n}_R{n}_P}{N{\delta }_A+\alpha n_R+\overline{\chi }{n}_P}-{\overline{\delta }}_P{n}_P+O\left(ϵ\right),\hfill \end{array}$

where the effective interference rate is given by $\overline{\chi }=\frac{\chi }{1+\kappa }$, the effective attack rate is $\overline{\alpha }=\frac{2\alpha \tilde{\nu }}{1+\kappa }$, and the effective death rate of predators is given by ${\overline{\delta }}_P=\frac{\tilde{\nu }}{1+\kappa }$. However, we did not find a straightforward way to determine $\kappa$ in terms of the parameters of the original dynamics given by Equations (20)–(23). In addition, we had to assume that densities $n_{\left[AR\right]}$ and $n_{\left[ARA\right]}$ are proportional along the temporal dynamics to obtain (55)–(56).

5. Stochastic Feeding Rates

A rigorous definition of a predator functional response can be precisely given as the number of resource units an average individual predator is able to consume per unit time. Since predator consumption rates respond in principle to resource density, any empirical approach intending to measure these rates will require fixing both resource and predator densities. We call these conditions chemostatic, in analogy with experiments done in chemostats, i.e., reactors where nutrient concentrations and density of microorganisms can be controlled [49]. Given the processes that define the stochastic dynamics in Section 2, we can now ask what type of functional response emerges under chemostatic conditions.

Therefore, let us assume that we maintain both the number of resource units $n_R^0$ and the total number of consumers $n_A^0=n_A+n_{\left[AR\right]}$ constant. Let us also assume first no formation of triplets and focus only on the feeding process. In this case, the dynamics can be simply described by the following two processes:

(57)$\begin{array}{ccc}\hfill X_A+{\varphi }_R& \stackrel{\hspace{1em}\alpha \hspace{1em}}{\to }& X_{\left[AR\right]}+\varnothing ,\hfill \end{array}$

(58)$\begin{array}{c}\hfill X_{\left[AR\right]}\stackrel{\hspace{1em}\nu \hspace{1em}}{\to }{X}_A,\end{array}$

where the first equation is tightly coupled to the addition of a new resource unit in order to keep resources fixed to exactly the same level $n_R^0$ all the time.

In order to calculate the functional response, i.e., the number of resource units captured by an individual predator per unit time, we introduced a new stochastic variable $R_T$, as the number of resource units consumed by a total population of consumers $n_A^0$, in a time interval, T. Then, the functional response, as a per capita feeding rate, over a certain period T and a total number of consumers $n_A^0$, becomes:

(59)$f(n_R^0,n_A^0;T)=\frac{1}{T}\frac{R_T}{n_A^0}.$

This functional response is a stochastic per capita feeding rate. In principle, it is assumed to be both prey and predator dependent. It is described by the probability distribution of cumulative feeding events n, realized by the total number of consumers $n_A^0$, between Time 0 and T. If we characterize the configuration of the system by n, the number of accumulated feeding events, and $n_A$, the number of free predators ready to attack resources at any given time, the stochastic dynamics of the feeding process is governed by the following total transition rates:

(60)$\begin{array}{cc}\hfill & r_{n,n_A}:=T\left[(n+1,n_A-1)\right|(n,n_A)]=\alpha \frac{n_R^0}{N}{n}_A,\hfill \end{array}$

(61)$\begin{array}{cc}\hfill & g_{n,n_A}:=T\left[(n,n_A+1)\right|(n,n_A)]=\nu (n_A^0-n_A).\hfill \end{array}$

Equations (60) and (61) represent the rate at which the individual-based reactions (57)–(58) occur. Compare these to reactions (6)–(7). Notice that here we did not consider consumer growth nor resource depletion. When a consumer individual attacks a resource item, this item is instantaneously replenished in order to maintain the resource density constant. These rates can be then used to write a master equation (see Equation (A14) in Appendix B), which governs the temporal evolution for the probability of having an accumulated total of n resource items being consumed and $n_A$ free consumers at time t, $P(n,n_A;t)$. This one-step process is linear in $n_A$ since both $n_A^0$ and $n_R^0$ are kept constant. Therefore, the probability distribution $P(n,n_A;t)$, can be calculated exactly. The probability distribution characterizing the stochastic variable $R_T$ is, in fact, its marginal, i.e., the probability of having a given number n of accumulated feeding events at time T, regardless how many free consumers $n_A$ are around. It will be given by:

(62)$P(n,t)=\sum _{n_A=0}^{n_A^0}P(n,n_A,t).$

Before giving more details about the calculation of both the full probability distribution $P(n,n_A;t)$ and its marginal, $P(n,t)$, we first analyzed the average feeding rate, $⟨f(n_R^0,n_A^0;T)⟩$. It is clear that the average total number of resource units $R\equiv ⟨R_T⟩$, consumed by a population of consumers monotonically increases in time according to:

(63)$\frac{dR}{dt}=\alpha \frac{n_R^0}{N}{n}_A,$

which, by defining:

(64)$\theta :=\alpha \frac{n_R^0}{N},$

can also be written in integral form as:

(65)$R\left(T\right)=\theta {\int }_0^T{n}_A\left(t\right)dt.$

Therefore, the average feeding rate per individual consumer over a period T can be written as:

(66)$⟨f(n_R^0,n_A^0;T)⟩=\frac{1}{T}\frac{\theta {\int }_0^T{n}_A\left(t\right)dt}{n_A^0}.$

Since the total number of consumers is kept constant ($n_A^0=n_A+n_{\left[AR\right]}$), the average deterministic rate equations describing the feeding processes as defined in Reactions (57) and (58) can be reduced to a single equation for $n_A$,

(67)$\frac{d{n}_A}{dt}=\nu n_A^0-(\theta +\nu )n_A.$

If we assume as initial condition $n_A\left(0\right)=n_A^0$, i.e., all consumers are free and ready to feed at time $t=0$, then we can integrate the last equation from zero to T as:

(68)$n_A\left(t\right)=\frac{\nu }{\nu +\theta }\left[1+\frac{\theta }{\nu }{e}^{-(\theta +\nu )t}\right]n_A^0,$

and therefore:

(69)${\int }_0^T{n}_A\left(t\right)dt=\frac{\nu }{\nu +\theta }\left[T+\frac{\theta }{\nu (\nu +\theta )}(1-e^{-(\theta +\nu )T})\right]n_A^0.$

Now, going back to Equation (66), we can write:

(70)$⟨f(n_R^0,n_A^0;T)⟩=\frac{\theta \nu }{\nu +\theta }+\frac{\theta }{T\nu (\nu +\theta )}(1-e^{-(\theta +\nu )T}).$

This average rate per individual consumer tends to a stationary value as T tends to infinity, which is:

(71)$⟨f\left(n_R^0\right)⟩=\frac{\theta \nu }{\nu +\theta },$

where we removed the time (to denote the asymptotic limit) and the consumer dependence. After introducing again $\theta$, as defined in Equation (78), it is clear that the average consumption rate per individual consumer corresponds exactly to the typical Holling type II functional response, where $\nu$ is the inverse of the handling time and $\alpha$ is the attack rate:

(72)$⟨f\left(n_R^0\right)⟩=\frac{\alpha \frac{n_R^0}{N}}{1+\frac{\alpha }{\nu }\frac{n_R^0}{N}}.$

See also the red curve in Figure 3. In an empirical setting, where we monitored the total number of resource units consumed by a controlled population of $n_A^0$ consumers under chemostatic conditions over a period T, we would obtain a different total number for each experimental replicate. The scheme given by the reactions (57) and (58) predicts a theoretical distribution that can then be compared to data. We give details about the calculation of this distribution in Appendix B. In short, the transition rates given by Equations (60) and (61) are first used to write a master equation from which an equation for the probability generating function can be derived:

(73)$\frac{\partial G(x,y;t)}{\partial t}=\nu n_A^0(y-1)G(x,y;t)+\left[\theta (x-y)-\nu y(y-1)\right]\frac{\partial G(x,y;t)}{\partial y},$

where the precise definition of $G(x,y,t)$ is given in Equation (58) in Appendix B. This PDE is then solved under the initial condition $G(x,y,0)=y^{n_A^0}$ by the method of the characteristics, with the normalization condition $G(0,0,t)=1$, which yields the functional form:

(74)$G(x,y;t)=e^{\nu n_A^0(y-1)t}{\left(\frac{y_0^{+}\left(x\right)e^{-\Delta \left(x\right)t}-y_0^{-}\left(x\right)\frac{y-y_0^{+}\left(x\right)}{y-y_0^{-}\left(x\right)}}{e^{-\Delta \left(x\right)t}-\frac{y-y_0^{+}\left(x\right)}{y-y_0^{-}\left(x\right)}}\right)}^{n_A^0},$

where $y_0^{-}\left(x\right)$, $y_0^{+}\left(x\right)$, and $\Delta \left(x\right)$ are functions of x (see Appendix B). It can be checked that this expression satisfies $G(0,0,t)=1$, for any t, and $G(x,y;0)=y^{n_A^0}$ as required. In addition, this allows the immediate calculation of the probability generating function of the marginal distribution $P(n;t)$ by simply evaluating the full probability generating function at $y=1$:

(75)$G(x,1;t)={\left(\frac{y_0^{+}\left(x\right)e^{-\Delta \left(x\right)t}-y_0^{-}\left(x\right)\frac{1-y_0^{+}\left(x\right)}{1-y_0^{-}\left(x\right)}}{e^{-\Delta \left(x\right)t}-\frac{1-y_0^{+}\left(x\right)}{1-y_0^{-}\left(x\right)}}\right)}^{n_A^0}.$

The distribution $P(n;t)$ gives the probability of n prey items disappearing under the pressure of a constant total number of consumers $n_A^0$, from Time 0 until time t. Under controlled experimental conditions [24], this probability can be used as the exact likelihood function for inference purposes. See also [50] and the references therein, for a broader description of the challenges of the experiments and the inference of functional response parameters. Figure 3 shows stochastic simulations of the per capita rate compared to the Holling type II average provided by (72).

Note that the Holling type II functional response in Equation (72) naturally arises as a per capita average asymptotic feeding rate after reaching stationarity under chemostatic conditions. If the feeding mechanism differs from the simple one assumed by the reactions (57) and (58), other functional responses will be obtained. For instance, if consumer feeding rates are strongly affected by the density of the resources, we can describe feeding dynamics by the following two processes:

(76)$\begin{array}{ccc}\hfill X_A+n{\varphi }_R& \stackrel{\hspace{1em}\alpha \hspace{1em}}{\to }& X_{\left[AR\right]}+(n-1){\varphi }_R,\hfill \end{array}$

(77)$\begin{array}{c}\hfill X_{\left[AR\right]}\stackrel{\hspace{1em}\nu \hspace{1em}}{\to }{X}_A.\end{array}$

If we repeat an analogous derivation as for Holling type II, we arrived formally at the same Equation (67), but instead, parameter $\theta$ should now be defined as:

(78)$\theta :=\alpha {\left(\frac{n_R^0}{N}\right)}^n$

which leads to a Holling type III functional response:

(79)$⟨f\left(n_R^0\right)⟩=\frac{\alpha {\left(\frac{n_R^0}{N}\right)}^n}{1+\frac{\alpha }{\nu }{\left(\frac{n_R^0}{N}\right)}^n}.$

Note that Reaction (76) can be regarded as if the presence of $n-1$ extra resources facilitated the attack by consumers. Resources would have a self-catalytic effect on their own loss. In nature, it is plausible that when resources are clumped together, consumers encounter them better and faster. However, a Holling type III of exact order n (Equation (79) represents here of course an idealization.

Furthermore, let us now assume consumer interference through the formation of triplets. Then, feeding dynamics is specified by the following scheme:

(80)$\begin{array}{ccc}\hfill X_A+{\varphi }_R& \stackrel{\hspace{1em}\alpha \hspace{1em}}{\to }& X_{\left[AR\right]}+\varnothing ,\hfill \end{array}$

(81)$\begin{array}{c}\hfill X_{\left[AR\right]}\stackrel{\hspace{1em}\nu \hspace{1em}}{\to }{X}_A.\end{array}$

(82)$\begin{array}{c}\hfill X_{\left[AR\right]}+X_A\stackrel{\hspace{1em}\chi \hspace{1em}}{\to }{X}_{\left[ARA\right]}+\varnothing ,\end{array}$

(83)$\begin{array}{c}\hfill X_{\left[ARA\right]}+\varnothing \stackrel{\hspace{1em}\eta \hspace{1em}}{\to }{X}_{\left[AR\right]}+X_A,\end{array}$

Again, under chemostatic conditions, the resource level $n_R^0$ and the total number of consumers $n_A^0$ are both kept constant. In order to calculate the functional response as a per capita feeding rate over certain period T, $⟨f(n_R^0,n_A^0;T)⟩$, we made use again of the definition given by Equation (66). The total consumer population is kept constant ($n_A^0=n_A+n_{\left[AR\right]}+2n_{\left[ARA\right]}$) and is distributed among the three types: free consumers $n_A$, handling/feeding consumers $n_{\left[AR\right]}$, and interfering consumers engaged in triplet formation $n_{\left[ARA\right]}$. Once this distribution reaches the steady state, there will be a steady number of consumers $n_A^{\star }$, free to attack and deplete resources. Therefore, at stationarity, Equation (66) becomes:

(84)$⟨f(n_R^0,n_A^0)⟩=\theta \frac{n_A^{\star }}{n_A^0},$

where we dropped the dependence of the averaging time interval T, because this is an asymptotic rate and $\theta$, defined again as in Equation (78), is a constant parameter, since $n_R^0$ is assumed constant. The value of $n_A^{\star }$ is defined by the steady state that emerges from the rate equations describing the feeding mechanism depicted in the reaction scheme (80)–(83). One can check that this system should read as:

(85)$\begin{array}{cc}\hfill \frac{d{n}_A}{dt}& =\nu n_{\left[AR\right]}+\eta n_{\left[ARA\right]}-\alpha \frac{n_R^0}{N}{n}_A-\chi \frac{n_{\left[AR\right]}}{N}{n}_A,\hfill \end{array}$

(86)$\begin{array}{cc}\hfill \frac{d{n}_{\left[AR\right]}}{dt}& =\eta n_{\left[ARA\right]}-\nu n_{\left[AR\right]}+\alpha \frac{n_R^0}{N}{n}_A-\chi \frac{n_{\left[AR\right]}}{N}{n}_A,\hfill \end{array}$

(87)$\begin{array}{cc}\hfill \frac{d{n}_{\left[ARA\right]}}{dt}& =\chi \frac{n_{\left[AR\right]}}{N}{n}_A-\eta n_{\left[ARA\right]}.\hfill \end{array}$

In order to determine the steady state of these feeding dynamics, we better work with densities $f_A=n_A/N$, $f_{\left[AR\right]}=n_{\left[AR\right]}/N$, and $f_{\left[ARA\right]}=n_{\left[ARA\right]}/N$ and make use of the constraint of a constant consumer population ($n_A^0=n_A+n_{\left[AR\right]}+2n_{\left[ARA\right]}$), which reduces the ODE system to the following two equations:

(88)$\begin{array}{cc}\hfill \frac{d{f}_A}{dt}& =\nu f_{\left[AR\right]}+\frac{\eta }{2}(f_A^0-f_A-f_{\left[AR\right]})-\theta f_A-\chi f_{\left[AR\right]}{f}_A,\hfill \end{array}$

(89)$\begin{array}{cc}\hfill \frac{d{f}_{\left[AR\right]}}{dt}& =\frac{\eta }{2}(f_A^0-f_A-f_{\left[AR\right]})-\nu f_{\left[AR\right]}+\theta f_A-\chi f_{\left[AR\right]}{f}_A.\hfill \end{array}$

By making these two equations equal to zero, we obtained steady state densities. In particular, the density of free consumers $f_A^{\star }$ can be written as:

(90)$f_A^{\star }=\frac{1}{4}\left[-\frac{\eta }{\chi }\left(1+\frac{\nu }{\theta }\right)+\sqrt{{\left(\frac{\eta }{\chi }\left(1+\frac{\nu }{\theta }\right)\right)}^2+8\frac{\eta }{\chi }\frac{\nu }{\theta }{f}_A^0}\right].$

Using this expression in Equation (84) and rearranging, we finally obtained the following per capita feeding rate of consumers at steady state given by:

(91)$⟨f(n_R^0,n_A^0)⟩=\frac{2\alpha \frac{n_R^0}{N}}{1+\frac{\alpha }{\nu }\frac{n_R^0}{N}+\sqrt{{\left(1+\frac{\alpha }{\nu }\frac{n_R^0}{N}\right)}^2+8\frac{\chi }{\eta }\frac{\alpha }{\nu }\frac{n_R^0}{N}\frac{n_A^0}{N}}}.$

Note that this is a predator-dependent functional response that generalizes Holling type II. If no triplet formation is considered ($\chi =0$), we recovered Equation (72). Here, we considered consumers interfering with other consumers that have already caught a resource item, perhaps in the hope of getting a share. However, other possible forms of predator interference are possible. Figure 4 shows the comparison with stochastic simulations.

The functional response given by (91) reduces almost to a Beddington–DeAngelis form in the limit $\frac{\chi }{\eta }\ll 1$ and $\frac{\alpha }{\nu }\ll 1$. In this case, (91) reduces to:

(92)$⟨f(n_R^0,n_A^0)⟩\approx \frac{\alpha \frac{n_R^0}{N}}{1+\frac{2\alpha }{\nu }\frac{n_R^0}{N}+\frac{2\chi \alpha }{\eta \nu }\frac{n_R^0}{N}\frac{n_R^A}{N}}.$

If we introduce the dimensionless rates ${\tau }_{\left[AR\right]}:=\frac{\alpha }{\nu }$, which can be interpreted as the average attack rate in units of handling time, and ${\tau }_{\left[ARA\right]}:=\frac{\chi }{\eta }$, standing for the ratio between the pace of triplet formation and triplet degradation, respectively, we found that only in the limit of ${\tau }_{\left[ARA\right]}\ll 1$ and ${\tau }_{\left[AR\right]}\ll 1$, Equation (92) is a suitable approximation of Equation (91). The equation obtained is a predator-dependent functional response of the Beddington–DeAngelis type except for the product $\frac{n_R^0}{N}\frac{n_R^A}{N}$ in the denominator. Therefore, we note that only if the system is at a very high resource density ($n_R^0\approx N$), we would tend to observe the exact Beddington–DeAngelis functional form.

6. Discussion

The aim of statistical physics is to develop phenomenological macroscopic results from a probabilistic examination of the underlying microscopic processes. This allows connecting processes at different levels of coarse-grained description. We applied this approach to consumer-resource dynamics, connecting elemental processes of growth and consumption at the individual level to a more coarse-grained description that leads, in the limit of large populations, to a deterministic description in terms of the two-dimensional typical predator-prey model (the macroscopic law, sensu Van Kampen [35]).

We set up individual dynamics in a way that makes contact with classic models in ecology. For instance, if we consider only resource dynamics, we have a population that invades an area characterized by N sites with a given immigration rate per site, ${\lambda }_R$. Then, local growth occurs, that is individuals send a number of propagules per unit time that only settle down successfully in proportion to the number of free, still non-colonized sites. This means that the system can only pack a maximum number of individuals, N. Finally, individuals die at a density-independent per capita rate, ${\delta }_R$. Under the usual assumption of one individual per site, this model corresponds to the open Levins model [51,52,53,54]. The extension of this model to S resources that compete for the N sites has been also carefully analyzed in the literature [55,56], both using one-step stochastic process and in the deterministic limit. Furthermore, if we drop local growth, then we have a system of S species undergoing an independent immigration-birth-death process, which, under ecological equivalence, leads to a typical neutral model for biodiversity [57,58]. Interestingly, when we sample a total of $N_S$ individuals from such a neutral community dynamics, which are necessarily distributed in a vector of abundances $\overrightarrow{n}=(n_1,\dots ,n_S)$, the probability of obtaining a given configuration vector, $\overrightarrow{n}$, follows exactly the Etienne likelihood function [59,60], the corner-stone of neutral biodiversity theory.

Consumer-resource dynamics was also set up in a general way although it clearly allows for further generalizations. First, we considered also that consumers immigrate from outside the N-site system at a per site rate, ${\lambda }_A$. While the system can only pack N resource individuals, in principle, there is no limit for the number of consumers. This number will be dynamically controlled by the resource level. The size of the system only directly limits resources. This differs from most approaches (e.g., [37,46,61,62]), but as we showed, it is not a problem to expand the system in terms of the size parameter N and recover, in this way, the corresponding deterministic rate equations in the large N limit. Some of the results we covered here, such as the existence of a Hopf bifurcations, can only be recovered when there is no external immigration of consumers. Since natural systems are most of the time open and subject to external inputs, this means that nice and stable consumer-resource dynamics (limit cycles) should be very difficult to observe in nature due to the stabilizing role of external immigration [63]. This does not preclude recovering stable cycles under experimental, close-to-immigration, controlled conditions [64,65].

The way we set up the subset of reactions controlling feeding dynamics has also the potential to clarify much of the discussion about predator-dependent functional responses of the Beddington–DeAngelis type [20,66,67,68]. We defined interference between consumers through structures that require multiple individuals to interact and stay together for an amount of time. Under this assumption, incidentally, we showed that, by considering only interference between handing consumers and free consumers, the exact Beddington–DeAngelis functional form cannot be fully recovered. We believe this point requires careful analysis that we will develop in a further publication.

The minimal set of reactions that we conceived reproduces consumer-resource dynamics in terms of four ODEs, Equations (20)–(23), in the deterministic limit. This involves the description of the different processes driving dynamics at the individual level. Our main result in this work has to do with the comparison between two alternative ways of deriving functional responses. First, an asymptotic approximation for the mean field version of the complete set of reactions, borrowing slaving arguments from [46], gives rise to a consumer-resource model, which predicts a functional form for the average per capita rate at which consumers deplete resources. Secondly, in a rather classic way [14,19], we derived the functional responses by describing the setting of a typical feeding experiment, that is the same system in chemostatic conditions, by only considering the subset of reactions characterizing the feeding interaction. In the first case, when going from the four-equation to the two-equation system, the dynamics of the usual kind, i.e., predator-prey equations with Holling type I (Equations (47) and (48)), type II (Equations (50) and (51)), and Beddington–DeAngelis (Equations (55) and (56)), are recovered by carefully assuming separation of time scales in the processes involved. Moreover, depending on the particular assumptions behind the separation of time scales, different functional forms may arise. In general, our analysis showed that predator feeding rates in the 2D consumer-resource models do not exactly match the typical functional responses presented in the literature, which we mostly recovered here under chemostatic conditions. Therefore, our results clearly demonstrated a discrepancy between the functional responses that emerge from the analysis of the full consumer-resource system through asymptotic derivations and separation of time scales and the functional responses derived under chemostatic conditions. For example, the attack rate and the handling time of type II functional response obtained under chemostatic conditions differ from the same parameters of the type II functional response that results from asymptotic derivations.

Although finding effective low-dimensional representations from high-dimensional dynamical systems is not an easy task, which formally requires methods from singular perturbation theory [69,70], we believe that the origins of the discrepancy we report is not related to the separation of time scales we assumed in order to go from the full, four-equation ODE system to the typical predator-prey 2D model. We rather think that it has to do with the simple way in which consumers transform resources into new consumers, the so-called numerical response of consumers in the ecological literature [48]. The instantaneous coupling between consumption and reproduction we assumed here, which is also assumed in simple predator-prey models, is not realistic and might be at the basis of this inconsistency. Structured population models, which assume that consumers accumulate mass or energy as they feed and reproduce at a later stage [71], may shed new light on the relation between individual feeding dynamics and macroscopic functional responses. Other possible generalizations, such as the explicit consideration of spatial dynamics [67] (see also Appendix C), may also help explain the relation between individual feeding dynamics and coarse-grained descriptions leading to effective functional responses. Future avenues of research include developing more robust methods to make inference from experiments of consumer-resource systems based on stochastic dynamics [26,72,73] rather than on simple ODE systems.

In sum, we demonstrated that the use of simple predator-prey models to analyze both experiments and natural predator-prey interactions should be done with caution. This includes food-web theory and applications. Models are simple representations of reality. As such, it is not rare to find inconsistencies when we try different types of models to analyze the same system. The art comes when theory is able to clearly establish the range of situations in which a certain mathematical model matches reality, which allows for both predictions and further experiments. However, a successful model of this kind should never claim that the processes considered occur in exactly the same way as they do in nature. Most of the time, process representation is a strong idealization. In this sense, we believe that ecological theories, and probably theories in general, are always effective theories.

Author Contributions

All authors contributed equally to the conceptualization, methodology, software, formal analysis, and writing. All authors read and agreed to the published version of the manuscript.

Funding

This work was funded by the Spanish Ministerio de Ciencia, Inovacion y Universidades and the the European Regional Development Fund (ERDF), under the project CRISIS (PGC2018-096577-B-I00 to D.A.).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The code to reproduce simulated data and figures can be found at https://github.com/Gpalam/The_Stochastic_Nature_of_Functional_Responses.

Acknowledgments

G.M.P. acknowledges the members of the D&M journal club at the Department of Systems Analysis, Integrated Assessment and Modelling (SIAM) of the Swiss Federal Institute of Aquatic Science and Technology (Eawag) for stimulating discussions. We also acknowledge the Theoretical and Computational Ecology group at the Center for Advanced Studies of Blanes (CEAB-CSIC) for providing a nice and supportive working environment. We thank an anonymous reviewer for useful comments that helped improve the clarity of our work.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Van Kampen Expansion of the Master Equation

In the absence of triplets, the master Equation (18) reduces to:

(A1)

$\begin{array}{cc}\hfill \frac{\partial P}{\partial t}& =({\mathbb{E}}_R^{-1}-\mathbb{I})\left(T_R^{+}P\right)+({\mathbb{E}}_R-\mathbb{I})\left(T_R^{-}P\right)+({\mathbb{E}}_A^{-1}-\mathbb{I})\left(T_A^{+}P\right)\hfill \\ \hfill & +({\mathbb{E}}_A-\mathbb{I})\left(T_A^{-}P\right)+({\mathbb{E}}_R{\mathbb{E}}_A{\mathbb{E}}_{AR}^{-1}-\mathbb{I})\left(T_{AR}^{+}P\right)\hfill \\ \hfill & +({\mathbb{E}}_A^{-2}{\mathbb{E}}_{AR}-\mathbb{I})\left(T_{AR}^{-}P\right).\hfill \end{array}$

In this Appendix, we denote $\mathbf{n}=(n_R,n_A,n_{\left[AR\right]})$. The expansion proceeds by decomposing each abundance as a term associated with a macroscopic density plus a term associated with noise,

(A2)

$\begin{array}{cc}\hfill & n_R=N{f}_R+N^{1/2}{\xi }_R,\hfill \\ \hfill & n_A=N{f}_A+N^{1/2}{\xi }_A,\hfill \\ \hfill & n_{\left[AR\right]}=N{f}_{\left[AR\right]}+N^{1/2}{\xi }_{\left[AR\right]},\hfill \end{array}$

where $\mathbf{\xi }=({\xi }_R,{\xi }_A,{\xi }_{\left[AR\right]})$ stands for the vector of noise terms. The expansion is based on the change of variables $\mathit{n}\to \mathbf{\xi }$. In the new variables, we denote $\Pi (\mathbf{\xi },t)=P(\mathit{n},t)$. Here, we use the shorthand X to denote each one of the labels in the set $\mathcal{X}:=\{R,A,[AR\left]\right\}$.

Observe that ${\mathbb{E}}_X$ changes $n_X$ into $n_X+1=N{f}_X+N^{1/2}{\xi }_X+N^{1/2}{N}^{-1/2}$. Therefore, after applying ${\mathbb{E}}_X$, the noise ${\xi }_X$ changes into ${\xi }_X+N^{-1/2}$. This yields the following asymptotic expansions for the operators ${\mathbb{E}}_X^{\pm }$ in the limit of large N,

(A3)

${\mathbb{E}}_X^{\pm }=\mathbb{I}\pm N^{-1/2}\frac{\partial }{\partial {\xi }_X}+\frac{1}{2}{N}^{-1}\frac{{\partial }^2}{\partial {\xi }_X^2}+\dots$

The time derivative in (A1) is taken at constant $\mathit{n}$. Therefore, taking time derivatives in (A2), we obtain:

(A4)

$\frac{d{\xi }_X}{dt}=-N^{1/2}\frac{d{f}_X}{dt}.$

Now, we apply the change of variables $\mathit{n}\to \mathbf{\xi }$. The chain rule yields:

(A5)

$\frac{\partial P}{\partial t}=\frac{\partial \Pi }{\partial t}+\sum _{X\in \mathcal{X}}\frac{\partial \Pi }{\partial {\xi }_X}\frac{d{\xi }_X}{dt}=\frac{\partial \Pi }{\partial t}-N^{1/2}\sum _{X\in \mathcal{X}}\frac{\partial \Pi }{\partial {\xi }_X}\frac{d{f}_X}{dt}.$

The l.h.s. of (A1) contains powers of $N^{1/2}$ and $N^0$. We expand the r.h.s. of (A1) up to the leading and next-to-leading terms in N. We proceed term by term. The first term, $({\mathbb{E}}_R^{-1}-\mathbb{I})\left(T_R^{+}P\right)$, reduces to:

(A6)

$\begin{array}{c}({\mathbb{E}}_R^{-1}-\mathbb{I})\left(T_R^{+}P\right)\approx \left(-N^{-1/2}\frac{\partial }{\partial {\xi }_R}+\frac{1}{2}{N}^{-1}\frac{{\partial }^2}{\partial {\xi }_R^2}\right)(N(1-f_R)-N^{1/2}{\xi }_R)\times \hfill \\ \hfill \left({\lambda }_R+\frac{\beta }{N-1}(N{f}_R+N^{1/2}{\xi }_R)\right)\Pi (\mathbf{\xi },t).\end{array}$

Up to order $N^0$, this yields:

$\begin{array}{cc}\hfill ({\mathbb{E}}_R^{-1}-\mathbb{I})\left(T_R^{+}P\right)& \approx N^{1/2}\left(-{\lambda }_R(1-f_R)-\beta f_R(1-f_R)\right)\frac{\partial \Pi }{\partial {\xi }_R}\hfill \\ \hfill & -\frac{\partial }{\partial {\xi }_R}\left(-{\lambda }_R{\xi }_R-\beta f_R{\xi }_R+\beta (1-f_R){\xi }_R\right)\Pi \hfill \\ \hfill & +\frac{1}{2}({\lambda }_R(1-f_R)+\beta f_R(1-f_R))\frac{{\partial }^2\Pi }{\partial {\xi }_R^2}+\mathcal{O}\left(N^{-1/2}\right).\hfill \end{array}$

Expanding up to order $N^0$ the following terms of the master equation, we obtain:

$\begin{array}{ccc}& & \begin{array}{cc}\hfill ({\mathbb{E}}_R-\mathbb{I})\left(T_R^{-}P\right)& \approx N^{1/2}\left({\delta }_R{f}_R\right)\frac{\partial \Pi }{\partial {\xi }_R}\hfill \\ & +\frac{\partial }{\partial {\xi }_R}({\delta }_R{\xi }_R\Pi )+\frac{1}{2}\left({\delta }_R{f}_R\right)\frac{{\partial }^2\Pi }{\partial {\xi }_R^2}+\mathcal{O}\left(N^{-1/2}\right),\hfill \end{array}\hfill \\ & & \begin{array}{cc}\hfill ({\mathbb{E}}_A^{-1}-\mathbb{I})\left(T_A^{+}P\right)& \approx N^{1/2}{\lambda }_A\frac{\partial \Pi }{\partial {\xi }_A}+\frac{1}{2}{\lambda }_A\frac{{\partial }^2\Pi }{\partial {\xi }_A^2}+\mathcal{O}\left(N^{-1/2}\right),\hfill \end{array}\hfill \\ & & \begin{array}{cc}\hfill ({\mathbb{E}}_A-\mathbb{I})\left(T_A^{-}P\right)& \approx N^{1/2}{\delta }_A{f}_A\frac{\partial \Pi }{\partial {\xi }_A}+\frac{\partial }{\partial {\xi }_A}({\delta }_A{\xi }_A\Pi )\hfill \\ & +\frac{1}{2}{\delta }_A{f}_A\frac{{\partial }^2\Pi }{\partial {\xi }_A^2}+\mathcal{O}\left(N^{-1/2}\right).\hfill \end{array}\hfill \end{array}$

The remaining terms are a bit convoluted,

$\begin{array}{ccc}& & \begin{array}{cc}\hfill ({\mathbb{E}}_R{\mathbb{E}}_A{\mathbb{E}}_{AR}^{-1}-\mathbb{I})\left(T_{AR}^{+}P\right)& \approx N^{1/2}\left(\alpha f_R{f}_A\right)\left(\frac{\partial }{\partial {\xi }_R}+\frac{\partial }{\partial {\xi }_A}-\frac{\partial }{\partial {\xi }_{\left[AR\right]}}\right)\Pi \hfill \\ & +\left(\frac{\partial }{\partial {\xi }_R}+\frac{\partial }{\partial {\xi }_A}-\frac{\partial }{\partial {\xi }_{\left[AR\right]}}\right)(\alpha f_A{\xi }_R+\alpha f_R{\xi }_A)\Pi \hfill \\ & +\frac{1}{2}\left(\alpha f_R{f}_A\right)\left(\frac{{\partial }^2\Pi }{\partial {\xi }_R^2}+\frac{{\partial }^2\Pi }{\partial {\xi }_A^2}+\frac{{\partial }^2\Pi }{\partial {\xi }_{\left[AR\right]}^2}\right)\hfill \\ & +\left(\alpha f_R{f}_A\right)\left(\frac{{\partial }^2\Pi }{\partial {\xi }_R\partial {\xi }_A}-\frac{{\partial }^2\Pi }{\partial {\xi }_R\partial {\xi }_{\left[AR\right]}}-\frac{{\partial }^2\Pi }{\partial {\xi }_A\partial {\xi }_{\left[AR\right]}}\right)\hfill \\ & +\mathcal{O}\left(N^{-1/2}\right),\hfill \end{array}\hfill \\ & & \begin{array}{cc}\hfill ({\mathbb{E}}_A^{-2}{\mathbb{E}}_{AR}-\mathbb{I})\left(T_{AR}^{-}P\right)& \approx N^{1/2}\left(\nu f_{\left[AR\right]}\right)\left(-2\frac{\partial }{\partial {\xi }_A}+\frac{\partial }{\partial {\xi }_{\left[AR\right]}}\right)\Pi \hfill \\ & +\left(-2\frac{\partial }{\partial {\xi }_A}+\frac{\partial }{\partial {\xi }_{\left[AR\right]}}\right)(\nu {\xi }_{\left[AR\right]}\Pi )\hfill \\ & +\frac{1}{2}\nu f_{\left[AR\right]}\left(4\frac{{\partial }^2\Pi }{\partial {\xi }_A^2}+\frac{{\partial }^2\Pi }{\partial {\xi }_{\left[AR\right]}^2}-4\frac{{\partial }^2\Pi }{\partial {\xi }_A\partial {\xi }_{\left[AR\right]}}\right)\hfill \\ & +\mathcal{O}\left(N^{-1/2}\right).\hfill \end{array}\hfill \end{array}$

Leading terms that multiply $\frac{\partial \Pi }{\partial {\xi }_R}$, $\frac{\partial \Pi }{\partial {\xi }_A}$ and $\frac{\partial \Pi }{\partial {\xi }_{\left[AR\right]}}$, respectively, yield a differential equation, which together form the following macroscopic, mean field system of ODEs:

(A7)

$\begin{array}{cc}\hfill -\frac{d{f}_R}{dt}& =-{\lambda }_R(1-f_R)-\beta f_R(1-f_R)+{\delta }_R{f}_R+\alpha f_R{f}_A,\hfill \\ \hfill -\frac{d{f}_A}{dt}& =-{\lambda }_A+{\delta }_A{f}_A+\alpha f_R{f}_A-2\nu f_{\left[AR\right]},\hfill \\ \hfill -\frac{d{f}_{\left[AR\right]}}{dt}& =-\alpha f_R{f}_A+\nu f_{\left[AR\right]},\hfill \end{array}$

which coincides with (20)–(23) for $\chi =0=\eta$ and changing back to the original variables $n_X=N{f}_X$ ($X\in \{R,A,[AR\left]\right\}$) in the limit of large N (i.e., ignoring the noise).

Next-to-leading terms yield the Fokker–Plank equation for the noise distribution $\Pi (\mathbf{\xi },t)$,

(A8)

$\frac{\partial \Pi }{\partial t}=-\sum _{i\in \mathcal{X}}\frac{\partial }{\partial {\xi }_i}({\mu }_i\Pi )+\frac{1}{2}\sum _{i,j\in \mathcal{X}}\frac{{\partial }^2}{\partial {\xi }_i\partial {\xi }_j}(D_{i,j}\Pi ),$

with a drift vector $\mathbf{\mu }$ whose components are given by:

(A9)

$\begin{array}{cc}\hfill {\mu }_R& =-({\lambda }_R+{\delta }_R-\beta +\alpha f_A){\xi }_R-\alpha f_R{\xi }_A,\hfill \\ \hfill {\mu }_A& =-\alpha f_A{\xi }_R-({\delta }_A+\alpha f_R){\xi }_A+2\nu {\xi }_{\left[AR\right]},\hfill \\ \hfill {\mu }_{\left[AR\right]}& =\alpha f_A{\xi }_R+\alpha f_R{\xi }_A-\nu {\xi }_{\left[AR\right]},\hfill \end{array}$

and diffusion matrix D whose entries are defined as:

(A10)

$\begin{array}{cc}\hfill D_{R,R}& ={\lambda }_R(1-f_R)+\beta f_R(1-f_R)+{\delta }_R{f}_R+\alpha f_R{f}_A,\hfill \\ \hfill D_{A,A}& ={\lambda }_A+{\delta }_A{f}_A+\alpha f_R{f}_A+4\nu f_{\left[AR\right]},\hfill \\ \hfill D_{\left[AR\right],\left[AR\right]}& =\alpha f_R{f}_A+\nu f_{\left[AR\right]},\hfill \\ \hfill D_{R,A}& =\alpha f_R{f}_A,\hfill \\ \hfill D_{R,\left[AR\right]}& =-\alpha f_R{f}_A,\hfill \\ \hfill D_{A,\left[AR\right]}& =-\alpha f_R{f}_A-2\nu f_{\left[AR\right]}.\hfill \end{array}$

Appendix B. The Method of Characteristics

Here, we give details about the derivation of probability generating function associated with the Holling type II consumer-resource dynamics. We characterized the configuration of the system by n, the number of accumulated feeding events, and the number of free predators, $n_A$, ready to attack resources, at any given time. We were interested in the probability distribution of cumulative feeding events, n, realized by the total number of consumers, $n_A^0$, between Time 0 and t. We calculated this distribution under the assumption that, initially, at $t=0$, there were a total of $n_A^0$ free consumers and, by definition, no feeding events yet. If we write this probability distribution as $P(n,n_A;t)$, the initial condition should be specified as:

(A11)

$P(n,m,0)=\left\{\begin{array}{cc}1& ifn=0,m=n_A^0,\hfill \\ 0& \mathrm{Otherwise}.\hfill \end{array}\right.$

According to the feeding process, the stochastic dynamics is governed by the following total probability rates:

(A12)

$\begin{array}{cc}\hfill & r_{n,m}:=T\left[(n+1,m-1)\right|(n,m)]=\theta m,\hfill \end{array}$

(A13)

$\begin{array}{cc}\hfill & g_{n,m}:=T\left[(n,m+1)\right|(n,m)]=\nu (m_0-m),\hfill \end{array}$

where, for simplicity, we defined $m:=n_A$, $m_0:=n_A^0$. Recall also that $\theta =\alpha n_R^0/N$, as usual.

This leads to the master equation, which should be always seen as a system of coupled ODEs, in our case, for $m=0,\dots ,m_0$, and n, in principle, taking integer values from −∞ to ∞:

(A14)

$\begin{array}{c}\frac{dP(n,m;t)}{dt}=\theta (m+1)P(n-1,m+1;t)+\nu (m_0-m+1)P(n,m-1;t)\hfill \\ \hfill -\left[(\theta -\nu )m+\nu m_0\right]P(n,m;t).\end{array}$

By definition, the probability generating function is:

(A15)

$G(x,y;t)=\sum _{n=\infty }^{\infty }\sum _{m=0}^{m_0}P(n,m;t)x^n{y}^m.$

Multiplying each equation in the ODE system (A14) by the corresponding factor $x^n{y}^m$ and summing over n and m, we obtained an equation in partial derivatives for the probability generating function,

(A16)

$\frac{\partial G(x,y;t)}{\partial t}=\nu m_0(y-1)G(x,y;t)+\left[\theta (x-y)-\nu y(y-1)\right]\frac{\partial G(x,y;t)}{\partial y}.$

The initial condition can be written in terms of the probability generating function in a very compact form, $G(x,y,0)=y^{m_0}$. Our objective was to find the function that fulfills both Equation (A16) and this initial condition.

We first note that x plays the role of a parameter. The relevant variables are y and t. Therefore, we looked for an integral surface $U_x(y,t)\equiv G(x,y,t)$ satisfying Equation (A16). Notice that our solution $z=U_x(t,y)$ can be regarded as a surface embedded in the 3D space defined by the variables t, y, and z. This method emphasizes that any surface in a 3D space can be characterized by a family of curves. These would be curves that, as a real constant changes, cover and define the whole surface. These curves are called characteristic curves. An equivalent way of writing our integral solution would be $Z(t,y,z)=0$ with the function $Z(t,y,z):=U_x(t,y)-z$. Furthermore, we could write the initial PDE in a very compact form, as the following dot product:

(A17)

$\overrightarrow{v}·\overrightarrow{n}=0,$

where $\overrightarrow{v}$ is defined as:

(A18)

$\overrightarrow{v}:=(1,\theta [x-y]-\nu y[y-1],\nu m_0[y-1]G(x,y;t))$

and $\overrightarrow{n}$ is the gradient of $Z(t,y,z)$, that is $\overrightarrow{n}=(Z_t,Z_y,Z_z)$, where:

(A19)

$Z_t:=\frac{\partial Z}{\partial t}=\frac{\partial G(x,y;t)}{\partial t},Z_y:=\frac{\partial Z}{\partial y}=\frac{\partial G(x,y;t)}{\partial y},Z_z:=\frac{\partial Z}{\partial z}=-1.$

Since $\overrightarrow{n}$ is the gradient of Z, it is normal to the integral surface at each and every point, $(t,y,z)$. Therefore, any curve on the surface should have a tangential direction, defined by the vector $\overrightarrow{v}$, that should be perpendicular to $\overrightarrow{n}$ at every point. Now, we write a curve in parametric form:

$\begin{array}{cc}\hfill t& =t\left(s\right),\hfill \\ \hfill y& =y\left(s\right),\hfill \\ \hfill z& =z\left(s\right),\hfill \end{array}$

and its tangent vector, $\overrightarrow{a}$, is therefore written as:

$\begin{array}{cc}\hfill a_t& =\frac{dt}{ds},\hfill \\ \hfill a_y& =\frac{dy}{ds},\hfill \\ \hfill a_z& =\frac{dz}{ds}.\hfill \end{array}$

As said, this curve is a characteristic curve of our integral solution only if its tangent vector $\overrightarrow{a}$ is parallel to $\overrightarrow{v}$, as defined in Equation (A18). Therefore, a relation of proportionality component a component should be fulfilled, i.e.,

(A20)

$\frac{\frac{dt}{ds}}{1}=\frac{\frac{dy}{ds}}{\theta [x-y]-\nu y[y-1]}=\frac{\frac{dz}{ds}}{\nu m_0[y-1]G(x,y;t)}$

which, for convenience, can be written as the following two ODEs to solve. Notice that $dz=d{U}_x=dG$:

(A21)

$\begin{array}{cc}\hfill \frac{dt}{1}& =\frac{dy}{\theta [x-y]-\nu y[y-1]}\hfill \end{array}$

(A22)

$\begin{array}{cc}\hfill \frac{dt}{1}& =\frac{dG}{\nu m_0[y-1]G}\hfill \end{array}$

After some algebra, respectively, the solutions of the first and second ODE become:

(A23)

$\begin{array}{cc}\hfill e^{-\Delta \left(x\right)t}\frac{y-y_0^{-}\left(x\right)}{y-y_0^{+}\left(x\right)}& =C_1,\hfill \end{array}$

(A24)

$\begin{array}{cc}\hfill \frac{e^{\nu m_0(y-1)t}}{G(x,y;t)}& =C_2,\hfill \end{array}$

where $C_1$ and $C_2$ are integration constants, and:

$\begin{array}{cc}\hfill \Delta \left(x\right)& =\sqrt{{(\theta -\nu )}^2+4\theta \nu x},\hfill \\ \hfill y_0^{+}\left(x\right)& =\frac{1}{2\nu }\left[-(\theta -\nu )+\Delta \left(x\right)\right],\hfill \\ \hfill y_0^{-}\left(x\right)& =\frac{1}{2\nu }\left[-(\theta -\nu )-\Delta \left(x\right)\right].\hfill \end{array}$

Recall that Equations (A23) and (A24) define a family of characteristic curves on our integral surface. This means that these curves draw this surface as a real parameter continuously changes, which implies that $C_1$ and $C_2$ cannot be freely chosen, but should depend on each other. This dependence can be written in terms of an arbitrary function, $C_1=\Psi \left(C_2\right)$. In fact, one can check that the function:

(A25)

$e^{-\Delta \left(x\right)t}\frac{y-y_0^{-}\left(x\right)}{y-y_0^{+}\left(x\right)}=\Psi \left(\frac{e^{\nu m_0(y-1)t}}{G(x,y;t)}\right)$

is a general solution of the initial PDE given by Equation (A16) whatever the function $\Psi \left(w\right)$ is. It is the initial condition that fully determines this function. We require that $G(x,y;0)=y^{m_0}$:

(A26)

$\frac{y-y_0^{-}\left(x\right)}{y-y_0^{+}\left(x\right)}=\Psi \left(\frac{1}{y^{m_0}}\right).$

One can check that this function should be:

(A27)

$\Psi \left(w\right)=\left(1-\frac{y_0^{-}\left(x\right)}{y_0^{+}\left(x\right)}\right)\frac{1}{1-w^{m_0}{y}_0^{+}\left(x\right)}+\frac{y_0^{-}\left(x\right)}{y_0^{+}\left(x\right)}.$

Plugging this expression into the general solution Equation (A25), we obtain:

(A28)

$e^{-\Delta \left(x\right)t}\frac{y-y_0^{-}\left(x\right)}{y-y_0^{+}\left(x\right)}=\left(1-\frac{y_0^{-}\left(x\right)}{y_0^{+}\left(x\right)}\right)\frac{1}{1-{\left(\frac{e^{\nu m_0(y-1)t}}{G(x,y;t)}\right)}^{m_0}{y}_0^{+}\left(x\right)}+\frac{y_0^{-}\left(x\right)}{y_0^{+}\left(x\right)}.$

After bringing $G(x,y;t)$ to the left-hand side of the equation, we obtain:

(A29)

$G(x,y;t)=e^{\nu m_0(y-1)t}{\left(\frac{y_0^{+}\left(x\right)e^{-\Delta \left(x\right)t}-y_0^{-}\left(x\right)\frac{y-y_0^{+}\left(x\right)}{y-y_0^{-}\left(x\right)}}{e^{-\Delta \left(x\right)t}-\frac{y-y_0^{+}\left(x\right)}{y-y_0^{-}\left(x\right)}}\right)}^{m_0},$

which is Equation (74) of the main text. One can check that this function is a particular solution of Equation (A16) satisfying the initial condition $G(x,y;0)=y^m$ and the normalization requirement $G(1,1,t)=1$ for all t.

Appendix C. Generalization to Spatial Dynamics

In this Appendix, we describe the dynamics of discrete individuals of a given species X moving in a two-dimensional discrete spatial lattice. The number of individuals in each site of the lattice $(x,y)$ at time t is given by $n(x,y;t)$ with $n\ge 0,n\in \mathbb{N}$, $x\in [1,\cdots ,L]$, $y\in [1,\cdots ,L]$, and $t\in \mathbb{R}$, where the total number of nodes is $M\equiv L\times L$.

The process is completely defined by the lattice edge $L>0$, $L\in \mathbb{N}$, the initial number of individuals in each site $n(x,y;0)$, ($x=1,\dots ,L$, $y=1,\dots ,L$), the possible movements, and the per capita movement rate of each individual into neighboring sites across the lattice. Assuming the same rate $\mu$ for every movement, the total rate of the process at time t is given by:

(A30)

$R\left(t\right)=\mu \sum _{i=1}^{M}n(x_i,y_i;t),$

This rate defines the distribution of times between movement events. This Markovian process, in the limit of continuum densities, can be described by the following coupled ODE system:

(A31)

$\frac{dN(x_i,y_i;t)}{dt}=-{\widehat{\mu }}_{i}N(x_i,y_i;t)+\mu \sum _{j\in \mathcal{N}(x_i,y_i)}N(x_j,y_j;t),$

for $i=1,\dots ,M$, where, here, we define:

$\begin{array}{cc}\hfill {\widehat{\mu }}_i& =\sum _{j\in \mathcal{N}(x_i,y_i)}\mu ,\hfill \\ \hfill N(x_i,y_i;t)& =\frac{n(x_i,y_i;t)}{A},\hfill \end{array}$

for $i=1,\dots ,M$, and A would be a physical area associated with every patch, which, assuming von Neumann neighbors, can be simply written as:

(A32)

$\frac{dN(x_i,y_i;t)}{dt}=4\mu \left(\frac{1}{4}\left[\sum _{j\in \mathcal{N}(x_i,y_i)}N(x_j,y_j;t)\right]-N(x_i,y_i;t)\right)$

for $i=1,\dots ,M$. Moreover, the definition of $\mu$ as a per capita random movement rate between patches ensures that the discrete model so defined also converges to the diffusion equation in the continuum limit, in the limit of the lattice step tending to zero.

In addition, the model can also be regarded as a metapopulation model where populations are located at every node of a network. General network structures can be implemented without much trouble. For instance, assume a network of patches at given physical distances, $d_{ij}$, from each other. By using the Hanski convention, that is that migration decays exponentially as between-patch Euclidean distance increases, then the per capita movement rate between patch i and j is defined as:

${\mu }_{ij}={\mu }_0{e}^{-\lambda d_{ij}}.$

In this case, the total rate given by Equation (A30) is now written as:

(A33)

$R\left(t\right)=\sum _{i=1}^M\widehat{\mu }(x_i,y_i)n(x_i,y_i;t),$

where $\widehat{\mu }(x_i,y_i)$ is defined, as in the case for ${\widehat{\mu }}_i$, analogously, as:

(A34)

$\widehat{\mu }(x_i,y_i)=\sum _{j\in \mathcal{N}(x_i,y_i)}{\mu }_{ij}.$

Throughout the main text, we described processes occurring in a single site, i.e., without spatially explicit dynamics. This is usually called the reaction part of the dynamics. However, in this Appendix, we showed that the diffusion part can be easily added by prescribing a metapopulation network topology in which individuals move from node to node at a certain per capita movement rate, where each node represents a local community.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures

Figure 1. Transcritical bifurcation for λR=2β, λA=0, and δA=3α/2. Here, the vertical axis stands for the scaled resource abundance fR, and the horizontal axis represents the control parameter q=1−δR/β. For q<q⋆=13/6, the (upper) solution fR(2)=3/2 is stable and for q>q⋆ stability changes and fR(1)=12q−2+(2−q)2+8 becomes the stable resource scaled abundance. Observe that for q>q⋆, consumer, pair, and triplet abundances become equal to zero in the stable branch. Stable (unstable) solutions are marked with full (dashed) lines.

Figure 2. Parameter space showing transitions between different stability regimes for λR=0, λA=0, and δR=0. In this case, the transcritical bifurcation threshold (32) reduces to the vertical line α=δA. The threshold αc (cf. Equation (42)) for the Hopf bifurcation is the limit (for β→0) of the line that separates stable oscillations and limit cycles. The remaining lines were computed by numerical evaluation of the Jacobian eigenvalues.

Figure 3. Plot of the per capita feeding rate of the stochastic process defined by the reactions (57)–(58). Simulations were carried out using the Gillespie algorithm with 10,000 steps, while the rates (black dots) were calculated, after a transient of 5000 steps, for different values of resource concentrations and compared with the functional response given by Equation (72) (red line). The simulation parameters are β=1.5, δA=1, α=2.5, ν=1, and N= 10,000. The total number of consumers is fixed to nA0=nA+n[AR]=200.

Figure 4. Plot of the per capita feeding rate of the stochastic process defined by the reactions (80)–(83). Simulations were carried out using the Gillespie algorithm with 15,000 steps, while the rates (black dots) were calculated, after a transient of 7500 steps, for different values of consumers concentrations, and compared with the functional response given by Equation (91) (red line). The simulation parameters are β=1.5, δA=1, α=2.5, ν=1, η=1, χ=100, and N= 10,000. The total number of resources is fixed to nR0=5000.