**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, *α* in this article, also appeared in predator–prey dynamics with prey grouped in herd. In particular, the first studies comprised only the case *α* = 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–14]. In particular, analogously to Metz and Diekmann [15], we assume the predators in two states, searching and handling, and model

**Citation:** Berardo, C.; Bulai, I.M.; Venturino, E. Interactions Obtained from Basic Mechanistic Principles: Prey Herds and Predators. *Mathematics* **2021**, *9*, 2555. https://doi.org/10.3390/math9202555

Academic Editors: Maria Lucia Sampoliand MonicaBianchini

Received: 31 August 2021 Accepted: 1 October 2021 Published: 12 October 2021

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

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

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 *yS* and handling individuals *<sup>y</sup><sup>H</sup>*. 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


Therefore, the number of individuals exposed to predation is proportional to *<sup>x</sup>α*, where the exponent *α* is defined as *herding index* [1] and depends only on the characteristic shape of the herd. Examples are *α* = 12 for a circle in two dimensions or *α* = 23 for a sphere in three dimensions. Intermediate values of *α* 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 *axα*, 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 1*ax<sup>α</sup>* . When we exclude handling, the predator functional response is linear and takes the form *f*(*x*) = *ax<sup>α</sup>* (*herd-linear functional response*). Alternatively, let *h* denote the handling time, then the time between two successive catches by the same predator is 1 *ax<sup>α</sup>*+*h*.

 All in all, consider the following fast processes

$$\begin{aligned} \widehat{\bigotimes^{S}} & \stackrel{ax^{a}}{\longrightarrow} \bigotimes^{H} \\ & \stackrel{h^{n-1}}{\bigotimes^{H}} \stackrel{h^{-1}}{\longrightarrow} \bigotimes^{S} \end{aligned} \quad \text{the product meets the prey and the prey is caught,} \tag{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 *α* in the encounter rate *axα*, 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 *yS* and *<sup>y</sup><sup>H</sup>*, the equations for the fast-time dynamics reduce to

$$\frac{dy^S}{dt} \quad = \quad -a\mathbf{x}^a y^S + h^{-1}y^H,\tag{2}$$

$$\frac{dy^H}{dt} \quad = \quad +ax^a y^S - h^{-1}y^H. \tag{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 *yS* and *y<sup>H</sup>* verify *dydt* = *dy<sup>S</sup> dt* + *dyHdt* = 0), therefore, we can reduce the system to one equation by using the conservation law *y* = *yS* + *<sup>y</sup><sup>H</sup>*. By using this law and solving the equilibrium equation for *<sup>y</sup>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 + ahx^{a'}} \quad y^H = y - y^S. \tag{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*(*x*) = *ax<sup>α</sup>y<sup>S</sup> y* , and, plugging the expression for *yS* at the quasi-equilibrium in (4), we obtain

$$f(\mathbf{x}) = \frac{a\mathbf{x}^a}{1 + a\mathbf{l}\mathbf{x}^a}.\tag{5}$$

We name the functional response in (5) the *herd-Holling type II functional response*. Note that when *α* = 1 we recover the Holling type II functional response and when *α* = 2 we obtain the Holling type III functional response. However, in the present context with prey herd geometry, we necessarily have *α* ∈ (0, <sup>1</sup>). 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.
