**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:


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. Sections 3.1–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 herdlinear response *f*(*x*) = *axα*, conversion factor *e* of captured prey into new predators, per capita natural mortality rate *m* for the predators and logistic growth *g*(*x*) = *rx*- 1 − *x K* for the prey, where *r* denotes their net growth rate and *K* their carrying capacity,

$$\frac{d\mathbf{x}}{dt}\_{\mathbf{y}} = -r\mathbf{x}\left(1 - \frac{\mathbf{x}}{K}\right) - a\mathbf{x}^a y\_{\mathbf{y}} \tag{6}$$

$$\frac{dy}{dt} \quad = \quad eax^{\kappa}y - my. \tag{7}$$

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

$$\frac{d\psi(t)}{dt} = r\chi\left(1 - \frac{\chi}{K}\right) - (1 - \varepsilon)ax^{a}y - m\psi(t) + m\chi. \tag{8}$$

We collect *dψ*(*t*) *dt* + *mψ*(*t*) on the left-hand side and drop the term (1 − *e*)*ax<sup>α</sup>y* > 0 to obtain

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

The value of max*xrx*-1 − *xK* + *mx* is at *x* = (*m*+*r*)*<sup>K</sup>* 2*r* and, by substituting this, we obtain

$$\frac{d\psi(t)}{dt} + m\psi(t) \le \frac{K(r+m)^2}{4r} \equiv \bar{M}.\tag{10}$$

We solve the equation for *ψ*(*t*) and ge<sup>t</sup> the upper bound for *ψ*(*t*), as well as *x*(*t*) and *y*(*t*)

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

To reduce the number of parameters, we introduce the dimensionless quantities *x*˜ = *xK* , *y* ˜ = *yK* , ˜*t* = *rt*, *a*˜ = *aKα r* , *m*˜ = *mr* . Applying the substitutions and dropping the tildes, we obtain the non-dimensional system

$$\frac{dx}{dt}\_{\text{m}} = \left. x(1-x) - ax^{a}y \right|\_{\text{m}} \tag{12}$$

$$\frac{dy}{dt} \quad = \quad \epsilon \mathbf{a} \mathbf{x}^a y - my,\tag{13}$$

with *g*(*x*) = *x*(<sup>1</sup> − *x*) and *f*(*x*) = *axα*.

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), \quad E\_1 = (1,0) \tag{14}$$

and the interior equilibrium *E*∗ = (*x*<sup>∗</sup>, *y*<sup>∗</sup>) with full expression below

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

Note that the interior equilibrium is positive if *mae*< 1.

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

$$J(x,y) = \begin{bmatrix} 1 - 2x - ax^{a-1}ay & -ax^{a} \\ aex^{a-1}ay & aex^{a} - m \end{bmatrix} . \tag{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* ≈ *rx*, so that *x* grows. Conversely, on *x* = 0 the second equation is *y* ≈ −*my* and *y* → 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 *mae* > 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*∗ = (*x*<sup>∗</sup>, *y*<sup>∗</sup>) in terms of the functions *f*(*x*) and *g*(*x*)

$$J(\mathbf{x}^\*, y^\*) = \begin{bmatrix} y'(\mathbf{x}^\*) - f'(\mathbf{x}^\*)y^\* & -f(\mathbf{x}^\*) \\ ef'(\mathbf{x}^\*)y^\* & ef(\mathbf{x}^\*) - m \end{bmatrix} = \begin{bmatrix} f(\mathbf{x}^\*) \begin{pmatrix} \frac{g'(\mathbf{x}^\*)}{f(\mathbf{x}^\*)} \end{pmatrix}' & -\frac{m}{\epsilon} \\ ef'(\mathbf{x}^\*)y^\* & 0 \end{bmatrix}. \tag{17}$$

The determinant of the matrix in (17) is *m f* (*x*<sup>∗</sup>)*y*<sup>∗</sup> > 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 *<sup>f</sup>*(*x*<sup>∗</sup>) *g*(*x*<sup>∗</sup>) *<sup>f</sup>*(*x*<sup>∗</sup>), and, in particular, on the slope of the prey zero-growth isocline *g*(*x*<sup>∗</sup>) *<sup>f</sup>*(*x*<sup>∗</sup>) = *<sup>x</sup>*<sup>−</sup>*<sup>α</sup>*[*x*(*<sup>α</sup>*−<sup>2</sup>)+<sup>1</sup>−*<sup>α</sup>*] *a <sup>x</sup>*=*x*<sup>∗</sup> (see also Gause [17] and Gause et al. [18]). We obtain that the equilibrium *E*∗ is stable if *mae* > *<sup>α</sup>*−<sup>1</sup> *<sup>α</sup>*−<sup>2</sup>*<sup>α</sup>* ∈ (0, 1) (check Table 1 for a summary of the feasibility and stability conditions).

We conclude that a transcritical bifurcation occurs at *mae* = 1, where the interior equilibrium exchanges stability with the predator-free equilibrium. A Hopf bifurcation appears at *mae* = *<sup>α</sup>*−<sup>1</sup> *<sup>α</sup>*−<sup>2</sup>*<sup>α</sup>*, as the eigenvalues of the community matrix become purely imaginary and the system converges to a stable limit cycle.

**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.


*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 herdlinear response *f*(*x*) = *axα*, 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, *gx*(*x*) = *rx*<sup>1</sup> − *xKx* and *gy*(*x*) = *sx*<sup>1</sup> − *xKy* , respectively, where *r* is the net growth rate of the prey and *Kx* their carrying capacity, while *s* denotes the predators' reproduction rate, i.e., not discounted by deaths,

$$\frac{d\mathbf{x}}{dt} = -r\mathbf{x}\left(1 - \frac{\mathbf{x}}{K\_{\mathbf{x}}}\right) - a\mathbf{x}^{a}y\_{\prime} \tag{18}$$

$$\frac{dy}{dt} = -sy\left(1 - \frac{y}{K\_y}\right) + \epsilon ax^a y - my. \tag{19}$$

In this way, the predators are subjected to intraspecific competition, which occurs at rate *sKy*.

To check that the populations do not grow unbounded, we set *ψ*(*t*) = *x*(*t*) + *y*(*t*) and, by repeating the steps in Section 3.1, we ge<sup>t</sup> the differential equation for the total population

$$\frac{d\psi(t)}{dt} + m\psi(t) \le \max\_{\mathbf{x}, \mathbf{y}} \left\{ rx \left( 1 - \frac{\mathbf{x}}{K\_{\mathbf{x}}} \right) + sy \left( 1 - \frac{\mathbf{y}}{K\_{\mathbf{y}}} \right) + m\mathbf{x} \right\}. \tag{20}$$

We differentiate the right-hand term with respect to *x* and *y* to ge<sup>t</sup> the local maximum (*r*+*m*) 2*r Kx*, *Ky*2 . By substitution in the equation above, we obtain

$$\frac{d\psi(t)}{dt} + m\psi(t) \le \frac{(r+m)^2}{4r} \mathcal{K}\_x + \frac{s}{4} \mathcal{K}\_y \equiv \check{\mathcal{M}}.\tag{21}$$

Therefore, the solution for the total population reads

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

where the upper bound is applicable also for *x*(*t*) and *y*(*t*).

To obtain the non-dimensional version of the system in (18) and (19), we consider the dimensionless variables and parameters *x*˜ = *xKx* , *y*˜ = *yKy* , ˜*t* = *rt*, *a*˜ = *Ky K*1−*<sup>α</sup> x r a*, *s*˜ = *sr*, *e*˜ = *Kx Ky e*, *m* ˜ = *m r* . We drop the tildes and obtain the dimensionless system

$$\frac{d\mathbf{x}}{dt} = \mathbf{x}(1-\mathbf{x}) - a\mathbf{x}^a y\_\prime \tag{23}$$

$$\frac{dy}{dt} = -sy(1-y) + eax^{a}y - my,\tag{24}$$

with *gx*(*x*) = *x*(<sup>1</sup> − *<sup>x</sup>*), *gy*(*y*) = *sy*(<sup>1</sup> − *y*) and *f*(*x*) = *axα*.

We proceed with computing the equilibria. The trivial equilibria are

$$E\_0 = (0,0), \quad E\_1 = (1,0), \quad E\_2 = \left(0, 1 - \frac{m}{s}\right),\tag{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 - a}}{a} \tag{26}$$

and

$$y = 1 - \frac{m}{s} + \frac{ea}{s} \mathbf{x}^a. \tag{27}$$

Note that the isocline in (26) intersects the *x*-axis at (0, 0) and (1, 0) and has a maximum at *x* = 1−*<sup>α</sup>* 2−*α* < 12 for 0 < *α* < 1, while the isocline in (27) intersects the vertical axis at (0, 1 − *ms* ) and is a root function translated by 1 − *ms* and dilated by *eas* . 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*∗<sup>1</sup> and *E*∗<sup>2</sup> if 1 − *ms* < 1 *<sup>a</sup>*(<sup>2</sup>−*<sup>α</sup>*)<sup>1</sup>−*<sup>α</sup>* <sup>2</sup>−*<sup>α</sup>*<sup>1</sup>−*<sup>α</sup>* − *eas* <sup>1</sup>−*<sup>α</sup>* <sup>2</sup>−*<sup>α</sup><sup>α</sup>*, or be tangent at *E*∗ = <sup>1</sup>−*<sup>α</sup>* 2−*α* , 1 − *ms* + *eas* <sup>1</sup>−*<sup>α</sup>* <sup>2</sup>−*<sup>α</sup><sup>α</sup>* when 1 − *ms* = 1 *<sup>a</sup>*(<sup>2</sup>−*<sup>α</sup>*)<sup>1</sup>−*<sup>α</sup>* <sup>2</sup>−*<sup>α</sup>*<sup>1</sup>−*<sup>α</sup>* − *eas* <sup>1</sup>−*<sup>α</sup>* <sup>2</sup>−*<sup>α</sup><sup>α</sup>* or never intersect in *x* ∈ (0, 1) when 1 − *ms* > 1 *<sup>a</sup>*(<sup>2</sup>−*<sup>α</sup>*)<sup>1</sup>−*<sup>α</sup>* <sup>2</sup>−*<sup>α</sup>*<sup>1</sup>−*<sup>α</sup>* − *eas* <sup>1</sup>−*<sup>α</sup>* <sup>2</sup>−*<sup>α</sup><sup>α</sup>*. The equilibria are obtained as the positive roots of the curve

$$\phi(\mathbf{x}) = 1 - \frac{m}{s} + \frac{\varepsilon a}{s} \mathbf{x}^{\alpha} - (1 - \mathbf{x}) \frac{\mathbf{x}^{1-\alpha}}{a} \tag{28}$$

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

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

When *m* > *s*, the isocline in (27) intersects the vertical axes at *y* = 1 − *m s* < 0 and we find at most one interior equilibrium *E*<sup>∗</sup>. We obtain the feasibility condition for *E*∗ 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

$$f(x,y) = \begin{bmatrix} 1 - 2x - ax^{a-1}ay & -ax^a \\ aex^{a-1}ay & s - 2sy + aex^a - m \end{bmatrix} . \tag{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 Tables 2 and 3.

We rewrite the Jacobian matrix evaluated at the interior equilibrium in terms of functions *f*(*x*), *gx*(*x*) and *gy*(*y*):

$$J(\mathbf{x}^\*, y^\*) = \begin{bmatrix} \mathbf{g}\_x'(\mathbf{x}^\*) - f'(\mathbf{x}^\*)y^\* & -f(\mathbf{x}^\*) \\ \varepsilon f'(\mathbf{x}^\*)y^\* & \mathbf{g}\_y'(y) + \varepsilon f(\mathbf{x}^\*) - m \end{bmatrix} = \begin{bmatrix} f(\mathbf{x}^\*) \begin{pmatrix} \frac{\mathbf{g}\_x(\mathbf{x}^\*)}{f(\mathbf{x}^\*)} \end{pmatrix}^\prime & -f(\mathbf{x}^\*) \\ \varepsilon f'(\mathbf{x}^\*)y^\* & -y^\* \end{bmatrix}.\tag{31}$$

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

$$\text{tr}f(\mathbf{x}^\*, y^\*) \quad = \quad f(\mathbf{x}^\*) \left( \frac{g\_{\mathbf{x}}(\mathbf{x}^\*)}{f(\mathbf{x}^\*)} \right)' - y^\*, \tag{32}$$

$$\det f(\mathbf{x}^\*, y^\*) \quad = \quad -y^\* f(\mathbf{x}^\*) \left( \frac{g\_{\mathbf{x}}(\mathbf{x}^\*)}{f(\mathbf{x}^\*)} \right)' + \epsilon f(\mathbf{x}^\*) f'(\mathbf{x}^\*) y^\*. \tag{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 tr*J*(*x*<sup>∗</sup>, *y*<sup>∗</sup>) < 0 and det *J*(*x*<sup>∗</sup>, *y*<sup>∗</sup>) > 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.

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


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


*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*(*x*) = *ax<sup>α</sup>* 1+*ahxα* 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*(*x*) = *rx*-1 − *xK* for the prey, where *r* denotes their net growth rate and *K* their carrying capacity,

$$\frac{d\mathbf{x}}{dt}\_{\cdot} = -r\mathbf{x}\left(1 - \frac{\mathbf{x}}{K}\right) - \frac{a\mathbf{x}^{a}}{1 + ah\mathbf{x}^{a}}y\_{\prime} \tag{34}$$

$$\frac{dy}{dt} \quad = \quad e \frac{a\mathbf{x}^a}{1 + ah\mathbf{x}^a}y - my. \tag{35}$$

The total population *ψ*(*t*) = *x*(*t*) + *y*(*t*) verifies

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

with *M* ¯ = *<sup>K</sup>*(*r*+*m*)<sup>2</sup> 4*r*as in Section 3.1.

 We obtain the dimensionless version of the model by applying the substitutions *x*˜ = *xK* , *y* ˜ = *yK* , ˜*t* = *rt*, *a*˜ = *aKα r* , ˜*h* = *rh*, *m*˜ = *mr* . We drop the tildes and ge<sup>t</sup> the equations

$$\frac{d\mathbf{x}}{dt} = \mathbf{x}(1-\mathbf{x}) - \frac{a\mathbf{x}^a}{1+ah\mathbf{x}^a}y\_{,} \tag{37}$$

$$\frac{dy}{dt} \quad = \quad \varepsilon \frac{a\mathbf{x}^a}{1 + ah\mathbf{x}^a}y - my,\tag{38}$$

with *g*(*x*) = *x*(<sup>1</sup> − *x*) and *f*(*x*) = *ax<sup>α</sup>* 1+*ahxα*

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, <sup>0</sup>), while the unique interior equilibrium *E*∗ = (*x*<sup>∗</sup>, *y*<sup>∗</sup>) has full expression

.

$$E\_\* = \left( \left( \frac{m}{ea - mah} \right)^{\frac{1}{a}}, \frac{e}{m} x^\* (1 - x^\*) \right) \tag{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) = \begin{bmatrix} 1 - 2x - \frac{a ax^{a-1}}{(1 + a b x^a)^2} y & \frac{a x^a}{1 + a b x^a} \\ \varepsilon \frac{a ax^{a-1}}{(1 + a b x^a)^2} y & \varepsilon \frac{a x^a}{1 + a b x^a} - m \end{bmatrix}.\tag{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(\mathbf{x}^\*, y^\*) = \begin{bmatrix} y'(\mathbf{x}^\*) - f'(\mathbf{x}^\*)y^\* & -f(\mathbf{x}^\*) \\ ef'(\mathbf{x}^\*)y^\* & ef(\mathbf{x}^\*) - m \end{bmatrix} = \begin{bmatrix} f(\mathbf{x}^\*) \begin{pmatrix} \frac{g'(\mathbf{x}^\*)}{f(\mathbf{x}^\*)} \end{pmatrix}' & -\frac{m}{\epsilon} \\ ef'(\mathbf{x}^\*)y^\* & 0 \end{bmatrix}. \tag{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, *g*(*x*<sup>∗</sup>) *<sup>f</sup>*(*x*<sup>∗</sup>) = *<sup>x</sup>*<sup>−</sup>*<sup>α</sup>*[*x*(*<sup>α</sup>*−<sup>2</sup>)+<sup>1</sup>−*α*+*ahx<sup>α</sup>*(<sup>1</sup>−2*x*)] *a <sup>x</sup>*=*x*<sup>∗</sup> . For the same reason as in Section 3.1, the system in (37) and (38) undergoes a Hopf bifurcation when *g*(*x*<sup>∗</sup>) *<sup>f</sup>*(*x*<sup>∗</sup>) = 0 and converges to a stable limit cycle for *g*(*x*<sup>∗</sup>) *<sup>f</sup>*(*x*<sup>∗</sup>) > 0, otherwise it converges to the interior equilibrium *E*<sup>∗</sup>. In Table 4 we give a summary of the feasibility and stability conditions of the equilibria.


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

*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*(*x*) = *ax<sup>α</sup>* 1+*ahxα* 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, *gx*(*x*) = *rx* 1 − *x Kx* and *gy*(*x*) = *sx* 1 − *x Ky* , respectively, where *r* is the net growth rate of the prey and *Kx* their carrying capacity, while *s* denotes the predators' reproduction rate,

$$\frac{dx}{dt} \quad = \quad rx\left(1 - \frac{x}{K\_x}\right) - \frac{ax^a}{1 + ahx^a}y\_{,} \tag{42}$$

$$\frac{dy}{dt} = -sy\left(1 - \frac{y}{K\_y}\right) + c\frac{a\mathbf{x}^a}{1 + ah\mathbf{x}^a}y - my. \tag{43}$$

Once again, note the second term in the predators' equation, whose coefficient *s Ky* models predators intraspecific competition.

The boundedness of the populations is ensured by the condition on the total population density *ψ*(*t*) = *x*(*t*) + *y*(*t*)

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

with *M*¯ = (*r*+*m*)<sup>2</sup> 4*r Kx* + *s* 4*Ky* as in Section 3.1.

We use the dimensionless quantities *x*˜ = *x Kx* , *y*˜ = *y Ky* , ˜*t* = *rt*, *a*˜ = *Ky K*1−*<sup>α</sup> x r a*, ˜ *h* = *r Kx Ky h*, *s*˜ = *s r*, *e*˜ = *Kx Ky e*, *m*˜ = *m r* to obtain the dimensionless system of equations

$$\frac{d\mathbf{x}}{dt} = \mathbf{x}(1-\mathbf{x}) - \frac{a\mathbf{x}^a}{1+ah\mathbf{x}^a}y\_{\prime} \tag{45}$$

$$\frac{dy}{dt} = -sy(1-y) + \varepsilon \frac{a\mathbf{x}^a}{1+ah\mathbf{x}^a}y - my,\tag{46}$$

with *gx*(*x*) = *x*(<sup>1</sup> − *<sup>x</sup>*), *gy*(*y*) = *sy*(<sup>1</sup> − *y*) and *f*(*x*) = *ax<sup>α</sup>* 1+*ahxα* .

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

$$E\_0 = (0,0), \quad E\_1 = (1,0), \quad E\_2 = \left(0, 1 - \frac{m}{s}\right),\tag{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+ahx^a)}{ax^a} \tag{48}$$

and

*y* = 1 − *ms* + *eax<sup>α</sup> s*(<sup>1</sup> + *ahxα*). (49)

The isocline in (48) intersects the horizontal axis at (0, 0) and (1, <sup>0</sup>), while the isocline in (49) intersects the vertical axis at (0, 1 − *ms* ). 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:

$$f(\mathbf{x}, y) = \begin{bmatrix} 1 - 2\mathbf{x} - \frac{a\mathbf{x}\mathbf{x}^{a-1}}{(1 + a\mathbf{h}\mathbf{x}^a)^2} y & \frac{a\mathbf{x}^a}{1 + a\mathbf{h}\mathbf{x}^a} \\\ e\frac{a\mathbf{x}\mathbf{x}^{a-1}}{(1 + a\mathbf{h}\mathbf{x}^a)^2} y & s - 2sy + e\frac{a\mathbf{x}^a}{1 + a\mathbf{h}\mathbf{x}^a} - m \end{bmatrix}. \tag{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* + *ea* 1+*ah* , an unstable saddle otherwise. We give a summary of these results in Table 5.

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


#### *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, *α*, when possible. Additionally, we obtain the two-parameter bifurcation diagrams with respect to the parameter pairs (*<sup>m</sup>*, -) or (*<sup>α</sup>*, -), 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 and we call it the nominal set (hypothetical values).

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


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.

**Figure 1.** One-parameter bifurcation diagram with respect to *m*. **Left panel**: the prey population dynamics. **Right panel**: the predator population dynamics. Thick lines: stable equilibria; dashed lines: unstable equilibria; HB: supercritical Hopf bifurcation point where a stable limit cycle arises with maximum amplitude given by the amplitude of the HB (circles); TB: transcritical bifurcation point where the coexistence equilibrium exchanges its stability with the prey-only equilibrium. The remaining parameter values are as in Table 6.

To complete the analysis, in Figure 2 we have plotted all the possible two-parameter bifurcation diagrams for (*<sup>m</sup>*, -), with - = *a*, *K*, *a*, *α* or *e*. Note that the HB curve appears in all two-parameter bifurcation diagrams, but only in the plots where we vary (*<sup>m</sup>*, *<sup>K</sup>*), (*<sup>m</sup>*, *a*) and (*<sup>m</sup>*,*<sup>e</sup>*) 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 *α* to change, the HB occurs only for *m* ≤ 0.5 (Figure 2, bottom left).

**Figure 2.** Two-parameter bifurcation diagrams with respect to *m* and, from **top-left** to **bottom-right**, *r*, *K*, *a*, *α* and *e*. Thick line: HB curve; dashed line: transcritical bifurcation curve. *E*∗ is the coexistence equilibrium and *E*1the prey-only equilibrium. The remaining parameter values are as in Table 6.

As a second example, we analyse the predator–prey dynamics with generalist predator and herd-linear functional response. In Table 7 we give the nominal set of parameter values for the model in (18) and (19).


**Table 7.** 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, *α*. Here a supercritical HB occurs at *α* = 0.5995 and a subcritical HB appears at *α* = 0.1476.

**Figure 3.** One-parameter bifurcation diagram with respect to *α*. **Left panel**: the prey population dynamics. **Right panel**: the predator population dynamics. Thick lines: stable equilibria; dashed lines: unstable equilibria; SHB: supercritical Hopf bifurcation point where a stable limit cycle arises with maximum amplitude given by the amplitude of the SHB (circles); UHB: subcritical Hopf bifurcation point, where an unstable limit cycle arises with maximum amplitude given by the amplitude of the UHB (pluses); TB: transcritical bifurcation point where the coexistence equilibrium exchanges its stability with the prey-only equilibrium. The remaining parameter values are as in Table 7.

The two-parameters bifurcation diagrams with respect to (*<sup>α</sup>*, -), with - = *r*, *Kx*, *a*, *s*, *Ky*, *e*, and *m* for the model with Equations (18) and (19) are given in Figure 4. When we vary the parameter pair (*<sup>α</sup>*, *<sup>a</sup>*), a HB appears for all values of *α* 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 (*<sup>α</sup>*, *<sup>m</sup>*), the transcritical bifurcation curve is present. Finally, if one of the parameters - = *Kx*, *a*, *s*, *Ky*, 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 we list the parameter values for the model in (34) and (35).


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

**Figure 4.** Two-parameter bifurcation diagram with respect to *α*, and, from top left to bottom right, *r*, *Kx*, *a*, *s*, *Ky*, *e* and *m*. Thick line: HB curve; dashed line: transcritical bifurcation curve. *E*∗ is the coexistence equilibrium and *E*1 the prey-only equilibrium. The remaining parameter values are as in Table 7.

In Figure 5 we obtain qualitatively similar results as for (6) and (7) for the oneparameter 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 (*<sup>m</sup>*, *K*) and (*<sup>m</sup>*, *a*) the HB curve is more concave; when we vary the parameter pair (*<sup>m</sup>*, *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.

**Figure 5.** One-parameter bifurcation diagram with respect to the natural mortality rate of the predators, *m*. **Left panel**: the prey population dynamics; **Right panel**: the predator population dynamics. Thick lines: stable equilibria; dashed lines: unstable equilibria. HB: supercritical Hopf bifurcation point where a stable limit cycle arises with maximum amplitude given by the amplitude of the HB (circles); TB: transcritical bifurcation point where the coexistence equilibrium exchanges its stability with the prey-only equilibrium. The remaining parameter values are as in Table 8.

**Figure 6.** Two-parameter bifurcation diagram with respect to *m* and, from top left to bottom right, *r*, *K*, *a*, *α*, *h* and *e*. Thick line: HB curve; dashed line: transcritical bifurcation curve. *E*∗ is the coexistence equilibrium and *E*1 the prey-only equilibrium. The remaining parameter values are as in Table 8.

Finally, we study the predator–prey dynamics with generalist predator and herd-Holling type II functional response. In Table 9 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 (*<sup>m</sup>*, -), with - = *r*, *Kx*, *a*, *α*, *h*, *s*, *Ky* and *e* show behaviours similar to the previous models, see Figures 7 and 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.

**Figure 7.** One-parameter bifurcation diagram with respect to *m*. **Left panel**: the prey population dynamics. **Right panel**: the predator population dynamics. Thick lines: stable equilibria; dashed lines: unstable equilibria; HB: supercritical Hopf bifurcation point where a stable limit cycle arises with maximum amplitude given by the amplitude of the HB (circles); TB: transcritical bifurcation point where the coexistence equilibrium exchanges its stability with the prey-only equilibrium. The remaining parameter values are as in Table 9.

**Figure 8.** Two-parameters bifurcation diagram with respect to *m* and, from top left to bottom right, *r*, *Kx*, *a*, *α*, *h*, *s*, *Ky*, and *e*. Thick line: HB curve; dashed line: transcritical bifurcation curve. *E*∗ is the coexistence equilibrium and *E*1 the prey-only equilibrium. The remaining parameter values are as in Table 9.


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

## **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–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 α* (see [5–10]), we model the density of the prey *x* on the edge of the herd as *x<sup>α</sup>* 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–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.

**Author Contributions:** Formal analysis, C.B.; investigation, C.B. and I.M.B.; methodology, C.B., I.M.B. and E.V.; software, I.M.B.; writing—review and editing, C.B., I.M.B. and E.V. All authors have read and agreed to the published version of the manuscript.

**Funding:** C. Berardo was supported by the Academy of Finland, Centre of Excellence in Analysis and Dynamics Research. I. M. Bulai was supported by MIUR through PON-AIM Linea 1 (AIM1852570-1).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** This research has been accomplished within the UMI Group TAA "Approximation Theory and Application".

**Conflicts of Interest:** The authors declare no conflict of interest.
