How do predator interference, prey herding and their possible retaliation affect prey-predator coexistence?

1.   Introduction

Different dynamics of prey-predator interaction have been considered over the years, starting with the classical Lotka-Volterra model. Of interest to us is the case in which prey with herd behavior and individualistic generalist predators coexist in a common area. For example, we can refer to a savanna, where large herbivores live alongside some predators that feed on them. A savanna is a mainly subtropical and tropical terrestrial biome that is mainly characterized by predominantly grassy vegetation with fairly widely spaced shrubs and trees; it is found in many transition zones between rainforest and desert or steppe regions in central Africa, Central and South America, India, Indochina, and Australia ^[1]. In this context, the choice to live in groups is one of the strategies that herbivores use to defend themselves against predator attacks and is therefore classified as a behavioral defense mechanism ^[2,3,4].

Predator hunting on prey herds can be modeled by observing that the individuals most likely to be affected by the predator attack are those on the herd perimeter. Therefore, denoting by $N$ the number of prey per spatial unit as a function of time, the number of possible prey caught should be proportional to $\sqrt{N}$; see ^[5]; for more general case, it should be proportional to $N^{\alpha}$, as described in ^[6], with the herd shape index being constrained as follows: $\frac{1}{2} \leq \alpha < 1$. Based on this assumption, several models that also incorporate the predator satiation phenomenon, ^{[7,8,9,10,11]}, and the presence of disease, e.g., those detailed in^[12,13,14], have been proposed. In addition, the prey response to predator attacks was considered in ^[15]. This is in agreement with the fact that, not only is aggregation advantageous for prey in terms of defense, the benefits tend to increase with growing herd size. Indeed, a large herd of herbivores can discourage hunting attempts by predators, who can become frightened, fearing that they in turn will be attacked in defense and injured, especially if the first ones have passive defense structural elements such as being large in size or having strong horns ^[16]. This prey-induced fear of predators is not to be confused with predator-induced fear of prey, in ^{[17,18,19,20,21]}.

In this paper, we propose a generalization of the two-population model presented in ^[15]. In particular, still focusing on individualistic generalist predators and prey that live in herds, we propose a new predator-prey interaction term. This new term has a denominator of the Beddington-DeAngelis form, ^[22,23], which contains a binary-value parameter that allows to take into account the prey response to predator attacks or not, depending on the species considered. In addition, this term presents a function $\Phi (N, \alpha)$ in the numerator, instead of $N^{\alpha}$, that behaves as $N$ for small values of $N$, and as $N^{\alpha}$ for large values of $N$. This function has been carefully constructed in such a way as to prevent singularity problems, similar to what was done in ^[24,25].

The article is organized as follows. The next section contains the details of the new model and its analytical study. In particular, the possible equilibrium points are derived with their feasibility and local stability conditions. In Section 3, on the other hand, the possible equilibrium bifurcations and the bistability of the model are numerically analyzed. We used MATLAB to identify the transcritical bifurcations from the predator-only equilibrium to the coexistence equilibrium and saddle-node bifurcations of the coexistence equilibrium. In addition, by using the MATLAB toolbox bSTAB, introduced in ^[26], bistability was identified and analyzed, as presented here with particular emphasis on the basin stability values; see ^[27,28] for further details. The paper ends with a discussion in which the new system is compared with some models in the literature. Of particular importance are Figures 10, 11, 12, and 13 which compare the population values at coexistence in the $N$-$P$ plane all of these models.

2.   Materials and methods

2.1.   Model equations

We propose here a simple model, particularly, a generalization of the two-population system presented in ^[15], to describe the interaction dynamics of generalist predators, which behave individualistically, and prey herds. In this case, prey may or may not be able to respond to predator attacks. This depends on the particular prey species that is being considered.

Herd behavior is generally represented by a power function, as described in ^[5,6], that can possibly be modified to account for the predator satiation phenomenon. This function can be modeled, as in ^[15], by the classical Holling type Ⅱ (HTⅡ) functional response ^[30], with the prey power function modified in the denominator to account for the prey response to predator attacks.

Here, at first, we choose to generalize this fraction by replacing the prey power function in the numerator with a function $\Phi (N, \alpha)$, with $\frac{1}{2} \leq \alpha < 1$, such that

and

Specifically, the function we select is given by

Since $\Phi(0, \alpha) = 0$ and

$\Phi(N, \alpha)$ is a nonnegative increasing function. Its behavior with $\alpha = \frac{1}{2}$ is shown in Figure 1. Using the function $\Phi(N, \alpha)$ instead of $N^{\alpha}$ presents advantages from a mathematical point of view because it avoids the singularity problems associated with $N^{\alpha}$. These problems incur in the study of local stability as a result of evaluating the Jacobian at equilibria under the condition that the first component is equal to zero. Moreover, from a modeling point of view, this function, $\Phi(N, \alpha)$, better describes the real situation than does $N^{\alpha}$. Indeed, for large values of $N$, it behaves as $N^{\alpha}$. For small values of $N$, $\Phi(N, \alpha)$ behaves like $N$, which is in agreement with the fact that as the members of a herd decrease, they each tend to interact individually with the predator population. The function proposed here is not the only function that exhibits the desired features; see ^[24,25], in which similar considerations for $\alpha = \frac{1}{2}$ are introduced.

Secondly, we chose to use a fraction denominator of the Beddington-DeAngelis form; see ^[22,23]. In addition, it contains a binary-valued parameter $q$ which, for $q = 1$, may account, or, for $q = 0$, may not account for prey response to predator attacks.

In summary, the model equations we propose are as follows:

where $N$ and $P$ respectively represent prey and predators, with the function $\Phi (N, \alpha)$ as defined in (2.1) and the parameter restrictions as shown below:

In this system, the predator's individualistic behavior is represented by the classical unit exponent in the numerator and the denominator of the fractional interaction terms. Furthermore, it is assumed that prey grow logistically with a reproduction rate $r$ and carrying capacity $K$, while predators behave similarly, with a growth rate $s$ and carrying capacity $H$. In particular, the logistic growth in (2.3) is due to the assumption that the predators are generalists. Hence, they have additional available resources that are not explicitly modeled. The remaining parameters indicate the predation rate $a$, the assimilation coefficient $e$, the herd shape index $\alpha$, the prey defense parameter $q$, and the Beddington-DeAngelis coefficients for prey and predators, $b$ and $c$, respectively. Table 1 summarizes the model parameters. We have assumed the restrictions stated in (2.4) and that all other parameters are positive.

Note that, if $q = 0$ as the number of prey increases, keeping fixed the number of predators $P$, the second term of (2.2) tends to $-ab^{-1}P$, and that of (2.3) tends to $eab^{-1}P$. Instead, if $q = 1$, both terms in (2.2) and (2.3) go to zero. This follows from the fact that, if $P$ is fixed, the function

has the following limit:

On the other hand, if the numbers of prey and predators increase simultaneously, $g(N, P)$ diverges positively if $q = 0$, and it tends to $b^{-1}$ if $q = 1$. The behavior of the function $g(N, P)$, for $q = 0$ and $q = 1$, is shown in Figure 2, for fixed parameter values $\alpha = \frac{1}{2}$ and $b = c = 1$. Thus, if $q = 1$, the model accounts for the prey response to predator attacks since, as the number of individuals in the prey herd increases, predation decreases, whereas, if $q = 0$, this is not the case.

2.2.   Model analysis

2.2.1.   Equilibrium points and their feasibility

By solving the system of equilibrium equations obtained by setting to zero the right-hand sides of (2.2) and (2.3), three unconditionally feasible equilibria can be immediately identified. These are the origin, $E_0 = (0, 0)$, and two equilibria with only one vanishing population, i.e., $E_N = (K, 0)$ and $E_P = (0, H)$.

The process of analytically determining the coexistence equilibrium of the form $E_{NP} = \left(N^*, P^*\right)$, with $N^* > 0$ and $P^* > 0$, as based on the equilibrium equations, turns out to be less straightforward and far more difficult. This amounts to solving a highly nonlinear algebraic system. Thus, we can only numerically determine the explicit coexistence equilibrium coordinates. However, we can also rely on a graphic interpretation that will provide conditions for feasibility. The general idea is to interpret the two equilibrium equations as curves in the phase plane and look for sufficient conditions that ensure their possible intersections, which represent the coexistence equilibria. We now proceed to detail the construction of these curves, which, in what follows, will be studied analytically to determine when they may have intersection points in the first quadrant.

By summing the two equilibrium equations, after multiplying the first equation by $e$, we obtain

which is a conic section. Moreover, from the first equilibrium equation, it is possible to explicitly get $P$ as a function of $N$, which is $P = \psi (N)$, with

The first curve, $\Gamma(N, P) = 0$, is given implicity, whereas the second one, $P = \psi(N)$, turns out to be an explicit function. If intersections between $\Gamma$, from (2.5), and $\psi$, from (2.6), exist inside of the first quadrant of the $N$-$P$ plane, these are coexistence equilibria of the model. We now seek to indentify these intersections or, at least, the conditions guaranteeing them, so as to ensure the feasibility of the coexistence point $E_{NP}$.

Note that the choice to consider the conic $\Gamma$ and the first nullcline, i.e, $P = \psi(N)$, that is instead of the two nullclines, is dictated by the difficulty of the problem and simplifies the analysis. Indeed, both nullclines have highly nonlinear algebraic expressions. In the first nullcline, however, it is possible to obtain $P$ as a function of $N$, which is not true for the second nullcline, nor can one find $N$ as a function of $P$. Consequently, the second nullcline turns out to be particularly problematic. Determining its qualitative behavior analytically is very complicated. Moreover, numerically obtaining its graphical representation is also difficult. The MATLAB built-in functions for graphical representation of its implicit expression also fail in this case. For these reasons, we found that it is possible to get the simpler conic $\Gamma$ from the sum of the two equilibrium equations and consider it in place of the second nullcline.

We begin by studying the conic $\Gamma$. To assess its nature, we need to determine its invariants.

Note that for each conic in non-homogeneous coordinates of the general form

two relevant matrices can be considered: the matrix of coefficients, $M_C = (a_{ij})$, which is a symmetric matrix of order 3, and the matrix of quadratic terms, which is the submatrix $S_C$ of $M_C$ that is obtained by deleting the third row and third column. See, for example, ^[31,32]. Further, $\det(M_C)$ and $\det(S_C)$ are respectively called cubic and quadratic invariants of $C$.

The two matrices associated with the conic (2.5) are given by

They have the following determinants, respectively:

Consequently, using the classification scheme of conics in non-homogeneous coordinates ^[31,32], we can say that the conic is non-degenerate and it is an ellipse. The axes of this ellipse are the two perpendicular lines $N = \frac{K}{2}$ and $P = \frac{H}{2}$. Then, its center is $C = \left(\frac{K}{2}, \frac{H}{2}\right)$ and its vertices are given by

with

Furthermore, the conic intersects the coordinate axes at the origin and two other points, i.e., $W_1 = (K, 0)$ and $W_2 = (0, H)$. Note also that the vertices of $\Gamma$ lie outside of the rectangle formed by the points $W_1$ and $W_2$, the origin, and $(K, H)$.

As for the function $P = \psi(N)$, note that it goes through $W_1$, as does $\Gamma$. Now, $\psi(N)$ is defined only for the values $N > 0$ such that

Observe also that

Assume that (2.7) holds. In view of the fact that $\psi(N)$ does not have zeros other than $N = K$, in the interval $[0, K]$, it is positive only if

Further, $\psi(N)$ itself will intersect $\Gamma(N, P)$, if at some point is lower than $\Gamma(N, P)$ and at other points is greater. Since both vanish at $N = K$, it is enough to require that the height of $\psi$ at the origin is above $W_2$, i.e., the condition $\psi(0) > H$ must be satisfied. Hence, the first set of sufficient conditions to ensure at least one coexistence equilibrium is defined by (2.8), together with

If the above condition (2.9) is not satisfied, $\psi$ is negative in $[0, K]$ and no feasible $E_{N, P}$ exists.

We now look into the situation in which (2.7) fails to hold. In such case, $\psi$ admits vertical asymptotes.

Note that, if they exist, they must lie to the left of $W_1$, i.e., for $0 < N < K$. Indeed, $\Phi(N, \alpha) > 0$ and the function $\psi$ must be positive on the right-hand side of (2.7) to have an equality. In the same interval, we also have $N^{(\psi)}(N) > 0$. Moreover, since, for $N > K$, the numerator and denominator of $\psi(N)$ have constant signs, $N^{(\psi)}(N) < 0$ and $D^{(\psi)}(N) > 0$, it follows that $\psi(N) < 0$ for $N > K$, that is, it is negative to the right of $W_1$, with

In addition, we can observe that

Let us consider separately the two functions in (2.7). The former is the increasing nonnegative function $\Phi(N, \alpha)$, while the latter, i.e.,

is a concave parabola through the origin, and $W_1$. Their graphs certainly intersect at the origin and possibly at other points within the first quadrant, which correspond to as many vertical asymptotes of $\psi(N)$. Since

a sufficient condition to have at least one vertical asymptote of $\psi(N)$ is the opposite one of (2.8), i.e.,

Let us say that the vertical asymptote is located at some $0 < N_1 < K$. Since $\psi'(K) < 0$, the function $\psi$ in $(N_1, K]$ has a positive branch descending from $+\infty$ at $N = N_1$ down to zero at $N = K$. Considering the continuity of both $\psi$ and $\Gamma$ in $(N_1, K]$ and the fact that $\Gamma(N_1, \alpha)$ is finite, there must be an intersection in the positive stripe $(N_1, K]\times R^+$, ensuring a coexistence equilibrium $E_{NP}$ for the system.

We can elaborate more on this situation. No other intersection exists in $[0, N_1]$ in the case that (2.10) holds because $\psi$ in such an interval is negative and it is negative at the origin. Indeed, we have

If more asymptotes exist to the left of $N_1$, for each of these asymptotes, $N = N_k$, we must have

because, once again, the function $\psi$ does not have zeros in $[0, K)$. Hence, no other feasible intersections can exist between $\psi$ and $\Gamma$.

By combining the conditions (2.10) and (2.9), we can conclude that a sufficient condition for the existence of at least one coexistence equilibrium is given by

2.2.2.   Local stability

We use the Jacobian of the system of equations (2.2) and (2.3) to study the local stability of its equilibria. This matrix is given by

with

and

where

Evaluating $J$ at the three equilibria with at least one nonvanishing population, we find that

The origin and $E_N$ are unconditionally unstable since both eigenvalues of $J(0, 0)$ are positive and the first and the second eigenvalues of $J(K, 0)$ are are respectively negative and positive. The second eigenvalue on the diagonal of the matrix $J(0, H)$, instead, is always negative, while the first one is negative if and only if

The point $E_P$ is asymptotically stable if and only if this condition is verified.

Finally, to discuss the local stability of the coexistence equilibrium, whose coordinates are not known, we can calculate the trace and the determinant of the Jacobian. The trace is given by

and the determinant is given by

with

Then, using the Routh-Hurwitz criterion, we can say that the equilibria of the form $E_{NP}$, when they are feasible, are asymptotically stable if and only if

The fact that these two stability conditions can be simultaneously satisfied for a feasible coexistence equilibrium has been verified through numerical simulations. For instance, by setting the parameter values

we can obtain an asymptotically stable coexistence equilibrium point by applying the initial conditions $N(0) = 60$ and $P(0) = 15$; see the left panel of Figure 3.

Table 2 summarizes the identified information concerning the model equilibrium points with their feasibility and local stability conditions.

3.   Results

From Table 2, the local stability condition (2.12) for the predator-only equilibrium and the feasibility condition (2.11) for coexistence are complementary. However, by performing numerical simulations, as the parameter $a$ changes, transcritical bifurcations can be identified at $a = a^{**}$, with $a^{**} \geq a^* = cr +\frac{r}{H}$, whose value changes depending on the initial conditions. Sotomayor's theorem constitutes a useful tool to analytically prove the presence of transcritical, saddle-node, and pitchfork bifurcations, provided that the coordinates of one equilibrium point are explicitly known ^[33]. In addition, it is necessary to know the exact value of the threshold for the bifurcation parameter. Thus, in this case, since we cannot exactly identify the critical value $a^{**}$, this theorem cannot be used. Consequently, we can only perform numerical simulations to visualize such transcritical bifurcations and estimate the value of $a^{**}$. For example, in Figure 4, using the parameter values (2.14) and the initial conditions $N(0) = P(0) = 30$, we observe a transcritical bifurcation from $E_{NP}$ to $E_P$ at $a = a^{**}$, with $a^{**} \approx 1.08834 > a^* = 0.20800$.

This is in agreement with the fact that the condition (2.11) is only sufficient for the feasibility of the coexistence point. Hence, $E_{NP}$ could be feasible even if (2.11) is not satisfied. Thus, there may exist regions in the parameter space for which the equilibria $E_P$ and $E_{NP}$ turn out to be simultaneously feasible and asymptotically stable. In these regions we have bistability. For instance, as illustrated in Figure 3, fixing the parameter values as detailed in (2.14), we can observe how the system can reach the two different equilibria under different initial conditions.

This bistability has been explored by using the MATLAB toolbox bSTAB, as introduced in ^[26]. In the square $Q = \left[0, 50 \right] \times \left[0, 50 \right]$, $n = 10^4$ uniformly distributed initial conditions were considered together with the relative tolerance of the numerical integration scheme $RelTol = 10^{-6}$. The basins of stability of the two equilibria mentioned above are illustrated in the right panel of the top row of Figure 5, while the left panel shows the basin stability values. Refer to ^[26,27,28] for their definitions. In particular, the basin stability value of an attractor $A$ in a $D$-dimensional state space is defined as follows:

where $\kappa \left(\mathbf {y} \right)$ is an indicator function

which determines whether $y$ belongs to the attraction basin $\mathcal {B}\left(A\right)$ and $\rho$ is a density distribution of states such that the system may be perturbed with $\int _{\mathbb {R}} \rho \left(\mathbf {y}\right) \mathrm {d} \mathbf {y} = 1$. The MATLAB toolbox bSTAB was used to estimate the basin stability values by performing Monte Carlo sampling from the density distribution $\rho$. Specifically, the two basin stability values found were $S_B (E_P) \simeq 0.13$ and $S_B (E_{NP}) \simeq 0.87$. By increasing the side of the square $Q$, $S_B (E_P)$ was found to decrease and, consequently, $S_B (E_{NP})$ increased, and vice versa.

In the bottom row of Figure 5, we can see how the values $S_B (E_P)$ and $S_B (E_{NP})$ vary as the parameter $a$ moves from 0.8 to 1.3. At first, we have that $S_B (E_P) \neq 0$ and $S_B(E_{NP}) \neq 0$ such that $S_B (E_P) + S_B(E_{NP}) = 1$; in this case, there is bistability of $E_P$ and $E_{NP}$. Then, we have that $S_B (E_P) = 1$ and $S_B (E_{NP}) = 0$; in this second case, only $E_P$ is feasible and asymptotically stable at the same time.

Figures 6 and 7 show how $S_B (E_P)$ and $S_B (E_{NP})$ change as the $c$ value varies. This parameter is of great interest since it controls the interference of predators. Recall that it appears in the denominator of the predator-prey interaction terms of the equations (2.2) and (2.3). To obtain the results shown in Figure 6, we did not consider the response of prey to predator attacks; hence, we set $q = 0$. Alternatively, Figure 7 shows the results of considering the prey response, which means that we set $q = 1$. Further, to see if there is a significant variation as the shape of the prey herd changes within the range of the shape index, we considered a low value, $\alpha = 0.5$, and a higher value, $\alpha = 0.9$, because technical constraints of bSTAB prevent to generate a figure for a continuous range of values of $\alpha$. These variations appear to be most significant in the case in which there is no prey response to predator attacks. Indeed, in Figure 6, we see that, for the lowest value of $\alpha$, the transition from having $E_P$ as the only attractor to having only $E_{NP}$ occurs earlier than in the case in which $\alpha = 0.9$. As can be seen in the right frame of the top row of Figure 6, for $\alpha = 0.5$, bistability occurs in the range of $0.95 \lessapprox c < 1.4$. As can be seen in the right frame of the bottom row of Figure 6, for $\alpha = 0.9$, either no bistability occurs, suggesting perhaps a transcritical bifurcation, or, there is a very tiny bistability region that is difficult to locate numerically given the very long processing time of bSTAB. Alternatively as shown in Figure 7, although different values of $S_B (E_P)$ and $S_B (E_{NP})$ can be observed in the intermediate phase, as $\alpha$ varies, the transitions from only stability of $E_P$ to bistability, as well as from bistability to only stability of $E_{NP}$, occur for very close values of $c$.

Finally, we can easily identify other bifurcations by examining the expressions (2.5) and (2.6). Indeed as $e$, $s$, and $H$ vary, the graphical representation of $P = \psi(N)$ remains fixed in the $N$-$P$ plane while $\Gamma(N, P) = 0$ moves; this is because only the expression of the ellipse contains these three parameters. Conversely, varying $\alpha$, $a$, $b$, $c$, and $q$ was found to keep the ellipse fixed and shift the graph of the function $\psi(N)$, whose expression is the only one that contains these five parameters. In both cases, the intersections between the two curves in the first quadrant of the $N$-$P$ plane can change, giving rise to bifurcations. For example, as shown in Figures 8 and 9, under the conditions of the fixed parameter values

we were able to observe saddle-node bifurcations of the coexistence equilibrium by varying $e$, $s$, and $H$ in the first case and $\alpha$, $a$, $b$, and $c$ in the second one. Tipically Sotomayor's theorem is also employed to analytically prove the presence of saddle-node bifurcations. However, in this case, since we cannot exactly identify the threshold values of the bifurcation parameters, the theorem is of no use. Thus, we investigated the occurrence of such bifurcations only through numerical simulations.

4.   Discussion

This new model exhibits just two possible steady states, i.e., the coexistence and the prey-free equilibrium. They can be linked by a transcritical bifurcation, as shown in Figure 4.

Figures 8 and 9 show that, in the interior of the first quadrant, there are two, one, or no intersections among $\Gamma(N, P) = 0$, from (2.5), and $P = \psi(N)$, from (2.6). For no intersection, the equilibrium attained by the system is $E_P$, i.e., the intersection of the conic section with the vertical axis. When there are two, these are coexistence equilibria, where the right one is stable and the left one is unstable, because it must separate the basins of attractions of $E_P$ and the stable branch of $E_{NP}$. Indeed, the bistability situation (the right intersection of the conic $\Gamma$ with $P = \psi(N)$ and the intersection of the conic section with the vertical axis) translates into the fact that $E_P$ remains stable and, topologically, that the left of the two intersections of (2.5) with (2.6) must be a saddle, while the right one must be an attractor. When the two curves are tangent to each other, a saddle-node bifurcation occurs.

Note that, for an increasing hunting rate, the predators wipe out the prey and continue to thrive because they are not specialists on $N$. Furthermore, Figure 3 shows that different initial conditions lead to different system outcomes with the same choice of model parameters, i.e., a bistability situation. In Figure 5, we show the basin stability values. In particular, the bottom frame shows that, if $a < \widehat{a}$, for a certain $1 < \widehat{a} < 1.1$, the size of the basin of attraction of the coexistence equilibrium is much larger than that of the prey-free equilibrium. For $a > \widehat{a}$, no prey survive and the basin of attraction of the prey-free equilibrium is the whole domain. This means that, in this case, for $a > \widehat{a}$, $E_P$ is a $100\%$ attractor. For the results shown in Figures 6 and 7, on the other hand, we consider variations of the basin stability values for $E_P$ and $E_{NP}$ as the parameters controlling predator interference, prey herding shape, and prey retaliation change.

4.1.   Comparison with simpler models

A direct analytical comparison with the two-population model in ^[15], of which the model proposed here is a generalization, is possible. The two systems have the same equilibrium points, i.e., the coexistence equilibrium with different components. However, for the model in ^[15], the unconditionally feasible equilibrium $E_P$ turns out to be always unstable. In that case, since coexistence is the only conditionally feasible and locally asymptotically stable equilibrium, there can be no transcritical bifurcations, but saddle-node bifurcations of the coexistence equilibrium have been numerically identified. A similar comparison, however, is not possible with many other previously mentioned prey-predator models, that have appeared in the literature because of different assumptions. In particular, here, we are considering generalist predators, while many of these models focus on specialist predators.

4.2.   Reformulation of some simpler models for comparison

In order to make a direct comparison by focusing only on predation terms, we chose to rewrite the models chosen for comparison, by inserting logistic growth terms into both differential equations we assume generalist predators in all cases. The general reference form is given by

where the function $f(N, P)$ is characteristic of each model. In particular, we were interested in comparing coexistence equilibria. As in ^[15] we choose to illustrate the differences graphically. We work for simplicity with $\alpha = 0.5$. There would be no substantial differences with a generic herd shape index $0.5 \leq \alpha < 1$. The four functions we consider concern predators that hunt prey that gather in herds, with or without predator satiation, predator interference and prey retaliation to predator attacks. These are respectively defined as follows:

The first one of these, as introduced in ^[5], is the simplest possible function in the case of individualistic predators and grouped prey. The second one, which appears, for example, in ^{[7,8,9,10,11,13,14]}, combines the HTⅡ functional response with prey herd behavior. The third one, proposed in ^[15], is a variant of the second function that takes into account the prey response to predator attacks. We observe that $f_1(N, P)$ is a special case of both $f_2(N, P)$ and $f_3(N, P)$ with $b = 0$. Finally, the fourth function, which is the one proposed in this paper, applies a denominator of Beddington-DeAngelis type ^[22,23], and has the function $\Phi(N, \alpha)$ in the numerator, is a generalization of the previous three. In particular, we are interested in both the case $q = 0$, in which there is no response of prey to predator attacks, and the case $q = 1$, in which this response is present.

Note that the comparison with models that apply

instead of $f(N, P)$ in the equations (4.1) and (4.2), i.e., models considering prey and predators both individualistic, particularly, the former with the classical quadratic or Holling type Ⅰ term and the latter with the Holling type Ⅱ functional response, has been already investigated. We refer the reader to the results in ^[15].

4.3.   Numerical results of comparison with the reformulation of simpler models

The graphical results obtained are shown in Figures 10, 11, 12, and 13. The parameters used were as follows:

To obtained the results shown in each subfigure, one of the parameters in (4.3) were varied within an interval, while all of the other ones remained fixed to see how the coordinates of the achieved coexistence equilibria change in the $N$-$P$ plane. The chosen initial conditions were applied as follows:

The model proposed here, with or without prey response to predator attacks, shows that, at each coexistence equilibrium, prey thrive in greater numbers and predators in smaller numbers than in the other three models considered. In particular, it can be observed that, by fixing all parameter values, the number of prey increases and the number of predators decreases from the first to the fifth case, i.e., the coexistence equilibrium moves rightward and downward. Note that, although there are four functional forms $f_i(N, P)$, $i = 1, \ldots, 4$, we split the analysis of $f_4$ into two, obtaining five cases. The fourth case corresponds to taking the function $f(N, P)$ in (4.1) and (4.2) equal to $f_4(N, P)$ with $q = 0$, i.e., without prey response to predator attacks. The fifth case corresponds to considering $f_4(N, P)$ with $q = 1$, indicating prey retaliation. This variation is more significant over the course of the first case to the third one and less so from the third case to the last one. The changes obtained for each model as the parameter values were varied also became significantly less pronounced over the course of the first to the fifth case, except for the variations of the parameters $H$ and $K$. In addition, we observed that, as $b$ and $c$ varied, the coordinates of the achieved coexistence equilibrium did not change in all systems because not all functions $f_i(N, P)$, for $i = 1, 2, 3, 4$, contained them. In particular, $f_1(N, P)$ does not contain $b$ and only $f_4(N, P)$ contains $c$. Consequently, in Figure 13, we did not observe a shifting of the coexistence equilibrium, as $b$ changed, in the left panel, and a variation in the position of this equilibrium was observed only for the last two cases, as $c$ value changed, in the right one.

5.   Conclusions

The new model has the following characteristics: the equilibrium points that it admits are the origin, the predator-only equilibrium, the predator-free equilibrium, and two, one, or no coexistence points. The first three are unconditionally feasible, while coexistence is conditionally feasible. Moreover, only the equilibrium $E_P$ and the coexistence points can be locally asymptotically stable. Consequently, these equilibria are the only two steady states to which the system can converge, regardless of the parameter values and the initial conditions.

Transcritical bifurcations from the predator-only equilibrium to the coexistence point are possible. In addition, saddle-node bifurcations of coexistence have also been numerically identified. Furthermore, there are combinations of parameter values for which, simultaneously, $E_P$ and $E_{NP}$ are admissible and locally asymptotically stable. In these cases, by fixing certain parameter values, the equilibrium toward which the system converges depends only on the initial conditions.

Situations of bistability, or, more generally, of multistability, have already been identified in the literature for different predator-prey dynamics; some examples can be found in ^{[19,34,35,36,37,38]}.

A direct analytical comparison with the two-population model in ^[15], of which this model is a generalization, has been done. Finally, the new model has been also compared graphically, in terms of coexistence equilibria, with several other prey-predator models that have appeared in the literature, considering generalist predators. The main evidence is that the new model always shows that prey and predators thrive in respectively larger and smaller numbers than in the other models considered for comparison. In particular, we have observed that by fixing all parameter values and refining the model (i.e., changing from the model that only takes into account grouped prey to the model that also takes into account the predator satiation, predator interference, and the possible prey response to predator attacks) the coexistence equilibrium moves rightward and downward. This is in agreement with the fact that the phenomenon of predator satiation and prey retaliation are advantageous for prey population growth and disadvantageous for predators. On the other hand, predator interference by itself should promote predator growth. However, in the new model, it is offset by the other two factors. Thus, the model seems to suggest that predator interference in the case of prey herds with or without response to predator attacks, does not significantly help the predator population to increase in size because the grouping of prey turns out to have a stronger influence.

Moreover, by combining the results of this study with those obtained in ^[15], we can generally remark prey that gather in herds suffer fewer losses from attacks of the predators than individualistic ones. This is also in agreement with the results of previous investigations, e.g., ^[38]. Similar results hold for herds of competing population. On the other hand, the advantages for symbiotic populations are reduced because only the outermost individuals can exploit the other populations for additional resources. A secondary but important result highlighted in ^[38] is the occurrence of tristability in the case of herd competition.

Overall, models with herds in which the response function takes different forms, among which the one presented in this paper is a further example, appear to exhibit novel, interesting features that are not observed with the classical quadratic population models and therefore should be given an appropriate amount of attention by the theoretical biologists, despite their simple conceptual formulations.

Use of AI tools declaration

The authors declare that they have not used artificial intelligence tools in the creation of this article.

Acknowledgments

The research work has been partially supported by the project "Approximation methods and models for the life sciences" of the Department of Mathematics of the University of Turin. The authors are members of the INdAM research group GNCS and the Research ITalian network on Approximation (RITA).

Author contributions

Francesca Acotto: Conceptualization, Data curation, Formal Analysis, Funding acquisition, Investigation, Methodology, Resources, Software, Validation, Visualization, Writing – original draft, Writing – review & editing; Ezio Venturino: Conceptualization, Formal Analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Supervision, Validation, Writing – original draft, Writing – review & editing. All authors have read and approved the final version of the manuscript for publication.

Conflict of interest

All authors declare no conflicts of interest in this paper.