Interactions Obtained from Basic Mechanistic Principles: Prey Herds and Predators

Abstract

We investigate four predator–prey Rosenzweig–MacArthur models in which the prey exhibit herd behaviour and only the individuals on the edge of the herd are subjected to the predators’ attacks. The key concept is the herding index, i.e., the parameter defining the characteristic shape of the herd. We derive the population equations from the individual state transitions using the mechanistic approach and time scale separation method. We consider one predator and one prey species, linear and hyperbolic responses and the occurrence of predators’ intraspecific competition. For all models, we study the equilibria and their stability and we give the bifurcation analysis. We use standard numerical methods and the software Xppaut to obtain the one-parameter and two-parameter bifurcation diagrams.

1. Introduction

Modelling herd behaviour of a population with ordinary differential equations, via a spatial factor with herding index (see Laurie et al. [1]), avoids an explicit spatial representation.

Already in 1987, Liu and colleagues [2] used the herding index, represented by an exponent, to give a non-linear incidence rate in epidemiological models, and later this concept found several applications in the same context [3,4].

The herding index, $\alpha$ in this article, also appeared in predator–prey dynamics with prey grouped in herd. In particular, the first studies comprised only the case $\alpha =1/2$ by assuming a two-dimensional herd representation [5,6], while later works allowed for a fixed herding index in the range $[1/2,1)$ [6,7]. The concept of herding index was further extended to models other than ordinary differential equations, such as time delay equations, [8], and fractional differential equations [9,10].

In predator–prey models, the assumption that the prey gather in a herd has a direct effect on the shape of the functional response, as the encounter rate of the predator and prey individuals increases with the prey density according to the herding index. Therefore, in the simple case of no prey handling, we obtain the herd-linear functional response. More complex is the herd-Holling type II functional response, already used by Djilali in [8] in a time delay differential model, which we derive in this paper from a system of fast time-scale state transitions.

The mechanistic derivation of the functional response follows from the time-scale separation method that has been recently formalised by Berardo et al. [11] and previously used in the works by Geritz and Gyllenberg [12,13,14]. In particular, analogously to Metz and Diekmann [15], we assume the predators in two states, searching and handling, and model the state transitions on a fast time-scale compared to other processes. The equilibrium distribution of the predators between the two states depends on the density of the prey available for capture on the edge of the herd and, as a consequence, the functional response varies with the prey density in a similar way as the Holling type II.

Finally, in this paper we study four different predator–prey–herd models. All the possible combinations arise from the following situations: first of all, assuming that the prey gather in herds, secondly, assuming that the predator can be specialist, i.e., feeding only on one species, or generalist, i.e., feeding on multiple resources, and lastly, considering two functional responses, the herd-linear and herd-Holling type II functional responses. The aim of the paper is to derive their mathematical formulation from the individual-level state transitions, and compare the models’ dynamics in terms of equilibria, stability and bifurcation diagrams.

The paper is organised as follows. In Section 2, we introduce the mechanistic derivation of the herd-Holling type II functional response for the predator–prey dynamics. In Section 3, we give partial results on the equilibrium and linear stability analysis, and the bifurcation diagrams obtained with numerical methods. In Section 4, we outline the main outcomes and draw our conclusions.

2. Materials and Methods

Mechanistic Derivation of the Functional Response for the Predator–Prey Dynamics

We model the scenario where a prey species x and a predator species y interact in the following way. The prey gather in herds and the predators are partitioned among searching individuals $y^S$ and handling individuals $y^H$. We define the handling time as the time necessary for killing, eating and digesting the prey, as well as resting and giving birth.

Suppose that there is one prey herd and multiple predator individuals. Assume further that

(i)

the geometric shape of the herd is fixed, e.g., a circle in two dimensions or a ball in three dimensions,

(ii)

the density of individuals inside the shape stays constant independently of the size of the herd,

(iii)

the predators can attack the prey on the edge of the herd or in a small boundary layer of the herd.

Therefore, the number of individuals exposed to predation is proportional to $x^{\alpha }$, where the exponent $\alpha$ is defined as herding index [1] and depends only on the characteristic shape of the herd. Examples are $\alpha =\frac{1}{2}$ for a circle in two dimensions or $\alpha =\frac{2}{3}$ for a sphere in three dimensions. Intermediate values of $\alpha$ would model fractal or multi connected sets; such exponents are investigated in [7].

Let us assume that the attack rate for a searching predator is $a{x}^{\alpha }$, i.e., the per capita attack rate is proportional to the number of available prey according to mass action. Thus, the time until prey capture becomes $\frac{1}{a{x}^{\alpha }}$. When we exclude handling, the predator functional response is linear and takes the form $f\left(x\right)=a{x}^{\alpha }$ (herd-linear functional response). Alternatively, let h denote the handling time, then the time between two successive catches by the same predator is $\frac{1}{a{x}^{\alpha }+h}$.

All in all, consider the following fast processes

$$\begin{array}{cc}\hfill & \overline{)y^S}\stackrel{a{x}^{\alpha }}{\to }\overline{)y^H}\mathit{the}\mathit{predator}\mathit{meets}\mathit{the}\mathit{prey}\mathit{and}\mathit{the}\mathit{prey}\mathit{is}\mathit{caught},\hfill \\ \hfill & \overline{)y^H}\stackrel{h^{-1}}{\to }\overline{)y^S}\mathit{the}\mathit{predator}\mathit{stops}\mathit{handling}.\hfill \end{array}$$

(1)

The above interactions were first used by Metz and Diekmann [15] to mechanistically derive the Holling type II functional response and what differs here is the exponent $\alpha$ in the encounter rate $a{x}^{\alpha }$, defining predation on the edge of the herd. We apply the time-scale separation method between fast and slow processes (for details on this modelling approach we refer to [11]). In particular, birth and death happen on a slower time scale. By converting the individual-level processes in (1) into differential equations, considering the rates at which individuals leave and enter the two subsets $y^S$ and $y^H$, the equations for the fast-time dynamics reduce to

$$\begin{array}{ccc}\hfill \frac{d{y}^S}{dt}& =& -a{x}^{\alpha }{y}^S+h^{-1}{y}^H,\hfill \end{array}$$

(2)

$$\begin{array}{ccc}\hfill \frac{d{y}^H}{dt}& =& +a{x}^{\alpha }{y}^S-h^{-1}{y}^H,\hfill \end{array}$$

(3)

where we recall that a is the predator’s attack rate and h denotes the handling time as defined in the individual-level processes above. The total population density y is constant on the fast time scale (as $y^S$ and $y^H$ verify $\frac{dy}{dt}=\frac{d{y}^S}{dt}+\frac{d{y}^H}{dt}=0$), therefore, we can reduce the system to one equation by using the conservation law $y=y^S+y^H$. By using this law and solving the equilibrium equation for $y^S$, i.e., setting the right hand side of (2) and (3) to zero, we obtain a unique quasi-steady state for the fast dynamics

$$y^S=\frac{y}{1+ah{x}^{\alpha }},y^H=y-y^S.$$

(4)

By definition, the corresponding functional response is given by the average number of prey caught by a searching predator per unit of time, i.e., $f\left(x\right)=\frac{a{x}^{\alpha }{y}^S}{y}$, and, plugging the expression for $y^S$ at the quasi-equilibrium in (4), we obtain

$$f\left(x\right)=\frac{a{x}^{\alpha }}{1+ah{x}^{\alpha }}.$$

(5)

We name the functional response in (5) the herd-Holling type II functional response. Note that when $\alpha =1$ we recover the Holling type II functional response and when $\alpha =2$ we obtain the Holling type III functional response. However, in the present context with prey herd geometry, we necessarily have $\alpha \in (0,1)$. In this case, the shape of the functional response is similar to the Holling type II functional response as it is concave and saturating, but the behaviour near the origin is different as it is infinitely steep.

3. Results

In this section we present the theoretical and numerical results for the following the Rosenzweig and MacArthur [16] models with functional responses derived in Section 2:

(i)

specialist predator and herd-linear functional response $f\left(x\right)=a{x}^{\alpha }$;

(ii)

generalist predator and herd-linear functional response $f\left(x\right)=a{x}^{\alpha }$;

(iii)

specialist predator and herd-Holling type II like functional response $f\left(x\right)=\frac{a{x}^{\alpha }}{1+ah{x}^{\alpha }}$;

(iv)

generalist predator and herd-Holling type II like functional response $f\left(x\right)=\frac{a{x}^{\alpha }}{1+ah{x}^{\alpha }}$.

The analysis is organised as follows. First, we check the unboundedness of the prey and predator populations, we derive the dimensionless version of the equations, we compute the equilibrium points and we study their stability by applying linear stability analysis. Section 3.1, Section 3.2, Section 3.3 and Section 3.4 cover this first part of the study and give the analytical results for each of the models above. Finally, in Section 3.5, we use standard numerical methods and the software Xppaut to give the one-parameter and two-parameter bifurcation diagrams.

3.1. Predator–Prey Dynamics with Specialist Predator and Herd-Linear Functional Response: Boundedness, Equilibrium Points and Stability Analysis

We study the dynamics of the model by Rosenzweig and MacArthur [16] with herd-linear response $f\left(x\right)=a{x}^{\alpha }$, conversion factor e of captured prey into new predators, per capita natural mortality rate m for the predators and logistic growth $g\left(x\right)=rx\left(1-\frac{x}{K}\right)$ for the prey, where r denotes their net growth rate and K their carrying capacity,

$$\begin{array}{ccc}\hfill \frac{dx}{dt}& =& rx\left(1-\frac{x}{K}\right)-a{x}^{\alpha }y,\hfill \end{array}$$

(6)

$$\begin{array}{ccc}\hfill \frac{dy}{dt}& =& ea{x}^{\alpha }y-my.\hfill \end{array}$$

(7)

We show that the populations do not grow unbounded (we refer to the work by Bulai and Venturino in [7]). We define with $\psi \left(t\right)=x\left(t\right)+y\left(t\right)$ the total population density and, summing up the equations for the prey and predator populations, we obtain

$$\begin{array}{c}\hfill \frac{d\psi \left(t\right)}{dt}=rx\left(1-\frac{x}{K}\right)-(1-e)a{x}^{\alpha }y-m\psi \left(t\right)+mx.\end{array}$$

(8)

We collect $\frac{d\psi \left(t\right)}{dt}+m\psi \left(t\right)$ on the left-hand side and drop the term $(1-e)a{x}^{\alpha }y>0$ to obtain

$$\frac{d\psi \left(t\right)}{dt}+m\psi \left(t\right)\le \underset{x}{max}\left\{rx\left(1-\frac{x}{K}\right)+mx\right\}.$$

(9)

The value of ${max}_x\left\{rx\left(1-\frac{x}{K}\right)+mx\right\}$ is at $x=\frac{(m+r)K}{2r}$ and, by substituting this, we obtain

$$\frac{d\psi \left(t\right)}{dt}+m\psi \left(t\right)\le \frac{K{(r+m)}^2}{4r}\equiv \overline{M}.$$

(10)

We solve the equation for $\psi \left(t\right)$ and get the upper bound for $\psi \left(t\right)$, as well as $x\left(t\right)$ and $y\left(t\right)$

$$\psi \left(t\right)=e^{-mt}\left(\psi \left(0\right)-\frac{\overline{M}}{m}\right)+\frac{\overline{M}}{m}\le max\left\{\psi \left(0\right),\frac{\overline{M}}{m}\right\}.$$

(11)

To reduce the number of parameters, we introduce the dimensionless quantities $\tilde{x}=\frac{x}{K}$, $\tilde{y}=\frac{y}{K}$, $\tilde{t}=rt$, $\tilde{a}=\frac{a{K}^{\alpha }}{r}$, $\tilde{m}=\frac{m}{r}$. Applying the substitutions and dropping the tildes, we obtain the non-dimensional system

$$\begin{array}{ccc}\hfill \frac{dx}{dt}& =& x(1-x)-a{x}^{\alpha }y,\hfill \end{array}$$

(12)

$$\begin{array}{ccc}\hfill \frac{dy}{dt}& =& ea{x}^{\alpha }y-my,\hfill \end{array}$$

(13)

with $g\left(x\right)=x(1-x)$ and $f\left(x\right)=a{x}^{\alpha }$.

The equilibria follow by setting the equations in (12) and (13) to zero. We obtain the trivial equilibrium points of the system

$$E_0=(0,0),E_1=(1,0)$$

(14)

and the interior equilibrium $E^{\ast }=(x^{\ast },y^{\ast })$ with full expression below

$$E^{\ast }=\left({\left(\frac{m}{ae}\right)}^{\frac{1}{\alpha }},\frac{1}{a}{\left(\frac{m}{ae}\right)}^{\frac{1-\alpha }{\alpha }}\left[1-{\left(\frac{m}{ae}\right)}^{\frac{1}{\alpha }}\right]\right).$$

(15)

Note that the interior equilibrium is positive if $\frac{m}{ae}<1$.

We use the Jacobian matrix of the system in (12) and (13) to study the stability of the equilibria

$$J(x,y)=\left[\begin{array}{cc}1-2x-a{x}^{\alpha -1}\alpha y& -a{x}^{\alpha }\\ ae{x}^{\alpha -1}\alpha y& ae{x}^{\alpha }-m\end{array}\right].$$

(16)

The Jacobian evaluated at $E_0$ has possibly a singularity, but the instability of this point can be assessed looking back at the original Equations (12) and (13). With $y=0$, and x near 0, the first equation behaves like $x^{\prime }\approx rx$, so that x grows. Conversely, on $x=0$ the second equation is $y^{\prime }\approx -my$ and $y\to 0$. Hence, $E_0$ is a saddle. When evaluated at the equilibrium $E_1$, the determinant of the Jacobian matrix is $m-ae$ and is positive if $\frac{m}{ae}>1$, that is, when the interior equilibrium is not feasible (i.e., does not exist or is negative). Under the same condition, the trace of the Jacobian at $E_1$, $ae-m-1$, is negative and the equilibrium is stable.

For simplicity, we rewrite the Jacobian matrix evaluated at the interior equilibrium $E_{\ast }=(x^{\ast },y^{\ast })$ in terms of the functions $f\left(x\right)$ and $g\left(x\right)$

$$J(x^{\ast },y^{\ast })=\left[\begin{array}{cc}{g}^{\prime }\left(x^{\ast }\right)-f^{\prime }\left(x^{\ast }\right)y^{\ast }& -f\left(x^{\ast }\right)\\ e{f}^{\prime }\left(x^{\ast }\right)y^{\ast }& ef\left(x^{\ast }\right)-m\end{array}\right]=\left[\begin{array}{cc}f\left(x^{\ast }\right){\left(\frac{g\left(x^{\ast }\right)}{f\left(x^{\ast }\right)}\right)}^{\prime }& -\frac{m}{e}\\ e{f}^{\prime }\left(x^{\ast }\right)y^{\ast }& 0\end{array}\right].$$

(17)

The determinant of the matrix in (17) is $m{f}^{\prime }\left(x^{\ast }\right)y^{\ast }>0$ (since the functional response is an increasing function of the prey density), therefore the stability of the interior equilibrium depends on the sign of the trace $f\left(x^{\ast }\right){\left(\frac{g\left(x^{\ast }\right)}{f\left(x^{\ast }\right)}\right)}^{\prime }$, and, in particular, on the slope of the prey zero-growth isocline ${\left(\frac{g\left(x^{\ast }\right)}{f\left(x^{\ast }\right)}\right)}^{\prime }=\frac{x^{-\alpha }[x(\alpha -2)+1-\alpha ]}{a}{|}_{x=x^{\ast }}$ (see also Gause [17] and Gause et al. [18]). We obtain that the equilibrium $E^{\ast }$ is stable if $\frac{m}{ae}>{\left(\frac{\alpha -1}{\alpha -2}\right)}^{\alpha }\in (0,1)$ (check

Table 1. Conditions for feasibility and stability of the equilibria of the system in (12) and (13). The trivial equilibrium $E_0$ is always unstable. TB: transcritical bifurcation. HB: Hopf bifurcation.

Condition	${\mathit{E}}_1$	${\mathit{E}}_{\ast }$	Bifurcation
$\frac{m}{ae}>1$	asympt. stable	not feasible
$\frac{m}{ae}=1$			TB
${\left(\frac{\alpha -1}{\alpha -2}\right)}^{\alpha }<\frac{m}{ae}<1$	unstable	asympt. stable
$\frac{m}{ae}={\left(\frac{\alpha -1}{\alpha -2}\right)}^{\alpha }$			HB
$0<\frac{m}{ae}<{\left(\frac{\alpha -1}{\alpha -2}\right)}^{\alpha }$	unstable	unstable

for a summary of the feasibility and stability conditions).

We conclude that a transcritical bifurcation occurs at $\frac{m}{ae}=1$, where the interior equilibrium exchanges stability with the predator-free equilibrium. A Hopf bifurcation appears at $\frac{m}{ae}={\left(\frac{\alpha -1}{\alpha -2}\right)}^{\alpha }$, as the eigenvalues of the community matrix become purely imaginary and the system converges to a stable limit cycle.

3.2. Predator–Prey Dynamics with Generalist Predator and Herd-Linear Functional Response: Boundedness, Equilibrium Points and Stability Analysis

We study the dynamics of the model by Rosenzweig and MacArthur [16] with herd-linear response $f\left(x\right)=a{x}^{\alpha }$, conversion factor e of captured prey into new predators and per capita natural mortality rate m for the predators. We assume logistic growth for both the prey and the predator species, $g_x\left(x\right)=rx\left(1-\frac{x}{K_x}\right)$ and $g_y\left(x\right)=sx\left(1-\frac{x}{K_y}\right)$, respectively, where r is the net growth rate of the prey and $K_x$ their carrying capacity, while s denotes the predators’ reproduction rate, i.e., not discounted by deaths,

$$\begin{array}{ccc}\hfill \frac{dx}{dt}& =& rx\left(1-\frac{x}{K_x}\right)-a{x}^{\alpha }y,\hfill \end{array}$$

(18)

$$\begin{array}{ccc}\hfill \frac{dy}{dt}& =& sy\left(1-\frac{y}{K_y}\right)+ea{x}^{\alpha }y-my.\hfill \end{array}$$

(19)

In this way, the predators are subjected to intraspecific competition, which occurs at rate $\frac{s}{K_y}$.

To check that the populations do not grow unbounded, we set $\psi \left(t\right)=x\left(t\right)+y\left(t\right)$ and, by repeating the steps in Section 3.1, we get the differential equation for the total population

$$\frac{d\psi \left(t\right)}{dt}+m\psi \left(t\right)\le \underset{x,y}{max}\left\{rx\left(1-\frac{x}{K_x}\right)+sy\left(1-\frac{y}{K_y}\right)+mx\right\}.$$

(20)

We differentiate the right-hand term with respect to x and y to get the local maximum $\left(\frac{(r+m)}{2r}{K}_x,\frac{K_y}{2}\right)$. By substitution in the equation above, we obtain

$$\frac{d\psi \left(t\right)}{dt}+m\psi \left(t\right)\le \frac{{(r+m)}^2}{4r}{K}_x+\frac{s}{4}{K}_y\equiv \overline{M}.$$

(21)

Therefore, the solution for the total population reads

$$\psi \left(t\right)=e^{-mt}\left(\psi \left(0\right)-\frac{\overline{M}}{m}\right)+\frac{\overline{M}}{m}\le max\left\{\psi \left(0\right),\frac{\overline{M}}{m}\right\},$$

(22)

where the upper bound is applicable also for $x\left(t\right)$ and $y\left(t\right)$.

To obtain the non-dimensional version of the system in (18) and (19), we consider the dimensionless variables and parameters $\tilde{x}=\frac{x}{K_x}$, $\tilde{y}=\frac{y}{K_y}$, $\tilde{t}=rt$, $\tilde{a}=\frac{K_y}{K_x^{1-\alpha }r}a$, $\tilde{s}=\frac{s}{r}$, $\tilde{e}=\frac{K_x}{K_y}e$, $\tilde{m}=\frac{m}{r}$. We drop the tildes and obtain the dimensionless system

$$\begin{array}{ccc}\hfill \frac{dx}{dt}& =& x(1-x)-a{x}^{\alpha }y,\hfill \end{array}$$

(23)

$$\begin{array}{ccc}\hfill \frac{dy}{dt}& =& sy(1-y)+ea{x}^{\alpha }y-my,\hfill \end{array}$$

(24)

with $g_x\left(x\right)=x(1-x)$, $g_y\left(y\right)=sy(1-y)$ and $f\left(x\right)=a{x}^{\alpha }$.

We proceed with computing the equilibria. The trivial equilibria are

$$E_0=(0,0),E_1=(1,0),E_2=\left(0,1-\frac{m}{s}\right),$$

(25)

with $E_2$ feasible if $m<s$. The interior equilibria are given by the intersection of the isoclines

$$y=\frac{(1-x)x^{1-\alpha }}{a}$$

(26)

and

$$y=1-\frac{m}{s}+\frac{ea}{s}{x}^{\alpha }.$$

(27)

Note that the isocline in (26) intersects the x-axis at $(0,0)$ and $(1,0)$ and has a maximum at $x=\frac{1-\alpha }{2-\alpha }<\frac{1}{2}$ for $0<\alpha <1$, while the isocline in (27) intersects the vertical axis at $(0,1-\frac{m}{s})$ and is a root function translated by $1-\frac{m}{s}$ and dilated by $\frac{ea}{s}$. Therefore, if the intersection point of the isocline in (27) lies in the positive quadrant, i.e., if $m<s$, we find three different configurations for the phase plane: the two isoclines can intersect at most twice at $E_{\ast 1}$ and $E_{\ast 2}$ if $1-\frac{m}{s}<\frac{1}{a(2-\alpha )}{\left(\frac{1-\alpha }{2-\alpha }\right)}^{1-\alpha }-\frac{ea}{s}{\left(\frac{1-\alpha }{2-\alpha }\right)}^{\alpha }$, or be tangent at $E_{\ast }=\left(\frac{1-\alpha }{2-\alpha },1-\frac{m}{s}+\frac{ea}{s}{\left(\frac{1-\alpha }{2-\alpha }\right)}^{\alpha }\right)$ when $1-\frac{m}{s}=\frac{1}{a(2-\alpha )}{\left(\frac{1-\alpha }{2-\alpha }\right)}^{1-\alpha }-\frac{ea}{s}{\left(\frac{1-\alpha }{2-\alpha }\right)}^{\alpha }$ or never intersect in $x\in (0,1)$ when $1-\frac{m}{s}>\frac{1}{a(2-\alpha )}{\left(\frac{1-\alpha }{2-\alpha }\right)}^{1-\alpha }-\frac{ea}{s}{\left(\frac{1-\alpha }{2-\alpha }\right)}^{\alpha }$. The equilibria are obtained as the positive roots of the curve

$$\varphi \left(x\right)=1-\frac{m}{s}+\frac{ea}{s}{x}^{\alpha }-(1-x)\frac{x^{1-\alpha }}{a}$$

(28)

and the non-negativity of y is ensured by the condition

$${\left(\frac{m-1}{ea}\right)}^{\frac{1}{\alpha }}<x<1.$$

(29)

When $m>s$, the isocline in (27) intersects the vertical axes at $y=1-\frac{m}{s}<0$ and we find at most one interior equilibrium $E_{\ast }$. We obtain the feasibility condition for $E_{\ast }$ by imposing that the curve in (27) takes positive values at $x=1$, that is, if $m<s+ea$.

The Jacobian matrix of the system in (23) and (24) is given by

$$J(x,y)=\left[\begin{array}{cc}1-2x-a{x}^{\alpha -1}\alpha y& -a{x}^{\alpha }\\ ae{x}^{\alpha -1}\alpha y& s-2sy+ae{x}^{\alpha }-m\end{array}\right].$$

(30)

The equilibrium $E_0$, restricting the analysis to the trajectories on the coordinate axes, is seen to be either an unstable node if $m<s$, or a saddle if $m>s$. The prey-only equilibrium $E_1$ is a stable node if $m>s+ae$, otherwise the steady state is an unstable saddle. Under its feasibility condition $m<s$, the equilibrium $E_2$ is always an unstable saddle. We summarise the feasibility and stability conditions studied above in

Table 2. Conditions for the feasibility and stability of the equilibria of the system in (23) and (24).

Condition	${\mathit{E}}_0$	${\mathit{E}}_1$	${\mathit{E}}_2$	${\mathit{E}}_{\ast 1}$	${\mathit{E}}_{\ast 2}$
$m<s$	unstable	unstable	unstable	See Table 3	See Table 3
$s<m<s+ae$	unstable	unstable	not feas.	stable	not feasible
$m>s+ae$	unstable	asympt. stable	not feasible	not feasible	not feasible

and

Table 3. Conditions for feasibility of the interior equilibria of the system in (23) and (24) when $m<s$.

Condition	${\mathit{E}}_{\ast 1}$	${\mathit{E}}_{\ast 2}$
$1-\frac{m}{s}<\frac{1}{a(2-\alpha )}{\left(\frac{1-\alpha }{2-\alpha }\right)}^{1-\alpha }-\frac{ea}{s}{\left(\frac{1-\alpha }{2-\alpha }\right)}^{\alpha }$	feasible	feasible
$1-\frac{m}{s}=\frac{1}{a(2-\alpha )}{\left(\frac{1-\alpha }{2-\alpha }\right)}^{1-\alpha }-\frac{ea}{s}{\left(\frac{1-\alpha }{2-\alpha }\right)}^{\alpha }$	feasible	not feasible
$1-\frac{m}{s}>\frac{1}{a(2-\alpha )}{\left(\frac{1-\alpha }{2-\alpha }\right)}^{1-\alpha }-\frac{ea}{s}{\left(\frac{1-\alpha }{2-\alpha }\right)}^{\alpha }$	not feasible	not feasible

.

We rewrite the Jacobian matrix evaluated at the interior equilibrium in terms of functions $f\left(x\right)$, $g_x\left(x\right)$ and $g_y\left(y\right)$:

$$J(x^{\ast },y^{\ast })=\left[\begin{array}{cc}{g}_x^{\prime }\left(x^{\ast }\right)-f^{\prime }\left(x^{\ast }\right)y^{\ast }& -f\left(x^{\ast }\right)\\ e{f}^{\prime }\left(x^{\ast }\right)y^{\ast }& g_y^{\prime }\left(y\right)+ef\left(x^{\ast }\right)-m\end{array}\right]=\left[\begin{array}{cc}f\left(x^{\ast }\right){\left(\frac{g_x\left(x^{\ast }\right)}{f\left(x^{\ast }\right)}\right)}^{\prime }& -f\left(x^{\ast }\right)\\ e{f}^{\prime }\left(x^{\ast }\right)y^{\ast }& -y^{\ast }\end{array}\right].$$

(31)

The trace and the determinant at the interior equilibrium are given by

$$\begin{array}{ccc}\hfill \mathrm{tr}J(x^{\ast },y^{\ast })& =& f\left(x^{\ast }\right){\left(\frac{g_x\left(x^{\ast }\right)}{f\left(x^{\ast }\right)}\right)}^{\prime }-y^{\ast },\hfill \end{array}$$

(32)

$$\begin{array}{ccc}\hfill detJ(x^{\ast },y^{\ast })& =& -y^{\ast }f\left(x^{\ast }\right){\left(\frac{g_x\left(x^{\ast }\right)}{f\left(x^{\ast }\right)}\right)}^{\prime }+ef\left(x^{\ast }\right)f^{\prime }\left(x^{\ast }\right)y^{\ast }.\hfill \end{array}$$

(33)

When only one interior equilibrium exists and is positive, the signs of the trace and the determinant determine its asymptotic stability, more specifically if $\mathrm{tr}J(x^{\ast },y^{\ast })<0$ and $detJ(x^{\ast },y^{\ast })>0$.

It seems rather difficult to obtain more analytical details on the stability of the equilibria and bifurcations for the model in (23) and (24). If possible, a more detailed analysis will be the topic of a future work.

3.3. Predator–Prey Dynamics with Specialist Predator and Herd-Holling Type II Functional Response: Boundedness, Equilibrium Points and Stability Analysis

We study the dynamics of the model by Rosenzweig and MacArthur [16] with the herd-Holling type II functional response $f\left(x\right)=\frac{a{x}^{\alpha }}{1+ah{x}^{\alpha }}$ derived in Section 2, conversion factor e of captured prey into new predators and per capita natural mortality rate m for the predators, logistic growth $g\left(x\right)=rx\left(1-\frac{x}{K}\right)$ for the prey, where r denotes their net growth rate and K their carrying capacity,

$$\begin{array}{ccc}\hfill \frac{dx}{dt}& =& rx\left(1-\frac{x}{K}\right)-\frac{a{x}^{\alpha }}{1+ah{x}^{\alpha }}y,\hfill \end{array}$$

(34)

$$\begin{array}{ccc}\hfill \frac{dy}{dt}& =& e\frac{a{x}^{\alpha }}{1+ah{x}^{\alpha }}y-my.\hfill \end{array}$$

(35)

The total population $\psi \left(t\right)=x\left(t\right)+y\left(t\right)$ verifies

$$\psi \left(t\right)\le max\left\{\psi \left(0\right),\frac{\overline{M}}{m}\right\},$$

(36)

with $\overline{M}=\frac{K{(r+m)}^2}{4r}$ as in Section 3.1.

We obtain the dimensionless version of the model by applying the substitutions $\tilde{x}=\frac{x}{K}$, $\tilde{y}=\frac{y}{K}$, $\tilde{t}=rt$, $\tilde{a}=\frac{a{K}^{\alpha }}{r}$, $\tilde{h}=rh$, $\tilde{m}=\frac{m}{r}$. We drop the tildes and get the equations

$$\begin{array}{ccc}\hfill \frac{dx}{dt}& =& x\left(1-x\right)-\frac{a{x}^{\alpha }}{1+ah{x}^{\alpha }}y,\hfill \end{array}$$

(37)

$$\begin{array}{ccc}\hfill \frac{dy}{dt}& =& e\frac{a{x}^{\alpha }}{1+ah{x}^{\alpha }}y-my,\hfill \end{array}$$

(38)

with $g\left(x\right)=x\left(1-x\right)$ and $f\left(x\right)=\frac{a{x}^{\alpha }}{1+ah{x}^{\alpha }}$.

We compute the equilibria by setting the equations in (37) and (38) to zero. The trivial equilibria are $E_0=(0,0)$ and $E_1=(1,0)$, while the unique interior equilibrium $E_{\ast }=(x^{\ast },y^{\ast })$ has full expression

$$E_{\ast }=\left({\left(\frac{m}{ea-mah}\right)}^{\frac{1}{\alpha }},\frac{e}{m}{x}^{\ast }(1-x^{\ast })\right)$$

(39)

and exists and is positive if and only if $mah+m<ea$.

The Jacobian matrix corresponding to the system in (37) and (38) is given by

$$J(x,y)=\left[\begin{array}{cc}1-2x-\frac{a\alpha x^{\alpha -1}}{{(1+ah{x}^{\alpha })}^2}y& \frac{a{x}^{\alpha }}{1+ah{x}^{\alpha }}\\ e\frac{a\alpha x^{\alpha -1}}{{(1+ah{x}^{\alpha })}^2}y& e\frac{a{x}^{\alpha }}{1+ah{x}^{\alpha }}-m\end{array}\right].$$

(40)

The origin is unstable, being a saddle, a fact that is shown restricting the system to the coordinate axes, as previously done for the system (12) and (13). The equilibrium $E_1$ is stable if and only if $ea<mah+m$ (under this condition the determinant of the Jacobian matrix at the equilibrium is positive and the trace is negative). The prey-only equilibrium changes its stability at $ea=mah+m$ when a transcritical bifurcation occurs. We can use the same formulation as in (17) for the Jacobian evaluated at the interior equilibrium, which for convenience we reproduce here

$$J(x^{\ast },y^{\ast })=\left[\begin{array}{cc}{g}^{\prime }\left(x^{\ast }\right)-f^{\prime }\left(x^{\ast }\right)y^{\ast }& -f\left(x^{\ast }\right)\\ e{f}^{\prime }\left(x^{\ast }\right)y^{\ast }& ef\left(x^{\ast }\right)-m\end{array}\right]=\left[\begin{array}{cc}f\left(x^{\ast }\right){\left(\frac{g\left(x^{\ast }\right)}{f\left(x^{\ast }\right)}\right)}^{\prime }& -\frac{m}{e}\\ e{f}^{\prime }\left(x^{\ast }\right)y^{\ast }& 0\end{array}\right].$$

(41)

As the determinant of the Jacobian matrix is always positive, the stability of the interior equilibrium depends on the sign of the trace, in particular on the slope of the predator isocline, ${\left(\frac{g\left(x^{\ast }\right)}{f\left(x^{\ast }\right)}\right)}^{\prime }=\frac{x^{-\alpha }[x(\alpha -2)+1-\alpha +ah{x}^{\alpha }(1-2x)]}{a}{|}_{x=x^{\ast }}$. For the same reason as in Section 3.1, the system in (37) and (38) undergoes a Hopf bifurcation when ${\left(\frac{g\left(x^{\ast }\right)}{f\left(x^{\ast }\right)}\right)}^{\prime }=0$ and converges to a stable limit cycle for ${\left(\frac{g\left(x^{\ast }\right)}{f\left(x^{\ast }\right)}\right)}^{\prime }>0$, otherwise it converges to the interior equilibrium $E_{\ast }$. In

Table 4. Conditions for feasibility and stability of the equilibria of the system in (37) and (38). TB: transcritical bifurcation. HB: Hopf bifurcation.

Condition	${\mathit{E}}_0$	${\mathit{E}}_1$	${\mathit{E}}_{\ast }$	Bifurcation
$ea<mah+m$	unstable	asympt. stable	not feasible
$ea=mah+m$				TB
$ea>mah+m$, ${\left(\frac{g\left(x^{\ast }\right)}{f\left(x^{\ast }\right)}\right)}^{\prime }<0$	unstable	unstable	asympt. stable
$ea>mah+m$, ${\left(\frac{g\left(x^{\ast }\right)}{f\left(x^{\ast }\right)}\right)}^{\prime }=0$				HB
$ea>mah+m$, ${\left(\frac{g\left(x^{\ast }\right)}{f\left(x^{\ast }\right)}\right)}^{\prime }>0$	unstable	unstable	unstable

we give a summary of the feasibility and stability conditions of the equilibria.

3.4. Predator–Prey Dynamics with Generalist Predator and Herd-Holling Type II Functional Response: Boundedness, Equilibrium Points and Stability Analysis

We study the dynamics of the model by Rosenzweig and MacArthur [16] with the herd-Holling type II functional response $f\left(x\right)=\frac{a{x}^{\alpha }}{1+ah{x}^{\alpha }}$ derived in Section 2, conversion factor e of captured prey into new predators and per capita natural mortality rate m for the predators. We assume logistic growth for both the prey and the predator species, $g_x\left(x\right)=rx\left(1-\frac{x}{K_x}\right)$ and $g_y\left(x\right)=sx\left(1-\frac{x}{K_y}\right)$, respectively, where r is the net growth rate of the prey and $K_x$ their carrying capacity, while s denotes the predators’ reproduction rate,

$$\begin{array}{ccc}\hfill \frac{dx}{dt}& =& rx\left(1-\frac{x}{K_x}\right)-\frac{a{x}^{\alpha }}{1+ah{x}^{\alpha }}y,\hfill \end{array}$$

(42)

$$\begin{array}{ccc}\hfill \frac{dy}{dt}& =& sy\left(1-\frac{y}{K_y}\right)+e\frac{a{x}^{\alpha }}{1+ah{x}^{\alpha }}y-my.\hfill \end{array}$$

(43)

Once again, note the second term in the predators’ equation, whose coefficient $\frac{s}{K_y}$ models predators intraspecific competition.

The boundedness of the populations is ensured by the condition on the total population density $\psi \left(t\right)=x\left(t\right)+y\left(t\right)$

$$\psi \left(t\right)\le max\left\{\psi \left(0\right),\frac{\overline{M}}{m}\right\},$$

(44)

with $\overline{M}=\frac{{(r+m)}^2}{4r}{K}_x+\frac{s}{4}{K}_y$ as in Section 3.1.

We use the dimensionless quantities $\tilde{x}=\frac{x}{K_x}$, $\tilde{y}=\frac{y}{K_y}$, $\tilde{t}=rt$, $\tilde{a}=\frac{K_y}{K_x^{1-\alpha }r}a$, $\tilde{h}=r\frac{K_x}{K_y}h$, $\tilde{s}=\frac{s}{r}$, $\tilde{e}=\frac{K_x}{K_y}e$, $\tilde{m}=\frac{m}{r}$ to obtain the dimensionless system of equations

$$\begin{array}{ccc}\hfill \frac{dx}{dt}& =& x\left(1-x\right)-\frac{a{x}^{\alpha }}{1+ah{x}^{\alpha }}y,\hfill \end{array}$$

(45)

$$\begin{array}{ccc}\hfill \frac{dy}{dt}& =& sy(1-y)+e\frac{a{x}^{\alpha }}{1+ah{x}^{\alpha }}y-my,\hfill \end{array}$$

(46)

with $g_x\left(x\right)=x\left(1-x\right)$, $g_y\left(y\right)=sy(1-y)$ and $f\left(x\right)=\frac{a{x}^{\alpha }}{1+ah{x}^{\alpha }}$.

The corresponding trivial equilibria correspond to the ones in Section 3.2 and are given by

$$E_0=(0,0),E_1=(1,0),E_2=\left(0,1-\frac{m}{s}\right),$$

(47)

with $E_2$ being feasible if $m<s$. The interior equilibria are given by the intersection of the isoclines

$$y=\frac{x(1-x)(1+ah{x}^{\alpha })}{a{x}^{\alpha }}$$

(48)

and

$$y=1-\frac{m}{s}+\frac{ea{x}^{\alpha }}{s(1+ah{x}^{\alpha })}.$$

(49)

The isocline in (48) intersects the horizontal axis at $(0,0)$ and $(1,0)$, while the isocline in (49) intersects the vertical axis at $(0,1-\frac{m}{s})$. As for the model with generalist predator and linear functional response in Section 3.2, we may expect that, under some conditions, the system admits two interior equilibria. However, given the formulation of the isoclines in (48) and (49), it seems difficult to find explicit analytical results and we refer to the next Section 3.5 for more details on the interior equilibria and their feasibility and stability conditions.

We obtain the Jacobian matrix to check the stability of the trivial equilibria:

$$J(x,y)=\left[\begin{array}{cc}1-2x-\frac{a\alpha x^{\alpha -1}}{{(1+ah{x}^{\alpha })}^2}y& \frac{a{x}^{\alpha }}{1+ah{x}^{\alpha }}\\ e\frac{a\alpha x^{\alpha -1}}{{(1+ah{x}^{\alpha })}^2}y& s-2sy+e\frac{a{x}^{\alpha }}{1+ah{x}^{\alpha }}-m\end{array}\right].$$

(50)

Again, restricting the analysis to the trajectories on the coordinate axes, we obtain that $E_0$ is an unstable node for $m<s$, or a saddle if $m>s$. We find that the origin $E_0$ is an unstable saddle, as well as the predator-only equilibrium $E_2$. The prey-only equilibrium $E_1$ is a stable node if $m>s+\frac{ea}{1+ah}$, an unstable saddle otherwise. We give a summary of these results in

Table 5. Conditions for feasibility and stability of the trivial equilibria of the system in (45) and (46).

Condition	${\mathit{E}}_0$	${\mathit{E}}_1$	${\mathit{E}}_2$
$m<s+\frac{ea}{1+ah}$	unstable	unstable	unstable
$m>s+\frac{ea}{1+ah}$, ${\left(\frac{g\left(x^{\ast }\right)}{f\left(x^{\ast }\right)}\right)}^{\prime }<0$	unstable	asympt. stable	unstable

.

3.5. One-Parameter and Two-Parameter Bifurcation Diagrams

In this section, we proceed with the bifurcation analysis. We give the one-parameter bifurcation diagrams and vary either the value of the predator mortality rate, m, or the herding index, $\alpha$, when possible. Additionally, we obtain the two-parameter bifurcation diagrams with respect to the parameter pairs $(m,🟉)$ or $(\alpha ,🟉)$, where 🟉 equals one of the other parameters in the model. Note that in the numerical simulations we use the model and parameters prior to non-dimensionalisation, to obtain a complete analysis with respect to all the model parameters.

We first study the predator–prey model with specialist predator and herd-linear functional response. In the one-parameter bifurcation plots, we fix the parameter values of the model in (6) and (7) as in

Table 6. Nominal set of parameter values for the model in (6) and (7).

Parameter	Description	Value
r	prey net reproduction rate	1
K	prey carrying capacity	5
a	predation rate	1
$\alpha$	herd exponent	$0.7$
e	conversion factor	$0.5$
m	natural mortality rate (predators)	1

and we call it the nominal set (hypothetical values).

In Figure 1 we give the one-parameter bifurcation diagrams with respect to the natural mortality rate of the predators, m. Note that when $m=0.5526$, the system undergoes a supercritical Hopf bifurcation (HB) and a stable limit cycle appears. At $m=1.543$, a transcritical bifurcation occurs, where the coexistence equilibrium loses its stability and the predator-free equilibrium becomes stable.

To complete the analysis, in Figure 2 we have plotted all the possible two-parameter bifurcation diagrams for $(m,🟉)$, with $🟉=a$, K, a, $\alpha$ or e. Note that the HB curve appears in all two-parameter bifurcation diagrams, but only in the plots where we vary $(m,K)$, $(m,a)$ and $(m,e)$ it occurs for every value of m. When we let r vary, the HB curve is present only at $m=0.5526$ (Figure 2, top left); similarly, when we allow $\alpha$ to change, the HB occurs only for $m\le 0.5$ (Figure 2, bottom left).

As a second example, we analyse the predator–prey dynamics with generalist predator and herd-linear functional response. In

Table 7. Nominal set of parameter values for the model in (18) and (19).

Parameter	Description	Value
r	prey net reproduction rate	$0.5$
$K_x$	prey carrying capacity	5
a	predation rate	1
$\alpha$	herd exponent	$0.7$
s	predator reproduction rate	$0.5$
$K_y$	predator carrying capacity	5
e	conversion factor	$0.5$
m	natural mortality rate (predators)	1

we give the nominal set of parameter values for the model in (18) and (19).

In Figure 3 we give the one-parameter bifurcation diagram with respect to the herd exponent, $\alpha$. Here a supercritical HB occurs at $\alpha =0.5995$ and a subcritical HB appears at $\alpha =0.1476$.

The two-parameters bifurcation diagrams with respect to $(\alpha ,🟉)$, with $🟉=r$, $K_x$, a, s, $K_y$, e, and m for the model with Equations (18) and (19) are given in Figure 4. When we vary the parameter pair $(\alpha ,a)$, a HB appears for all values of $\alpha$ independently of the value of a, while for the remaining cases the HB is present only for some parameter values. Moreover, we observe that only when we vary the parameter pair $(\alpha ,m)$, the transcritical bifurcation curve is present. Finally, if one of the parameters $🟉=K_x$, a, s, $K_y$, and e are below certain threshold values, with the other parameter values fixed as in Figure 4, the coexistence equilibrium is asymptotically stable; analogously, when r is above a certain threshold value the system converges to the interior equilibrium.

In this paragraph, we focus on the predator–prey dynamics with specialist predator and herd-Holling type II functional response. In

Table 8. Nominal set of parameter values for the model in (34) and (35).

Parameter	Description	Value
r	prey net reproduction rate	$0.5$
K	prey carrying capacity	5
a	predation rate	1
$\alpha$	herd exponent	$0.7$
h	handling time	$0.2$
e	conversion factor	$0.5$
m	natural mortality rate (predators)	1

we list the parameter values for the model in (34) and (35).

In Figure 5 we obtain qualitatively similar results as for (6) and (7) for the one-parameter bifurcation diagrams with respect to the natural mortality rate of the predators, m.

The dynamic described in Figure 6 is similar to the one in Figure 2. There are two main differences: in the two-parameter bifurcation diagrams with respect to $(m,K)$ and $(m,a)$ the HB curve is more concave; when we vary the parameter pair $(m,h)$ (this case is not present in Figure 2), one can see that the HB occurs for all values of h and for m smaller than a threshold value.

Finally, we study the predator–prey dynamics with generalist predator and herd-Holling type II functional response. In

Table 9. Nominal set of parameter values for the third model (42) and (43).

Parameter	Description	Value
r	prey net reproduction rate	$0.5$
$K_x$	prey carrying capacity	5
a	predation rate	1
$\alpha$	herd exponent	$0.7$
h	handling time	$0.2$
s	predator reproduction rate	$0.5$
$K_y$	predator carrying capacity	5
e	conversion factor	$0.5$
m	natural mortality rate (predators)	1

we list the nominal set of parameter values for model in (42) and (43). This is the most general model that we study which encompasses all the previously considered cases.

Both the one-parameter bifurcation diagram with respect to m, and two-parameter bifurcation diagrams with respect to $(m,🟉)$, with $🟉=r$, $K_x$, a, $\alpha$, h, s, $K_y$ and e show behaviours similar to the previous models, see Figure 7 and Figure 8, respectively. It is worth noting that the results for models (34) and (35) are different as we give the two-parameter bifurcation diagrams with respect to the herd exponent as first parameter, while we refer to the predator natural mortality rate in the other two-parameter bifurcation plots.

4. Conclusions

The aim of this paper is to formalise, by means of illustrative examples, how prey herd behaviour can be modelled with ordinary differential equations from first principles. We propose four different Rosenzweig–MacArthur predator–prey models where the prey gather in herds and only the individuals on the edge are subjected to the predators’ attacks.

We use the mechanistic approach and time-scale separation method to derive from the individual mechanisms a Holling type II–like functional response for the predators, the herd-Holling type II functional response, which takes into account the prey herd shape. As strongly emphasised in [11,12,13,14,15], the strength of the mechanistic approach lies in the possibility of interpreting the population equations in terms of individual state transition parameter rates.

Moreover, by introducing the so-called herding index $\alpha$ (see [5,6,7,8,9,10]), we model the density of the prey x on the edge of the herd as $x^{\alpha }$ and we are able to track information on the herd shape in a simple system of ordinary differential equations for the population dynamics.

We focus our attention on the ecological dynamics with one prey and one predator species, including specialist versus generalist predators, i.e., no predators’ intraspecific competition versus predators’ intraspecific competition. Such antagonistic behaviour has been widely observed in population ecology, especially in aquatic species and insects, and has been proven to deeply affect niche expansion and speciation (see, for example, [19,20,21,22,23]). We further assume either herd-linear or herd-Holling type II functional responses for the predators.

Unlike the two simple cases with specialist predator where both the analytical and numerical results on the attractors and bifurcations for the population dynamics given in this paper are exhaustive, for the two models assuming intraspecific competition there is room for a more detailed study. Indeed, we have found that such ecological dynamics can give rise to multiple stable attractors (cycles and steady states), while our current analysis comprises only the cases with one interior equilibrium.

Moreover, given the evidence of possible chaotic behaviour in models with prey herd behaviour, as found for instance in the recent paper [24], the analysis can be widened in both analytical and numerical directions, by introducing an adequate change of variable. This has the advantage of simplifying the exponential term in the population equations as in [7] and analysing two-bifurcation diagrams in the same way as performed in [25].

Further extensions of our models may include more complex types of predators and prey-predator dynamics, e.g., generalist predators with respect to prey switching and prey–predator–superpredator models, where the superpredator species feeds on the prey and alternative food sources, such as cannibalising the weaker and younger predator individuals (see, for example, the state transitions and time-scale separation in [26,27] for a system with one predator species with cannibalistic tendencies). However, the greater the number of states and state transitions, the more complicated the population dynamics. Extending the mechanistic derivation of the population equations from individual-level interactions to ecosystems with multiple prey and predator species can further increase the complexity of the model, which may become parameter-heavy and will require advanced numerical methods.