Interactions between different predator-prey states. A method for the derivation of the functional and numerical response

TL;DR

This paper presents a formal method to derive predator functional and numerical responses from fast state transition models, linking individual behavior to population-level responses.

Contribution

It introduces a novel systematic approach to derive functional and numerical responses based on fast state transitions, enhancing interpretability.

Findings

Explicit derivation of functional response from state transitions

Linking individual behavior to population dynamics

Method applicable to predator-prey systems

Abstract

In this paper we introduce a formal method for the derivation of a predator's functional response from a system of fast state transitions of the prey or predator on a time scale during which the total prey and predator densities remain constant. Such derivation permits an explicit interpretation of the structure and parameters of the functional response in terms of individual behaviour. The same method is also used here to derive the corresponding numerical response of the predator as well as of the prey.

Equations

$$\left\{\begin{array}[]{l}\frac{dX}{dt}=g(X)X-f(X)Y\\ \frac{dY}{dt}=\gamma f(X)Y-\delta Y\end{array}\right.$$

$$\left\{\begin{array}[]{l}\frac{dX}{dt}=g(X,Y)X-f(X,Y)Y\\ \frac{dY}{dt}=\gamma(X,Y)f(X,Y)Y-\delta(X,Y)Y\end{array}\right.$$

$$\left\{\begin{array}[]{l}\frac{dx_{k}}{dt}(t)=\sum_{i=1}^{m}A_{ki}x_{i}(t)+\sum_{i=1}^{m}\left(\sum_{j=1}^{n}B_{ij}^{(k)}y_{j}(t)\right)x_{i}(t),\quad k=1,...,m\\ \frac{dy_{k}}{dt}(t)=\sum_{i=1}^{n}\left(\sum_{j=1}^{m}C_{ij}^{(k)}x_{j}(t)\right)y_{i}(t)+\sum_{i=1}^{n}D_{ki}y_{i}(t),\quad k=1,...,n\end{array}\right.$$

$$-A_{ii}$$

$$-B_{ij}^{(i)}$$

$$-D_{ii}$$

$$-C_{ij}^{(i)}$$

$$f(X,Y)=\frac{{\sum }_{i=1}^m{\sum }_{j=1}^n{\beta }_{ij}{\hat{x}}_i{\hat{y}}_j}{Y}$$

$$g(X,Y)=\frac{{\sum }_{i=1}^m{\lambda }_i{\hat{x}}_i}{X}-\frac{{\sum }_{i=1}^m{\mu }_i{\hat{x}}_i}{X}$$

$$\gamma (X,Y)f(X,Y)-\delta (X,Y)=\frac{{\sum }_{i=1}^m{\sum }_{j=1}^n{\gamma }_{ij}{\beta }_{ij}{\hat{x}}_i{\hat{y}}_j}{Y}-\frac{{\sum }_{j=1}^n{\delta }_j{\hat{y}}_j}{Y}$$

$$\left\{\begin{array}[]{l}\bf\dot{x}=(A+B(y))x\\ \bf\dot{y}=(C(x)+D)y\end{array}\right.$$

$$\left\{\begin{array}[]{l}\bf\dot{x}=Ax\\ \bf\dot{y}=Dy\end{array}\right.$$

$$\left\{\begin{array}[]{l}\frac{dx_{k}}{dt}=\sum_{i=1}^{m}A_{ki}x_{i}(t),\quad k=1,...,m\\ \frac{dy_{k}}{dt}=\sum_{i=1}^{n}\left(\sum_{j=1}^{m}C_{ij}^{(k)}x_{j}(t)\right)y_{i}(t)+\sum_{i=1}^{n}D_{ki}y_{i}(t),\quad k=1,...,n\end{array}\right.$$

$${\bf A}=\left[\begin{array}[]{cccccc}-A_{11}&A_{12}&0&0&\dots&0\\ A_{21}&-A_{22}&A_{23}&0&\dots&0\\ 0&A_{32}&-A_{33}&\dots&\dots&0\\ \vdots&&&\ddots&&\vdots\\ 0&\dots&\dots&\dots&-A_{m-1,m-1}&A_{m-1,m}\\ 0&\dots&\dots&\dots&A_{m,m-1}&-A_{m,m}\end{array}\right]$$

$$x_{k-1}$$

$$\frac{d{x}_k}{dt}=A_{k,k-1}{x}_{k-1}-A_{kk}{x}_k+A_{k,k+1}{x}_{k+1}$$

$$\frac{d{x}_1}{dt}$$

$$\frac{d{x}_m}{dt}$$

$${\hat{x}}_k={\hat{x}}_1\prod _{i=1}^{k-1}\frac{A_{i+1,i}}{A_{i,i+1}}.$$

$${\hat{x}}_1=\frac{X}{1+{\sum }_{k=2}^m{\prod }_{i=1}^{k-1}\frac{A_{i+1,i}}{A_{i,i+1}}}.$$

$${\hat{x}}_1=\frac{X}{1+{\sum }_{k=2}^m{\mathcal{A}}_k{Y}^{k-1}}$$

$${\hat{x}}_k={\hat{x}}_1{\mathcal{A}}_k{Y}^{k-1}.$$

$$\frac{dS}{dt}=-\left(\sum _{i=1}^m{c}_i{\hat{x}}_i\right)S+dH$$

$$\hat{S}=\frac{1}{1+\frac{1}{d}{\sum }_{i=1}^m{c}_i{\hat{x}}_i}Y$$

$$\hat{H}=\frac{\frac{1}{d}{\sum }_{i=1}^m{c}_i{x}_i}{1+\frac{1}{d}{\sum }_{i=1}^m{c}_i{\hat{x}}_i}Y.$$

$$f(X,Y)=X\cdot \frac{c_1+{\sum }_{k=2}^m{c}_k{\mathcal{A}}_k{Y}^{k-1}}{1+{\sum }_{k=2}^m{\mathcal{A}}_k{Y}^{k-1}+\frac{1}{d}X\left(c_1+{\sum }_{k=2}^m{c}_k{\mathcal{A}}_k{Y}^{k-1}\right)}.$$

$$f_1\left(X,Y\right)=\frac{c_{1}X}{1+\frac{1}{d}{c}_{1}X}.$$

$$f(X,Y)=\frac{aX}{1+bX+cY}.$$

$${\hat{x}}_1=\frac{A_{12}X}{A_{12}+A_{1}Y},{\hat{x}}_2=\frac{A_{1}YX}{A_{12}+A_{1}Y}$$

$$f_2(X,Y)=\frac{c_1{A}_{12}X+c_2{A}_{1}XY}{1+\frac{A_1}{A_{12}}Y+c_1\frac{1}{d}X+\frac{c_2{A}_1}{A_{12}}\frac{1}{d}XY}.$$

$$S$$

$$H$$

Full text

K****eywords Predator-prey model  $\cdot$ Functional response  $\cdot$ Numerical response  $\cdot$ Mechanistic modelling  $\cdot$ Structured population

1 Introduction

The functional response is defined as the average number of prey captured per individual predator per unit of time as a function of the population densities of the prey, the predator or both. Well known examples are the Holling I, II and III [18] and the Beddington-DeAngelis [10, 5] functional responses. While the Holling type II functional response was derived using a time budget argument, the Holling type III as well as the functional response by Beddington and DeAngelis were introduced without an explicitly modelled underlying mechanism. However, as we show here, these responses (and many more) can be derived from a system of fast state transitions of the prey or predator during which the total prey and predator densities remain constant.

For example, Metz and Diekmann [25] derived the Holling type II functional response assuming two predator states, searching and handling, where the transition from the searching to the handling state is the result of an actual prey capture. As a consequence, the equilibrium distribution of predator densities over the two states depends on the prey density: the higher the prey density, the greater the proportion of the individual predators in the handling state and, since it is only the searching predators that capture the prey, the average number of prey captured per predator per unit of time varies with the prey density exactly as described by Holling.

Likewise, the Beddington-DeAngelis functional response, whose traditional interpretation is in terms of predator interference, was derived in a different context by Geritz and Gyllenberg [15]. They assumed two prey states, exposed and hiding, in addition to the two predator states of searching and handling. The transition from the exposed to the hiding state is mediated by the encounter with the predator. The equilibrium density distributions of both the prey and the predator over their respective states, therefore, depend on one another’s population density. As searching predators capture only exposed prey, the functional response now is not just a function of the prey density, as in the Holling type II functional response, but also of the density of the predator itself.

Equally important as the functional response are the numerical responses of the prey and the predator. One distinguishes between demographic and aggregative numerical responses. The latter is a consequence of individuals moving in space and will not be considered here. The demographic numerical response is the rate of change in population density due to birth and death as a function of the population densities of the prey, the predator or both. The same individual-level processes that determine the functional response can also determine the numerical response. Examples of how the numerical response of the prey can depend on the density of the predator (other than through prey capture) have been given by Geritz and Gyllenberg [16, 17].

Most predator-prey models in the literature are special cases of the model by Gause [13, 14]

[TABLE]

where $X$ and $Y$ are the densities of the prey and the predator, respectively, $g(X)$ is the per capita growth rate of the prey population if the predator is absent, $f(X)$ is the predator functional response, $\gamma>0$ is defined as the conversion factor of prey into predator offspring and $\delta>0$ is the per capita natural mortality rate of the predator.

The predator’s numerical response (through birth) in equation (1) is given by the term $\gamma f(X)$. However, the linear relationship between the predator’s numerical response and its functional response is lost if different prey-handling states have different fertility levels. In this paper we give an example with two predator states (starving and well-fed) where the former hunts to survive (but does not reproduce) while the latter hunts to reproduce. The proportion of predator individuals in each state depends on the prey density such that at low prey densities there are relatively more starving predators. This leads to a nonlinear relation between the numerical and functional responses, which on the population level can be described by a non-constant conversion factor $\gamma(X,Y)$.

Likewise, if searching and handling predator individuals have different death rates and their relative densities vary with the prey density, then the mean death rate will vary accordingly. Moreover, if two prey states (like exposed and hiding) have different birth or death rates, and if the relative densities of the states depend on the predator density (as in a previous example), then the numerical response of the prey will depend on the predator density in a way that is not directly linked to prey capture. All in all this leads to the more general model

[TABLE]

Once the states and state transitions have been specified, the expressions for $f$, $g$, $\gamma$ and $\delta$ as functions of both $X$ and $Y$ follow automatically.

Functional responses provide a connection between two levels of description of a biological population: the microscopic level, where the interactions between individual behavioural states are described, and a macroscopic level, where only the population size is tracked. The interplay between behavioural states and functional responses has been the focus of numerous research works, see e.g. [27, 1, 21, 2]. More recently, several teams are trying to understand the impact of stochasticity (due, for instance, to stochastic interactions between predators and prey, or to the limited number of individuals involved). For more details on stochastic models, we refer to [22], while the connection between stochastic and deterministic models is discussed in [9]. Moreover, in [6] the first deviation from the deterministic dynamics implied by stochastic effects has been analysed.

In the context of global warming and rapid changes of species range, the quantitative information provided by functional responses is a valuable tool to understand population dynamics. Recent studies on invasive species have adopted this approach, for example [11, 4, 33, 8], while other teams are using this notion to discuss the efficiency of biocontrol agents, as given in [7, 30]. In both cases, the precise description provided by the functional response proves key to understanding the effect of the antagonistic behaviour on biodiversity and species density. Other works investigate further the quantitative capabilities of functional responses: they develop fitting methods and algorithms to estimate the parameters of the models, see [28, 31, 29].

In this article, we consider that both the predator and the prey population are structured by behavioural states. We give a formal method for the mechanistic derivation of a predator’s functional response which provides a clear interpretation of the population dynamics in terms of the underlying individual level processes. In addition, we apply the same method to derive the corresponding numerical response of the predator and the prey as well.

Section 2.1 is focused on the time scale separation argument and the general method used to derive the functional and numerical responses. In Section 2.2, we discuss the existence and uniqueness issues for the fast dynamics that are necessary to apply the time scale separation idea. We illustrate these notions with a canonical example in Section 3.

In Sections 4 and 5, we consider two cases where either the prey population or the predator population is structured by behavioural states, while the other species is only described by its total density. An interesting outcome is that the conversion rate and the death rate of the predator population are no longer constant if we suppose that the behavioural states of the predator population impact its reproduction rate or mortality rate. In Section 6, we consider a case where both populations are structured by behavioural states. In spite of the more complex interaction structure, we show that it is still possible to understand the fast dynamics of the model and to derive explicitly the functional and numerical responses. In each section, we compare the functional responses that we obtain with the well known functions of Holling and Beddington and DeAngelis.

In Section 7, we discuss the advantages and drawbacks of our approach.

2 The general method

2.1 The model

Consider a predator-prey model with $x=(x_{i})_{i=1}^{m}$ and $y=(y_{i})_{i=1}^{n}$ where $x_{i}$ and $y_{i}$ denote the densities of the prey population and the predator population in the various states. By assuming that the transitions between the different states are fast, we can ignore slower processes such as birth and decay.

The ordinary differential equations which model the fast time scale scenario are given by

[TABLE]

with the consistency conditions on the parameter values

[TABLE]

The prey move from state $i$ to state $k$ with rate $A_{ki}$, with $k,i\in[1,m]$. They leave state $i$ at rate $A_{ii}$ and enter one of the $k\neq i$ states with $k\in[1,m]$ at rate $A_{ki}$ with $k\neq i,k\in[1,m]$, such that (4) follows. When interacting with predators in state $j$, prey individuals in state $i$ leave their state at rate $B_{ij}^{(i)}$ and move into one of the $k\neq i$ states, $k\in[1,m]$ at rate $B_{ij}^{(k)}$ with $k\neq i,k\in[1,m]$. Then, the consistency condition in (5) is required.

The predators move from state $i$ to state $k$, with $i,k\in[1,n]$, by interacting with the prey in state $j\in[1,m]$ with rate $C_{ij}^{(k)}$ or spontaneously with rate $D_{ki}$. As for the prey transitions, the consistency conditions on the parameter values which characterise the interactions between the predator states are given in (6) and (7).

Note however that the model is not fully general since it does not include prey-prey or predator-predator interactions, nor the formation of complexes of (possibly multiple) prey or predator individuals, such as the formation of groups (see e.g. [16]). For example, stalking of the prey and unsuccessful attacks can be modelled as fast reactions, but they involve mixed predator-prey states. These cases fall outside the general framework presented here, as we include only pure prey and pure predator states. However, while including prey-prey or predator-predator interactions as well as mixed predator-prey states is straightforward for concrete applications (see, for example, [21]), it becomes more difficult to give a general qualitative analysis of the fast dynamics. Moreover, we consider only two separate time scales. The approach can be readily extended to multiple time scales.

In order to derive the corresponding functional response on the slow time scale, it is necessary that the fast dynamics settles on a unique hyperbolically stable steady state $(\bf\hat{x},\bf\hat{y})$, where $\bf\hat{x}$ and $\bf\hat{y}$ denote the population column vectors.

We calculate the functional response by considering the total number of prey in state $i$ caught with capture rate $\beta_{ij}$ by an individual predator in state $j$ over the size $Y$ of the predator population

[TABLE]

with $X$ denoting the size of the prey population.

Likewise, the prey numerical response can be expressed as

[TABLE]

where $\lambda_{i}$ and $\mu_{i}$ are respectively the per capita birth and natural mortality rates corresponding to the prey state $i$. Next, for the predator, let $\gamma_{ij}$ and $\delta_{j}$ denote respectively the per capita fecundity of a predator in state $j$ that has captured a prey in state $i$ and the per capita natural mortality rate of a predator in state $j$. Then the predator’s numerical response is given by

[TABLE]

where $\gamma(X,Y)=\frac{\sum_{i=1}^{m}\sum_{j=1}^{n}\gamma_{ij}\beta_{ij}\hat{x}_{i}\hat{y}_{j}}{\sum_{i=1}^{m}\sum_{j=1}^{n}\beta_{ij}\hat{x}_{i}\hat{y}_{j}}$ and $\delta(X,Y)=\frac{\sum_{j=1}^{n}\delta_{j}\hat{y}_{j}}{Y}$ are respectively the density-dependent conversion factor and mortality rate.

The equilibrium on the short time scale gives the frequency distribution of individuals over the various states. This is the same as the distribution of the amount of time that a single individual spends in the various states. Therefore, the population level responses on the long time scale result from time-averaging on the short time scale.

Note that the function $\gamma$ does not have a direct interpretation in terms of the individual behaviour. It is a population level model component that we introduce here in order to give the equation in the form of the predator-prey model in (1). The functional form of the product of $\gamma$ and $f$, on the other hand, does have an individual level interpretation, which is given in equation (10).

2.2 Existence and uniqueness of the fast dynamics equilibrium

The system in (3) can be rewritten in matrix form as

[TABLE]

The matrices $\bf A+B(y)$, $\bf A$ and $\bf B(y)$ in $M_{m}(\mathbb{R})$, where we use the notation $M_{m}(\mathbb{R})$ to denote the $m\times m$-matrix space over $\mathbb{R}$, are non-negative off-diagonal matrices and have negative main diagonal entries. The same conditions apply to the matrices $\bf C(x)+D$, $\bf D$ and $\bf C(x)$ in $M_{n}(\mathbb{R})$.

In the linear case, when $\bf{B}=0,\bf{C}=0$, because of the consistency conditions, the matrices $\bf A$ and $\bf D$ correspond to the transition rate matrices of a continuous time Markov chain and the system in (11) becomes

[TABLE]

Under the assumption that this Markov chain is irreducible and aperiodic, there exists a unique stationary distribution ${\bf\pi}$, corresponding to the fast dynamics steady state we are looking for and which can be found by solving the system in (11). Furthermore, the convergence to the limit distribution is exponentially fast. As shown in the example in Section 3, a similar argument is used in the triangular case, when the transitions of one of the two species are not affected by the other population densities.

Alternatively, when we consider the non-linear case ($B,C\neq 0$), the existence of the equilibrium corresponding to the fast dynamics is guaranteed by the Perron-Frobenius Theorem and Shauder’s Fixed Point Theorem. We give the detailed proof in the Appendix A. However, the uniqueness of this equilibrium is more difficult to establish.

When we consider a model with a small number of states (typically 4 states in total), the uniqueness and hyperbolic stability of the steady states can often be checked directly. This is the case of the application that we will discuss in Section 6, where the hyperbolic stability of the fast dynamics equilibrium is verified, as we show in the Appendix C.

If the number of states is larger (more than 4 states in total), we have not been able to prove or refute the uniqueness and hyperbolic stability of the steady states under the assumptions presented in Section 2.1. If we relax these assumptions, however, we can show that it is possible to build examples where the uniqueness does not hold. In Appendix D, we construct matrices $\bf A,B,C,D$ that satisfy the assumptions except for some diagonal coefficients of $\bf A$ and $\bf D$ that are equal to [math] (specifically $A_{11}=A_{nn}=D_{11}=D_{nn}=0$). The fast dynamics has then at least two different steady states. This example, where the uniqueness is an issue, can be seen as a model for an actual biological system. Therefore the uniqueness problem appears not only as a mathematical challenge, but also as an important question for the general application of the method, which first of all requires a good understanding of the fast dynamics asymptotic behaviour.

The example discussed in Section 3 is a triangular case of (11), where $\bf B(y)=0$ for all $\bf y$, while the applications in Sections 4 and 5 model the scenario with $\bf A=0$ (i.e., only a single prey state) and $\bf B(y)=0$ for all $\bf y$ and Section 6 gives an application with the complete model form.

3 Application: when the transitions of the prey are not due to the interactions with the predator states.

Consider the following ODE system

[TABLE]

The equations in (13) do not take into account the movements between different prey states due to the interactions with the predator states, while the predator movements can be induced by the prey states.

The first equation can be rewritten as $\dot{\bf x}=\bf A\cdot x$, where $\bf x$ and $\dot{\bf x}$ are the column vectors of respectively the different prey states $\left(x_{i}\right)_{i=1}^{m}$ and their derivatives $\left(\frac{dx_{i}}{dt}\right)_{i=1}^{m}$ and $\bf A$ is the matrix in $M_{m}(\mathbb{R})$ with elements $A_{ij}$. The elements on the main diagonal are all strictly negative, while the elements off the main diagonal are non-negative. Furthermore, each column sums to zero. The negative diagonal entry represents the lifetime rate of the corresponding state, while the off-diagonal entries are proportional to the transition probabilities of the embedded Markov chain of the continuous time Markov jump process.

The equilibrium $\hat{\bf x}$ for the fast dynamics satisfies the equation $\bf A\cdot\hat{x}=0$. We denote by $\bf E_{m}$ the matrix of ones in $M_{m}(\mathbb{R})$. Let $X$ and $Y$ be the total population sizes. Then, $\hat{\bf x}$ satisfies the equation $\bf E_{m}\cdot\hat{\bf x}=\bar{X}$, where $\bar{\bf X}$ denotes the vector of length $m$ with all elements equal to $X$. Summing up the two equations, we obtain that $\bf\hat{x}=(A+E_{m})^{-1}\cdot\bar{X}$.

We repeat the same procedure with the system of equations given by $\bf(C(\hat{x})+D)\cdot\hat{y}=0$, where $\hat{\bf y}$ denotes the column vector of length $n$ of the predator states at equilibrium and $\bf(C(\hat{x})+D)$ the matrix in $M_{n}(\mathbb{R})$ evaluated at $\hat{\bf x}$. Furthermore, $\hat{\bf y}$ satisfies $\bf E_{n}\cdot\hat{y}=\bar{Y}$, where $\bar{\bf Y}$ denotes the vector of length $n$ with all elements equal to the total predator density $Y$. Then, the equilibrium is given by $\bf\hat{y}=(C(\hat{x})+D+E_{n})^{-1}\cdot\bar{Y}$.

Now we apply the above to the following matrix for the prey states

[TABLE]

The matrix is analogous to the transpose of the generator matrix corresponding to a birth-death process where the parameters $A_{k+1,k}$ and $A_{k,k+1}$ are respectively the birth rates and the death rates. In this case, the prey transitions are not predator induced and occur only between neighbouring states. The prey leave the class $x_{k}$ with rate $A_{kk}=A_{k-1,k}+A_{k+1,k}$. Moreover, the transitions from the $(k-1)$-state and $(k+1)$-state to the $k$-state happen at rates $A_{k,k-1}$ and $A_{k,k+1}$, as shown below:

[TABLE]

The ordinary differential equations for the $k=2,...,m-1$ states are

[TABLE]

while the fast dynamics of the $k=1,m$ states is modelled by

[TABLE]

We use the solution of the balance equation for the stationary distribution of the birth-death process. The equilibrium of the system $\bf\dot{x}=A\cdot x$ is of the form

[TABLE]

Given the normalisation condition $\sum_{k=1}^{m}x_{k}=X$ for the total prey population, we obtain

[TABLE]

We now set $A_{i+1,i}=A_{i}Y$, that is, we assume that the prey transitions from each state to the consecutive one are directly proportional to the total predator density. If, for example, consecutive prey states represent increasing levels of protection, then the prey are likely to move from the current state to the next one at higher predator densities. This applies to the context of prey defenses triggered by predator kairomones, e.g. chemo-signals which warn the prey of danger (see, for example, [26, 3, 32]). We denote $\prod_{i=1}^{k-1}\frac{A_{i}Y}{A_{i,i+1}}$ by $\mathcal{A}_{k}Y^{k-1}$. Then:

[TABLE]

[TABLE]

Additionally, we assume that the predator has two states: searching and handling. In particular, we assume that the handling state includes every action of the predator that occurs after prey capture, such as the actual killing of the prey, carrying the prey to the lair, eating, digesting, resting and giving birth. On the other hand, the searching state is considered as an highly active state. Therefore, births happen only while the mother is in the handling state, although rarely enough to be negligible on the short time scale in order not to violate the assumption of constant total prey size. The ODE for the searching predators with density $S$ and with attack rates $(c_{i})_{i=1}^{m}$ corresponding to each prey state is given by

[TABLE]

where $H$ is the density of the handling predators, $\frac{1}{d}$ is the average handling time and $S+H=Y$.

The fast dynamics equilibrium for the searching predators and the handling predators is

[TABLE]

[TABLE]

Since we include prey capture as a fast process, we need that the predator population size is much smaller than the prey population size, i.e. $Y\ll X$, so that the effect of the prey capture on the total prey population size $X$ is negligible. The prey capture is proportional to the predator population size. If we did not make the assumption of rare predators, then the prey would die out on the short timescale. The total population size $Y$ in this example and in the applications in Section 4, 5 and 6 is always the magnified or scaled-up predator population size. In Appendix B we give the technical details and assumptions about the individual behaviour in order to achieve time scale separation between the fast dynamics and the slow dynamics.

The predator functional response $f(X,Y)$ is calculated as in (8) and given by

[TABLE]

If $\frac{1}{d}X\left(c_{1}+\sum_{k=2}^{m}c_{k}\mathcal{A}_{k}Y^{k-1}\right)\gg 1+\sum_{k=2}^{m}\mathcal{A}_{k}Y^{k-1}$, the food source is superabundant and an increase in the prey density $X$ does not increase the feeding rate, which reaches a constant saturation level $d$, as in the Holling type II functional response. Then the function in (25) is increasing with the total prey population until it saturates at this value.

Furthermore if there is only one prey state the functional response in (25) simplifies to the Holling type II functional response

[TABLE]

We recall that in (25) we assume $A_{21}=A_{1}Y$, that is the rate which determines the transitions of the prey from the defended state $x_{1}$ to the exposed state $x_{2}$ is linearly depending on the total predator density. This interpretation agrees also with the individual behaviours modelled by Geritz and Gyllenberg in [15] for their mechanistic derivation of the Beddington-DeAngelis functional response, where the available prey are in state $x_{1}$, while the $x_{2}$ class denotes those individuals that found a refuge from the predators. In particular, in the literature the most common form of the function by Beddington and DeAngelis is given by

[TABLE]

Here we obtain a generalisation of the Beddington-DeAngelis functional response in (27), which differs from the one in [15] because we suppose that the prey in both states $x_{1}$ and $x_{2}$ can be captured but at different rates:

[TABLE]

[TABLE]

The graph of the function in (29) is illustrated in Fig. 1. When the prey population $X$ increases, the asymptotic behaviour is the same as described above for the class of functions in (25). At high predator density $Y$, the functional response in (29) with $c_{1}<c_{2}$ increases and tends to the Holling type II functional response (Fig. 1(a)). This is in contrast to the Beddington-DeAngelis functional response, which is a decreasing function of $Y$. On the other hand, if $c_{1}>c_{2}$ the functional response decreases with the total predator size $Y$ (Fig. 1(b)).

If $c_{2}=0$ (or sufficiently small), the functional response in (29) is a Beddington-DeAngelis functional response of the form given in (27) with $a=c_{1}A_{12}$, $b=c_{1}\frac{1}{d}$ and $c=\frac{A_{1}}{A_{12}}$.

If $c_{1}=0$ (or sufficiently small), then the slow predator dynamics will have an Allee effect, i.e., the predator population cannot grow when its density is below a given threshold. In particular, it would not be able to invade the predator-free population. This is because of the $XY$ term in the numerator of the functional response, which gives a squared $Y$ term for the population level birth rate of the predator. Then, at low predator densities, the predator birth term is dominated by the negative linear term describing predator death. Allee effects in predator-prey systems lead to homoclinic or heteroclinic bifurcations. At the individual level, we assumed that $A_{21}=A_{1}Y$, that is the rate at which the prey enter state $x_{2}$ is proportional to the total predator density. If $c_{1}=0$ and $c_{2}>0$, then the predators consume only the prey in state $x_{2}$. However, if the total predator density $Y$ is small, there are not enough prey in state $x_{2}$ for the predators to survive.

4 Application: a functional response with density dependent handling time

4.1 Individual level reactions, population equations and fast equilibrium

We analyse the same scenario presented in Section 3, that is, the predator population is structured in two classes, the searching predators $S$ and handling predators $H$. In this application we consider only one prey state. We define with $c_{1}$ the attack rate. We assume moreover that the handlers return to the searching state with prey density dependent rate $c_{2}X$ or spontaneously with rate $d$. This assumption is ecologically reasonable if the uptake of resources from the corpse of the killed prey declines with the handling time. The capture of a new prey becomes then worthwhile especially if the prey density is high and an indicator of the overall prey density is given by the average time until a new prey gets into the handling predator’s field of vision. Furthermore, the density dependent transition may be the result of an actual encounter with a prey individual or may be triggered by prey kairomones (as assumed in the previous section, but with the roles of prey and predator reversed). All these interactions are fast time processes in comparison to birth and natural death and are summarised below:

[TABLE]

Note that in the first transition the prey disappears due to prey capture, while in the second the prey acts merely as a catalyst.

The corresponding population level differential equations of the fast time dynamics are given by applying the law of mass action and the time scale separation is presented in details in Appendix B:

[TABLE]

The total predator population is constant and given by $Y=S+H$. Then, we can reduce the system of equations to only one equation and solve the steady state equation. The fast dynamics equilibrium is $(\hat{S},\hat{H})=\left(\frac{(d+c_{2}X)Y}{d+(c_{1}+c_{2})X},\frac{c_{1}XY}{d+(c_{1}+c_{2})X}\right)$.

4.2 Functional response

We can now derive the corresponding functional response

[TABLE]

The functional response in (32) is a two-parameter function, because only the ratios $\frac{c_{1}}{d}$ and $\frac{c_{2}}{d}$ matter and the shape of the functional response is therefore affected by these fractions, as shown in Fig. 2. Furthermore, by considering the formulation

[TABLE]

we observe that the functional response in (32) is like the Holling type II functional response, but with a density dependent handling time given by the ratio $\frac{1}{c_{2}X+d}$. The higher the prey density is, the faster the predator will quit handling and start searching for fresh food. Note that such behaviour is functional because if the prey is scarce, then the predators will tend to diligently consume the food source until it is completely exhausted and handle the prey longer than if the prey were abundant.

4.3 Predator numerical response

We assume that the predators in the handling state give birth. Therefore the per capita reproduction rate of the predators is proportional to

[TABLE]

The conversion factor $\gamma(X)=\frac{1}{c_{2}X+d}$ is then a decreasing function of the prey density $X$ and in particular it is proportional to the average time spent handling the prey. If we assume that in this particular scenario the predator per capita natural mortality rate is the same for the predators both in the searching state and in the handling state, then the mortality rate is constant as in (1).

5 Application: type III functional response and corresponding predator numerical response

5.1 Individual level reactions, population equations and fast equilibrium

We now assume that the searching predators are divided into two subclasses according to the level of starvation, well-fed ($S_{1}$) and starving ($S_{2}$). It is natural then to assume different capture rates for the two classes, e.g. starving predators have a lower capture rate than satiated predators, $c_{1}>c_{2}$. We suppose again the handling predators in class $H$. We denote with $H_{1}$ the predators that enter state $H$ from state $S_{1}$ and with $H_{2}$ the predators that enter $H$ from state $S_{2}$. The predators in $H_{1}$ and $H_{2}$ handle the prey for, on average, $\frac{1}{d_{1}}$ units of time, then they enter class $S_{1}$. If a well-fed predator does not capture a prey in on average $\frac{1}{d_{2}}$ units of time, it enters the state $S_{2}$. We assume that starving predators have very low per capita fecundity since they are hunting to survive and restock their basic energy reserve. On the contrary the well-fed predators invest part of the energy gained from the food source for reproduction. In this case, the per capita fecundities for the two types of consumers, namely $\Gamma_{1}$ and $\Gamma_{2}$, differ. In particular $\Gamma_{1}>\Gamma_{2}$, since in case of starvation the individuals are likely to cease energy allocation to reproduction, see [23]. We assume the offspring to be in state $S_{2}$. We consider the transitions between the different predator states to be fast processes with respect to birth and death. In addition, we consider the total predator population size considerably smaller than the total prey size. The time scale separation is achieved through the scaling given in Appendix B. We summarise below the individual level processes:

[TABLE]

The equations that describe the population level fast time dynamics are given by

[TABLE]

If the total predator density is constant, then the conservation law $\frac{dY}{dt}=\frac{dS_{1}}{dt}+\frac{dS_{2}}{dt}+\frac{dH_{1}}{dt}+\frac{dH_{2}}{dt}=0$ holds. The fast dynamics is settled on the asymptotically stable equilibrium

[TABLE]

5.2 Functional response

The functional response corresponding to the dynamics above is

[TABLE]

The functional response in (38) is an increasing function of $X$ up to a saturating level given by $d_{1}$, as shown in Fig. 3(a). It is a type III functional response of the form

[TABLE]

with $a=c_{2}$, $b=c_{1}c_{2}\frac{1}{d_{2}}$, $c=c_{2}\left(\frac{1}{d_{1}}+\frac{1}{d_{2}}\right)$ and $d=\frac{c_{1}c_{2}}{d_{1}d_{2}}$. A necessary condition for the function to be convex in the neighbourhood of [math] is that the second derivative has to be positive. This is true if and only if $\frac{c_{1}}{c_{2}}>1+d_{2}\frac{1}{d_{1}}$, that is the attack rate of the well-fed predators, $c_{1}$, is sufficiently higher than the attack rate of the starving predators, $c_{2}$. This result is consistent with the biological interpretation of the individual level dynamics that we have provided.

In the literature the most common form of the Holling type III functional response is

[TABLE]

The function in (39) can be mathematically reduced to the function in (40), if we let $c_{2}\rightarrow 0$ and $c_{1}\rightarrow\infty$, such that the product $c_{1}c_{2}$ stays constant. This can be interpreted as well-fed predators being extremely efficient hunters, while starving predators are very unsuccessful searchers.

5.3 Predator numerical response

On the slow time scale, we suppose that the reproduction rate of the predators is proportional to the density of the handling predators at the fast time equilibrium and that the energy intake from the consumption of the prey is partly allocated to reproduction. At an individual level, we consider the following reactions, which happen at a slow time scale with respect to the interactions modelled in (35):

[TABLE]

The per capita reproduction rate in this particular scenario is not a constant, but it is an increasing function of the total prey population size (see Fig. 3(c)) with saturating value given by the fecundity rate of the well-fed predator individuals, $\Gamma_{1}$, at high prey density:

[TABLE]

The conversion factor

[TABLE]

is a function of the prey density (see Fig. 3(b)), saturating on the value $\frac{\Gamma_{1}}{d_{1}}$ when the food source is abundant. In fact, $\frac{1}{d_{1}}$ is the average time spent handling the prey and $\Gamma_{1}$ denotes the per capita fecundity of the well-fed predators. Furthermore, the function $\gamma(X)$ is increasing if and only if $\Gamma_{1}>\Gamma_{2}$, which is consistent with the biological assumptions given in Section 5.1.

Suppose further that the predators in the two searching states differ not only in their capture rates, but also in their respective mortality rates $\delta_{1}$ and $\delta_{2}$, with $\delta_{2}>\delta_{1}$. The individual level interactions occurring at the slow time scale which model the natural mortality of the predators are then given by

[TABLE]

Under these assumptions we note that the average per capita mortality rate is given by

[TABLE]

Mortality is then no longer constant, but a decreasing function of the prey density, as shown in Fig. 3(d). In particular, in the absence of the prey, the function $\delta$ takes the value of the per capita mortality rate of the starving predators, $\delta_{2}$.

5.4 Alternative interpretation of the individual level processes

As far as the functional response is concerned, the model can be interpreted also in the context of predators structured according to their level of experience. This is the usual interpretation of the Holling type III functional response, although until now a rigorous derivation was missing. Here we show how the Holling type III functional response can be derived with our approach. Suppose that class $S_{2}$ contains the individuals that lack experience and are not well skilled in capturing the prey, while class $S_{1}$ includes those experienced individuals with success rate $c_{1}>c_{2}$. Note that at low prey density, almost all predators will be inexperienced. Predators in both classes, after interaction with the prey, enter the class of handling predators $H$. The average time spent handling the prey is $\frac{1}{d_{1}}$ units of time. Predators that have captured (and handled) a prey are considered experienced (class $S_{1}$), but they lose this status, and hence the ability to capture prey at high rate, after on average $\frac{1}{d_{2}}$ units of time (transition to class $S_{2}$). Furthermore, we assume that the individual level processes that determine the interactions between the predator and prey states are fast time reactions with respect to birth and death. The above scenario can be visualised in the following way:

[TABLE]

The ODE system which describes the population level fast time dynamics is the following

[TABLE]

The total predator density $Y$ is constant, such that $\frac{dY}{dt}=\frac{dS_{1}}{dt}+\frac{dS_{2}}{dt}+\frac{dH}{dt}=0$ Then the fast dynamics settles on the asymptotically stable equilibrium:

[TABLE]

The corresponding functional response has already been given in (38).

6 Application: a functional response that induces an Allee effect in the predator population dynamics

6.1 Individual level reactions, population equations and fast equilibrium

In the following model, we assume that the prey has a natural tendency to seek protection, but that searching predators are able to overcome the prey defenses, for example by causing them to panic and by attacking the isolated individuals. As shown in [15], this leads to fast-processes that are the opposite of those in the Beddington-DeAngelis model. The individual level reactions are:

[TABLE]

The corresponding differential equations for the fast time population dynamics are

[TABLE]

with the conservation laws $E+P=X$ and $S+H=Y$ (see Appendix B for details on the time scale separation between the fast and slow dynamics). The corresponding fast time equilibrium is given by

[TABLE]

where $p=\frac{a}{b}$, $q=\frac{d}{c}$ and $\Delta_{p,q}(X,Y)=q\left(p^{2}q+2p\left(q+2X\right)Y+qY^{2}\right)$. In the Appendix C, we give the phase portrait corresponding to the system in (50), for different values of the parameters, in order to show that the fast dynamics equilibrium is unique and hyperbolically stable.

6.2 Functional response

We derive the corresponding functional response (Fig. 4(a))

[TABLE]

When $p\rightarrow 0$, that is $a\rightarrow 0$ or $b\rightarrow\infty$, the functional response tends to the Holling type II functional response, since the prey is most of the time available for being captured:

[TABLE]

The case in which $q\rightarrow\infty$, that is $c\rightarrow 0$ or $d\rightarrow\infty$, corresponds to the scenario where the predators are almost all the time searching. Therefore, the predators searching and attacking the prey correspond to the total predators $Y$ and the portion of prey subjected to predation is given by $\frac{bY}{a+bY}X=\frac{Y}{p+Y}X$, with $\frac{bY}{a+bY}$ being the probability for the prey to be in the vulnerable state. Taylor expanding with respect to $\frac{1}{q}$ near zero and retaining only the lowest order term in $\frac{1}{q}$, we get (Fig. 4(b))

[TABLE]

We note that if both $q\rightarrow\infty$ and $p\rightarrow 0$, then the functional response in (55) becomes linear, as in the Holling type I functional response. In this case the predators handle the prey for an infinitely short time and the prey is most of the time exposed to the predators’ attacks. The function is increasing with the attack rate $c$. This is a typical functional response for filter feeders as shown in [jeschke2004consumer].

6.3 Prey and predator numerical responses

If we assume that only the handling predators give birth, then the predator per capita birth rate is proportional to

[TABLE]

i.e. it is proportional to the functional response, as usual. On the other hand, when we consider the searching predators $S$ and the handling predators $H$ having different death rates $\delta_{1}\neq\delta_{2}$, then the overall per capita death rate is

[TABLE]

which is no longer constant, as usually given in the literature, but depends on $X$ and $Y$.

We note that the product of $XY$ in the numerator of the functional response in (57) leads to a squared $Y$ term in the population equations for the prey and the predator. In the predator equation for the slow dynamics, this leads to an Allee effect (i. e. at low predator densities almost all prey are protected and cannot be captured), since the per capita birth rate at low predator densities is of order $Y$ (see (58)), while, when $q\rightarrow\infty$ and the density of the handling predators is very small, the per capita death rate is approximated by the constant value $\delta_{1}$ (see (59)). A similar situation was modelled from first principles in [16].

If we assume that only the protected prey in class $P$ give birth, i.e., if no other slow interactions between the individual prey or their resources occur, then the prey per capita birth rate is proportional to

[TABLE]

However, when the hiding prey and the available prey have different death rates $\mu_{1}\neq\mu_{2}$ and if there are no other sources of slow death, e.g. interference competition among the prey, then the overall per capita death rate is

[TABLE]

Furthermore, we note that by visualising the prey and predator numerical responses in terms of the functional response, it is possible to understand what the former look like in the limiting cases of the latter as treated above in (56) and (57).

7 Conclusions

In this paper, we proposed a method for the derivation of the functional response from a system of prey-predator interactions which occur on a fast time scale, with respect to birth and death. Many functional responses appear in the literature, but they often lack of interpretation at the individual level. The time scale separation argument that we use in this paper is a possible approach to link the macroscopic behaviour of the population to the microscopic dynamics of the state transitions of individuals. Such derivation permits an explicit interpretation of the structure and parameters of the functional response in terms of the individual behaviour.

Elements of the two time-scales are implicit in the traditional approach to deriving the Holling type II functional response, as the consumption rate of the predator instantaneously adjusts to the current prey density. Here we formalise the two-time scales approach in a systematic way. Specific instances of this method can be found in the literature, for example in the works by Metz and Diekmann [25] and Geritz and Gyllenberg [15, 16]. However, in this paper we embed these instances into a more general and formal framework that, in addition to the predator’s functional response, also gives a derivation for the numerical responses of the predator as well as of the prey.

In addition to a general outline of the method, we give several concrete applications, including an application that leads to a generalisation of the Holling type III functional response. The functional response Holling type III has been associated with switching between alternative prey depending on their relative abundance. An explicit derivation was given by Van Leeuwen et al. [24]. Alternatively, the Holling type III functional response can be associated with different hunger states of the searching predators instead of different experience levels. We give here a mechanistic explanation for the latter. The specific form, as found in the literature, is recovered as a limiting case and is easily understood in terms of the explicitly modelled underlying individual behaviour.

In another application, the handling predator may abandon its catch if it detects another live prey. This leads to a Holling type II functional response with density-dependent handling time. In particular, both the handling time and the conversion factor are decreasing functions of the prey density. Such behaviour is adaptive if the uptake of resources from the killed prey declines with the handling time.

Furthermore, we discuss the functional response corresponding to a simple non-linear system for the fast dynamics, where we consider two states for the predators and for the prey and we are able to compute explicitly the fast dynamics equilibrium. Here the predators may overcome the prey defenses by causing panic among the prey and by attacking the isolated individuals. We model the prey and predator numerical responses by assuming the two species structured by states with different birth and death rates. The results at the population level are consistent with the individual level reactions and show that at low predator densities an Allee effect is likely to appear.

The method presented here is not the most general method possible. For example, we did not consider interactions among the prey themselves or the predators themselves like the exchange of information about the presence of prey or predators leading to a change in the motivational state or the state of alertness. Neither did we include states involving more than one individual, such as two predators fighting over a kill, or several prey seeking protection in numbers, or a predator stalking or fighting a prey. It is not difficult to extend the method to include these cases (e.g. see [16] where prey groups of different sizes are modelled as different prey states), but it becomes more difficult to prove the existence and, in particular, the uniqueness of an equilibrium of the fast dynamics of the state transitions.

In order to apply slow-fast time scale separation, it is necessary that the fast dynamics is settled on a unique and hyperbolically stable steady state. We are not able to give a general result. In particular, the uniqueness of the equilibrium corresponding to the fast dynamics remains an open question. However, by relaxing the conditions on the parameter values, we have built an example in which the fast time steady state is not unique: this may set a limit to the assumptions that we can make on the coefficients of the matrices which model the interactions in order to get a unique hyperbolically stable equilibrium.

In our approach only individuals in some specific discrete states are able to reproduce. The proportion of time that an individual spends in these states is equal to the proportion of individuals in such states at the fast time equilibrium. In this way birth is limited by a time-budget. Another (and possibly more realistic) approach would be to model births as energy limited, like in the dynamic energy budget models (see, for example, [23, 17]).

A possible disadvantage of the approach is that to use it in a practical way, one must make assumptions about the transitions between microstates and the model may become parameter heavy. However, as a theoretical tool, the method has potential. One of the key questions in ecology today, raised by Durrett and Levin [12] among others, is how to scale up from the level of individual behaviours in a population to functional responses and dynamics equations at the population level. Deriving functional and numerical responses from the behaviour of the individual prey and predator is important if one wants to go beyond a mere description of the population dynamics to an understanding in terms of the underlying individual level processes. Also the other way around, that is, if one wants to know the effect of certain changes in the behaviour of the individual prey or predator, the derivation of the population model from first principles in terms of individual behaviour is a necessity.

Appendix A Appendix

Proposition A.1**.**

*Let $A,(B_{i,j}^{k})_{i,j\in\{1,\dots,m\}}\in M_{m}(\mathbb{R})$ and $(C_{i,j}^{k})_{i,j\in\{1,\dots,n\}},D\in M_{n}(\mathbb{R})$ be matrices with non-negative off-diagonal coefficients. Suppose that $A$ and $D$ are transition matrices such that the linear system in (12) has a unique stable equilibrium and $(B_{i,j}^{k})_{i,j\in\{1,\dots,m\}},(C_{i,j}^{k})_{i,j\in\{1,\dots,n\}}$ are irreducible matrices respectively for all $y>0,y\in\mathbb{R}^{n}$ and for all $x>0,x\in\mathbb{R}^{m}$. Assume, moreover, that all these matrices are transition matrices and that the conservation laws on the total population density $\sum_{i=1}^{m}x_{i}=X$ and $\sum_{i=1}^{n}y_{i}=Y$ hold, with $X$ and $Y$ constant.

Then, the system in (11) has at least one equilibrium point.*

Proof.

Consider the steady state equations

[TABLE]

In the first set of equations, $\bf A+B(y)$ $\in M_{m}(\mathbb{R})$ is an irreducible nonnegative off-diagonal matrix for all $\bf y>$ [math], $\bf y$ $\in\mathbb{R}^{n}$ and by the Perron-Frobenius Theorem it has a simple dominant nonnegative eigenvalue, that is [math]. Let $\psi(\bf y)$ be the corresponding eigenvector satisfying $\sum_{i=1}^{m}\psi(y)_{i}=X$.

In the same way, for the second set of equations we have that $\bf C(x)+D$ $\in M_{n}(\mathbb{R})$ is an irreducible nonnegative off-diagonal matrix for all $\bf x$ $>0$, $\bf x$ $\in\mathbb{R}^{m}$ and by the Perron-Frobenius Theorem it has a simple dominant nonnegative eigenvalue, that is [math]. Let $\bf\phi(x)$ be the corresponding eigenvector satisfying $\sum_{i=1}^{n}\phi(x)_{i}=Y$.

The continuous map $\Phi$ from the compact convex set $\{\vec{x}\in\mathbb{R}^{m}:\sum_{i=1}^{m}x_{i}=X\}$ in itself, $\Phi=\psi\circ\phi$ has at least one fixed point by the Shauder’s Fixed Point Theorem. This shows that the fast dynamics has at least one steady state. ∎

Appendix B Appendix

We give here the time scale separation for the system in Section 3 in details. The dynamical system of the interactions modelled in Section 3 is given by

[TABLE]

where $\lambda_{k}$ and $\mu_{k}$ are respectively the per capita birth and natural mortality rate for the prey in state $k$ as given in (9), $\Gamma$ is the conversion rate of prey into predators such that $\Gamma H=\gamma(X,Y)f(X,Y)Y$ and $\delta$ is the per capita mortality rate of the predators.

Let $\varepsilon>0$ be a small and dimensionless scaling parameter. In order to separate the fast and slow dynamics, we define the following scalings for the parameters of fast time interactions and the predator population that is assumed to be much smaller than the prey population: $\tilde{A}_{k,k-1}=\varepsilon^{-1}{A}_{k,k-1}$, $\tilde{A}_{kk}=\varepsilon^{-1}{A}_{kk}$, $\tilde{A}_{k,k+1}=\varepsilon^{-1}{A}_{k,k+1}$, $\tilde{c}_{k}=\varepsilon^{-1}{c}_{k}$, $\tilde{d}=\varepsilon^{-1}d$, $\tilde{Y}=\varepsilon Y$, $\tilde{H}=\varepsilon H$, $\tilde{S}=\varepsilon S$. We now give the slow-fast equations corresponding to the system in (63) using the scaled parameters:

[TABLE]

We introduce the scaled short time $t=\varepsilon^{-1}\tau$ and let $\varepsilon\rightarrow 0$. We give the equations for the dynamics on the fast time scale:

[TABLE]

The variables $X$ and $Y$ are constants on the fast time scale and from the equations in (65) for the fast variables $x_{k}$, $S$ and $H$ we can now derive the fast dynamics equilibria given in (20), (21), (23) and (24).

The time scale separations for the models in Sections 4 and 5 follow the passages and the scalings given above for the system in Section 3.

We consider now the dynamical system for the interactions modelled in Section 6:

[TABLE]

In this case, as in [15], we additionally assume that $\tilde{b}$ is large in comparison to the other parameters. In this way, the term $-bPS$ in the short time scale equations for the exposed and the protected prey is not negligible, but part of the fast dynamics. We introduce the scaled parameter $\tilde{b}=\varepsilon^{-2}b$. For an alternative scaling, we could assume that the parameters $\tilde{c}$ and $\tilde{d}$ are small in comparison to the other parameters (see also [15]). In this case, we would not need to assume that the total predator size is much smaller than the total prey size, but we would separate the dynamics in (66) into three separate time scales.

Appendix C Appendix

We compute the solutions for the system of equations in Section 6.1, for different values of the parameters $a$, $b$, $c$, $d$ and different initial conditions. The numerical simulations in Fig. 5 show that the uniqueness and hyperbolic stability of the fast dynamics steady state is verified for the chosen values of the parameters.

Appendix D Appendix

D.1 A counterexample to the uniqueness of the fast dynamics steady state

Consider the system in (3). Suppose $\vec{A},(\vec{B}_{i,j}^{k})_{i,j\in\{1,\dots,n\}},(\vec{C}_{i,j}^{k})_{i,j\in\{1,\dots,n\}},\vec{D}\in M_{n}(\mathbb{R})$.

Consider weaker assumptions on $(\vec{B}_{i,j}^{k})_{i,j\in\{1,\dots,n\}}$ and $(\vec{C}_{i,j}^{k})_{i,j\in\{1,\dots,n\}}$, such that they are not non-negative off-diagonal matrices, but verify for every $i,j\in\{1,\dots,n\}$

[TABLE]

In this particular case, we are able to give numerically a counterexample to the uniqueness of the steady state of the system.

We consider the following symmetric system, where $\vec{A}\equiv\vec{D}$ and $(\vec{B}_{i,j}^{k})_{i,j\in\{1,\dots,n\}}\equiv(\vec{C}_{i,j}^{k})_{i,j\in\{1,\dots,n\}}$:

[TABLE]

This model satisfies the relations necessary for the conservation of mass: one can check that $\frac{d(x_{1}+x_{2})}{dt}(t)=\frac{d(y_{1}+y_{2})}{dt}(t)=0$ (for any $(x_{1}(t),x_{2}(t),y_{1}(t),y_{2}(t))$). Moreover, the system has two steady states, namely

[TABLE]

and

[TABLE]

One may check that, thanks to the Implicit Function Theorem, perturbations of the steady states still exist if we make all the coefficients non zero, but close to the coefficients chosen above.

Under the stronger assumptions of non-negative off-diagonal matrices $({\bf B}_{i,j}^{k})_{i\in\{1,\dots,n\}}$ and $({\bf C}_{i,j}^{k})_{i\in\{1,\dots,n\}}$, we are also able to build a counterexample to the uniqueness of the steady state corresponding to the fast dynamics in (3). In particular, we still assume a symmetric situation where the matrices describing the dynamics for the first set and the second set of equations are equivalent. Moreover, we impose that the matrices corresponding to the linear part of the system of equations are transition matrices such that state $1$ and state $n$ are absorbing states, in a stochastic sense.

We construct a counterexample based on the following cross-diffusion system:

[TABLE]

with $n,m\in C^{2}([0,1))$ and Neumann boundary conditions. To show that this system may have several solutions, one may consider $a$ satisfying

[TABLE]

Then $n=\lambda+\cos(\pi x)$, $m=\lambda+\cos(\pi(1-x))$ is a solution and the symmetry of the system makes it non unique. We have developed a discrete version of this idea and we give in the following section the corresponding numerical example. We proceed in three steps:

**Step 1: Defining $\bf x,y,T^{x},T^{y}$ such that $\bf T^{y}x=0$ and $\bf T^{x}y=0$.

**

Definition of $\bf x$ and $\bf y$.

We consider $\vec{x},\vec{y}\in\mathbb{R}^{n}$ defined by

[TABLE]

for every $i\in\{1,\dots,n\}$. Note that $\bf x$ is different from $\bf y$.

Definition of $\bf T^{y}$.

Let $\vec{T}^{\vec{y}}\in M_{n}(\mathbb{R})$ a tridiagonal matrix: for any $i\in\{1,\dots,n\}$,

[TABLE]

except for the two terms $T^{y}_{1,1}=-\frac{t^{y}_{1}}{2}$, $T^{y}_{2,1}=\frac{t^{Y}_{1}}{2}$, $T^{y}_{n,n}=-\frac{t^{y}_{n}}{2}$, $T^{y}_{n-1,n}=\frac{t^{y}_{n}}{2}$ and all the other coefficients are [math]. In these expressions, we have denoted

[TABLE]

Thanks to the definitions (74) and (73), for $i\in\{2,\dots,n-1\}$

[TABLE]

while for $i=1$ (a similar computation can be made for $i=n$),

[TABLE]

We have then shown that

[TABLE]

Definition of $\bf T^{x}$.

Let $\vec{T}^{\vec{x}}\in M_{n}(\mathbb{R})$ a tridiagonal matrix: for any $i\in\{1,\dots,n\}$,

[TABLE]

except for the two terms $T^{x}_{1,1}=-\frac{t^{x}_{1}}{2}$, $T^{x}_{2,1}=\frac{t^{x}_{1}}{2}$, $T^{x}_{n,n}=-\frac{t^{x}_{n}}{2}$, $T^{x}_{n-1,n}=\frac{t^{x}_{n}}{2}$ and all the other coefficients are [math]. In these expressions, we have denoted

[TABLE]

Then,

[TABLE]

**Step 2: Defining a linear interpolation between $\bf T^{y}$ and $\bf T^{x}$.

**

We want to define $\alpha_{i}$, $\beta_{i}$ such that

[TABLE]

where $\mu(i)=1$ if $i\leq\frac{n}{2}$ and $\mu(i)=n$ otherwise.

Case where $i\leq\frac{n}{2}$.

[TABLE]

Taking the difference of those two equations leads to

[TABLE]

which shows that $\beta_{i}$ is positive:

[TABLE]

Next, we add up the two equations appearing in (111). Thanks to (73), we have $x_{i}+y_{i}=2\lambda$ and thanks to the definition (75), (80) of $t^{y}_{i}$ and $t^{x}_{i}$,

[TABLE]

The value of $\beta_{i}$ is given by (90), then

[TABLE]

Therefore

[TABLE]

for $i\leq\frac{n}{2}$ and $i\neq 1$.

Case where $i>\frac{n}{2}$.

[TABLE]

Taking the difference of those two equations leads to

[TABLE]

which shows that $\beta_{i}$ is positive:

[TABLE]

Next, we add up the two equations appearing in (111). Thanks to (73), we have $x_{i}+y_{i}=2\lambda$ and from the definitions in (75) and (80) of $t^{y}_{i}$ and $t^{x}_{i}$ we obtain

[TABLE]

The value of $\beta_{i}$ is given by (90), then

[TABLE]

Therefore

[TABLE]

for $i>\frac{n}{2}$ and $i\neq n$.

**Step 3: Conclusion.

**

Definition of $\bf A,B,C,D$.

Let $\vec{A},\vec{B},\vec{C},\vec{D}\in M_{n}(\mathbb{R})$ tridiagonal matrices. We define the matrix $\bf A$ as follows, using the coefficients $\alpha_{i}$ defined by (92) for $i\in\{1,\dots,n\}$:

[TABLE]

Moreover we denote $A_{n,1}=\frac{\alpha_{1}}{2}$, $A_{1,n}=\frac{\alpha_{n}}{2}$, while the other coefficients are [math]. Thanks to (92), $\bf A$ is then an off-diagonal non-negative matrix, and

[TABLE]

We define next the matrix $\bf\tilde{B}$, using the coefficients $\beta_{i}$ defined by (90) for $i\in\{1,\dots,n\}$:

[TABLE]

Moreover we obtain $\tilde{B}_{n,1}=\frac{\beta_{1}}{2}$, $\tilde{B}_{1,n}=\frac{\beta_{n}}{2}$, while the other coefficients are [math]. Thanks to (90), $\bf\tilde{B}$ is an off-diagonal non-negative matrix. We define the family of matrices $(\vec{B}_{i,j}^{k})_{k}$ by

[TABLE]

For any $i\in\{1,\dots,n\}$, $(\vec{B}_{i,j}^{k})_{k,j}$ is then an off-diagonal non-negative matrix, and it satisfies for $i\in\{1,\dots,n\}$

[TABLE]

Finally, we define

[TABLE]

Note that the matrices $\vec{A},\,(\vec{B}^{k})_{k},\,(\vec{C}^{k})_{k},\,\vec{D}$ we have constructed satisfy the assumptions given at the beginning of this section.

Showing that $\bf(x,y)$ and $\bf(y,x)$ are two steady states.

For all $k\in\{1,\dots,n\}$, we use the definition of $\bf A$ and $\bf B$ to obtain

[TABLE]

Similarly,

[TABLE]

Thanks to the symmetry of the coefficients (see (98)), we also obtain

[TABLE]

[TABLE]

Finally, we have constructed two steady-states $\bf(x,y)$ and $\bf(y,x)$ for the system of differential equations that we consider:

[TABLE]

and

[TABLE]

D.2 Numerical example

We construct a numerical example where we consider the prey and the predators structured into three states. The individual level reactions correspond to the following network and reaction rates

[TABLE]

Note that $x_{1}$ and $x_{3}$, $y_{1}$ and $y_{3}$ are absorbing states for the transition matrices $\bf A$ and $\bf D$, respectively. A possible biological interpretation of the interactions in (103) is given by assuming the prey population in two different locations, in particular the prey individuals in state $x_{1}$ are in the first location, while the prey individuals in state $x_{3}$ are in the second location. The searching predators are also divided into those individuals in state $y_{1}$, which are searching for a prey in the first location, and those in state $y_{3}$, which are hunting in the second location. When a prey in state $x_{1}$ meets a predator in state $y_{1}$ (and similarly for a prey in state $x_{3}$ meeting a predator in state $y_{3}$), it goes hiding with rate $\frac{1}{6}$ and with rate $\frac{1}{4+\sqrt{2}}$ it goes back either to state $x_{1}$ or to state $x_{3}$. After an encounter with a prey individual in state $x_{1}$, with probability per unit of time $\frac{1}{6}$ the predator in state $y_{1}$ starts handling the prey. The same interactions occur for the predators in the second location when they meet a prey state $x_{3}$.

Define $x_{i}$ and $y_{i}$, with $i=1,2,3$ and $n=3$ as follows

[TABLE]

[TABLE]

Note that $\bf x$ is different from $\bf y$.

Define the matrices $\bf T^{y}$ and $\bf T^{x}$ such that $\bf T^{y}x=0$ and $\bf T^{x}y=0$ as follows

[TABLE]

[TABLE]

Then

[TABLE]

and

[TABLE]

We can define $\bf T^{x}$ in the same way and show that

[TABLE]

Then we define a linear interpolation between $\bf T^{y}$ and $\bf T^{x}$. In particular, we define $\alpha_{i}$, $\beta_{i}$ such that

[TABLE]

where $\mu(i)=1$ if $i\leq\frac{n}{2}$ and $\mu(i)=n$ otherwise.

From now on, all the numerical values are obtained by taking $\lambda=2$. This choice of the parameter values does not have any particular biological justification.

Then, by solving the above system of equations for $\alpha_{i}$ and $\beta_{i}$, we get

[TABLE]

Then, we can define the following system of ordinary differential equations

[TABLE]

The system in (113) has at least two steady-states. Thanks to the symmetry of the coefficients, $(x_{1},x_{2},x_{3},y_{1},y_{2},y_{3})=(3,2+\frac{1}{\sqrt{2}},1,1,2+\frac{1}{\sqrt{2}},3)$ and $(x_{1},x_{2},x_{3},y_{1},y_{2},y_{3})=(1,2+\frac{1}{\sqrt{2}},3,3,2+\frac{1}{\sqrt{2}},1)$ are both equilibrium points. Furthermore, the positive solutions of the steady states equations are given by

[TABLE]

or

[TABLE]

By imposing the conservation law on the sums of the $x_{i}$ and $y_{i}$, i.e. $\sum_{i=1}^{3}x_{i}=X=6+\frac{1}{\sqrt{2}}$ and $\sum_{i=1}^{3}y_{i}=Y=6+\frac{1}{\sqrt{2}}$, we get either

[TABLE]

[TABLE]

[TABLE]

or the following formulation

[TABLE]

[TABLE]

[TABLE]

The Jacobian matrix of the system (not evaluated yet at the equilibrium and not considering the conservation law on the $x_{i}$ and $y_{i}$) is given by

[TABLE]

We get always three eigenvalues with negative real part $\lambda_{1}$, $\lambda_{5}$, $\lambda_{6}$ and three null eigenvalues $\lambda_{2}$, $\lambda_{3}$, $\lambda_{4}$.

The eigenvalues zero determine a space with dimension at most three in the six dimensional one. As for the other three eigenvalues, they take negative values for every positive value of $x_{1}$ and $y_{1}$ smaller then $X$ and $Y$. The dynamics of the system converges to the manifold generated by the eigenvectors corresponding to the zero value eigenvalues.

It is not clear what is the behaviour on the stable manifold, but one could check, by numerically simulating the solutions of the system above for many initial conditions, that each point in the space generated by the eigenvectors corresponding to the null eigenvalues is an equilibrium point of the system. Therefore, we could conclude that on this space we have infinite not isolated fixed points, which are stable, but not asymptotically stable.

**Acknowledgements

**This research was funded by the Academy of Finland, Centre of Excellence in Analysis and Dynamics Research and supported by a public grant as part of the Investissement d’avenir project, reference ANR-11-LABX-0056-LMH, LabEx LMH. Part of this research was performed during the participation of three of the authors at the thematic research program "Mathematical Biology" organised by the Mittag-Leffler Institute. We also thank the two anonymous reviewers for their constructive remarks.