This work is licensed under the Creative Commons Attribution 4.0 International License.
Introduction
In the study of ecological interactions in the framework of population dynamics, interactions of the Lotka-Volterra type [13] are a useful simplification. They have been used to analyze periodical fluctuations in populations [19], competition between species [12] and even to simulate evolving ecosystems [5]. Holling [9, 10] studied the data and hypotheses regulating the interactions between species, focusing on the predator-prey case, and proposed several types of “response functions” (types I, II and III) [6].
Traditional type I responses, such as those contained in the classical Lotka-Volterra system, are derived from the Mass Action Law, stating that the number of encounters between two populations is proportional to the product of both population sizes. In predators-prey interactions, the predation rate is constant and implies an unlimited growth. Indeed, if the prey population increases, the predators keep also on capturing prey without bounds. If the modeling either the handling time of the prey or the predator satiation effects is important in a given predator-prey interaction, a type II response may be used, because it includes a saturating effect. Indeed, for a fixed number of predators it bounds the predation rate. Finally, a type III response may represent a case where small prey populations that are hunted may fall to such low numbers that prevent reproduction, because animals cannot find themselves to mate, and therefore are bound to disappear. Another alternative interpretation concerning predator birds that rely on a “mental image” of the most common prey and therefore ignore similar species that are less frequent, appears in [7].
Effects dependent on the group size also may influence the predation rate, [4]. The Trafalgar Effect [18], where information is shared between members of the same group, may increase the range of predator detection of an individual that belongs to a group, when compared to solitary ones. Also the confusion effect may reduce a predator’s effectiveness when attacking prey [4]. Finally, simply by flocking together, the number of individuals that are actually exposed to an attack may reduced, in the sense that only those on the boundary of the herd may be more susceptible to an attack.
For this last form of gathering, in [1] a model is proposed that shows novel behaviour when compared to Lotka-Volterra type interactions. The proposed model is based on the observation that for non-fractal shapes lying in a two-dimensional space, the number of individuals on the boundary of the surface occupied by the herd may be approximated by a constant, related to the shape of the herd, times the square root of the total number of individuals occupying that surface. Hence the model functional response is proportional to the product of the square root of prey population and the total number of predators.
The model discussed in this work is an extension of a model proposed to describe herd behaviour, where predators cannot access the whole population [1]. In this sense, although the equations seem to represent the same saturating effects of predation as in Holling type II and III response functions, it has a quite diverse biological interpretation. The use of the square-root is associated with two-dimensional shapes, but a generalization of such term has already been proposed in [3], where three-dimensional shapes (of fish schools, for instance) or even fractal shapes can be modelled, yielding similar results.
Such approach carries a few problems. First, the proportional to the square root approximation of the number of individuals on the border gets less and less precise as the total number of individuals decreases. Secondly, small groups may not display herding behaviour, because group defense needs a certain threshold to become effective or for the animals to flock together. Thirdly, due to the shape of the square root function, such response function indicates that at very low populations, predation would be higher than in regular Lotka-Volterra type interactions. Thus grouping behaviour, below a certain threshold would increase predation. Such difference magnifies as the prey population tends to zero. Finally, using the square root leads to a certain technical issue related to a singularity in the Jacobian of the system, which, although not too problematic, leads to some difficulties in the interpretation of certain trajectories.
The goal of this paper is to propose and analyze a model that deals with those mathematical difficulties. In particular, we define a response function that behaves like a square root when prey abound and works approximately as a Holling type I response, i.e. Mass Action Law, for small prey populations. In this way, we correctly model the fact that group defense becomes less and less effective as group size decreases, tending asymptotically to an individualistic behaviour. A first attempt to deal with such mathematical difficulties has been done in [2], where piecewise continuous models were used to describe the predation when the prey population below or above a certain threshold. The results were mathematically interesting, leading to both supercritical and subcritical Hopf bifurcations in the system, with a rich variety of bifurcations in the system. Here both the model and analysis are much simpler and yet the model is able to grasp all the difficulties presented in the first approach by Ajiraldi and Venturino [1].
The paper is organized as follows. In Section 2 we present the proposed mathematical model and in Section 3 classical techniques are used for its analysis, such as linear stability analysis, the Grob-Hartman Theorem and the Poincaré-Bendixson Theorem (which is implicit in conditions for Hopf Bifurcations). In Section 4 for certain parameters combinations the system is shown to undergo a supercritical Hopf bifurcation, giving rise to sustained oscillations, and in Section 5 numerical simulations are reported to illustrate this behaviour. Section 6 presents some biological interpretations of the results obtained in the analysis of the model, we do not go deep in it, since the focus of our work is on the extension of model already known. Finally, in Section 7 comment the results obtained in the work as a whole.
The model
We propose a mathematical predator-prey model considering “large” populations, displaying herd behaviour. If the herd is too small it may not be possible to form an appropriate group defense or the boundary of the herd may be composed of the totality of the population. For such small groups it would be more reasonable to adopt a traditional mass action interaction term. Let F and R respectively denote the predators and prey populations. The interaction term can be described by
{a_2}\sqrt R F
, valid for large prey populations and by a1RF, in the case of a small number of prey.
The model is thus described via a piecewise function:
g(R) = \left\{ {\eqalign { {{a_1}R} \hfill & {{\rm{if}}\,0 \le R \le {h_1}} \hfill \cr {{a_2}\sqrt R } \hfill & {{\rm{if}}\,{h_{\rm{1}}} < R} } } \right.
where h1 represents a critical threshold of group size for effective defense. Such representation has at least two main defects. First, it represents a discountinuous transition in the benefits provided by each extra individual around the threshold value h1 (g′(R) is discontinuous in h1). This is probably not very realistic since a smooth transition between ineffective and effective group defense is expected, at least for some species. Second, the discontinuity in g′(R) affects the differentiability of the vector field that defines the dynamical system. This may interfere with the application of traditional theorems of dynamical systems such as the Existence and Uniqueness Theorem [17], which requires the continuity of the partial derivatives of the vector field.
Given the above considerations we propose the following response function:
f(R) = a{R \over {\sqrt {R + {h_2}} }}.
where a and h2 are parameters for which the biological interpretation will be explained after some calculations.
Since we are trying to create a smooth transition function f that approximates g in the extreme cases, that is, when R → 0 and when R → ∞, we impose that the relative error between g(R) and f(R), i.e. g(R)/f(R) − 1, goes to zero both when R → 0 and R → ∞. Through simple calculations, we obtain that a = a2 and h2 = (a2/a1)2. So f(R) can be understood as an approximation to g(R) with a smooth transition between the two extreme behaviours: one for the group defense
{a_2}\sqrt R
and the other for individual behaviour (a1R).
From the continuity of g(R) it follows also that h1 = (a2/a1)2, so that the value of for the threshold h2 is consistent, and we can simply write
{h_1} = {h_2} = \tilde h
, which we adopt for the remainder of this section. In Figure 1 a graph with the response functions for a2 = 1,
\tilde h = 20
and
{a_1} = {1 \over {\sqrt {20} }}
is shown.
Fig. 1
Response functions g(R) and f(R) of population R for the parameters a2 = 1,
\tilde h = 20
and
{a_1} = {1 \over {\sqrt {20} }}
.
With those observations, the interpretation of parameters a and
{h_2} = \tilde h
are quite simple, a = a2 is the coefficient of predation that takes into account the group defense/herd effects, which is affected by e.g. the shape of the herd, the agility and efficiency of the predators etc. The parameter
\tilde h
is a threshold for the transition between herd grouping and solitary behaviour, that it can be interpreted as an approximate number of individuals necessary for the formation of effective group defense.
The proposed model, with the smooth transition function (1) is given by:
\eqalign{ & {{{dR} \over {dt}} = r\left( {1 - {R \over K}} \right)R - f(R)F} \cr & {{{dF} \over {dt}} = - mF + ef(R)F.} }
The prey dynamics for R contains a logistic growth and a predation term, so that in the absence of the predators, the prey would grow at rate r to the environment carrying capacity K.
The predators are assumed to be specialist on the prey, so that in the absence of the latter, they die with mortality rate m. Their hunting is expressed by the second term, with conversion coefficient e < 1.
Boundedness of the trajectories
The ecological model is well-posed if its dependent variables cannot grow unboundedly. It is thus necessary to show that the system’s trajectories remain confined within a compact set. In the following, assume e ≤ 1, which is the case if both populations are described in biomass units.
Proposition 1.
Consider the system given by 2 and its trajectories R(t), F(t), if R(0) ≥ 0, F(0) ≥ 0 and e ≤ 1, then the solutions R(t) and F(t) are always non-negative and bounded.
Proof
Let φ(t) = R(t) + F(t), then φ(t) is a differentiable function, and taking m > 0, summing the equations in model (2) and observing that e ≤ 1 we find:
\eqalign{ & {{{d\varphi (t)} \over {dt}} + m\varphi (t) = rR\left( {1 - {R \over K}} \right) - mF + (e - 1){{aRF} \over {\sqrt {R + {\tilde h}} }} + mR + mF} \cr &\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, { \le rR\left( {1 - {R \over K} + {m \over r}} \right) = p(R).} }
The function p(R) is concave parabola, with maximum located at R⋆, and corresponding maximum value
M = p({R^*}) = {{rK} \over 4}{\left( {1 + {m \over r}} \right)^2}.
Thus,
{{d\varphi (t)} \over {dt}} + m\varphi (t) \le M.
Integrating the differential inequality, we find
\varphi (t) \le \left( {\varphi (0) - {M \over m}} \right){e^{ - mt}} + {M \over m} \le {\rm{max}}\left\{ {\varphi {\rm{(0),}}{M \over m}} \right\}
Thus, for model (2) the solution remains bounded. That trajectories remain non-negative follows directly from the facts that the coordinate axes are solutions of the homogeneous system and that the initial conditions are nonnegative, to make biological sense. The uniqueness of the solution trajectories implies that the axes cannot be crossed and therefore the claim.
Non-dimensional model
Model (2) can be written in a non-dimensional version via the rescaling
{R^*}(t) = {1 \over K}R(t)
,
{F^*}(t) = {1 \over {eK}}F(t)
, t*= mt and
\lambda = {r \over m}
,
\theta = {{ae\sqrt K } \over m}
:
\eqalign { & {{{d{R^*}} \over {d{t^*}}} = \lambda {R^*}(1 - {R^*}) - {{\theta {R^*}{F^*}} \over {\sqrt {{R^*} + (\tilde h/K)} }}} \cr & {{{d{F^*}} \over {d{t^*}}} = - {F^*} + {{\theta {R^*}{F^*}} \over {\sqrt {{R^*} + \left( {\tilde h/K} \right)} }}.} }
For simplicity, in (3) we drop the stars from the above quantities and write simply R, F and t. In this way, the system (3) becomes
\eqalign {& {{{dR} \over {dt}} = \lambda R(1 - R) - {{\theta RF} \over {\sqrt {R + h} }}} \cr & {{{dF} \over {dt}} = - F + {{\theta RF} \over {\sqrt {R + h} }}} }
where
h = \tilde h/K
represents the critical threshold of group size for effective defense in non-dimensional scaling Since we assume an observable group defense effect, it is reasonable to take
0 < \tilde h < K
, and hence 0 < h < 1.
Equilibria
To obtain the equilibria, we consider the equations
\dot R = 0
and
\dot F = 0
, respectively representing the vertical and horizontal isoclines, namely:
R = 0,\,\,F = {{\lambda (1 - R)\sqrt {R + h} } \over \theta }
and
F = 0,\,\,f(R) = {1 \over \theta } = \mu ,
which leads to
R = \mu \sqrt {R + h}
.
The intersection of the isoclines define the equilibrium points of the system. There are three such points, representing respectively ecosystem collapse E1 = (0,0), the prey-only equilibrium E2 = (1,0) and coexistence E3 = (R3,F3). Here the populations values are explicitly known:
{R_3} = {\mu \over 2}( {\mu + \sqrt {{\mu ^2} + 4h} } ),\,\,\,\,\,\,\,\,\,{F_3} = \lambda {R_3}(1 - {R_3}).
Stability and feasibility
In this section analyze the feasibility and stability of the equilibria. Since we are working with population models, we define that an equilibrium is feasible if and only if both of its coordinates are real and non-negative.
To analyze the stability, note that the Jacobian of (4) is given by
J = \left( {\matrix{ { - \lambda R + \lambda (1 - R) - {F \over {\mu \sqrt {R + h} }} + {{RF} \over {2\mu (R + h)\sqrt {R + h} }} } & - {{R \over {\mu \sqrt {R + h} }}} \cr {{F \over {\mu \sqrt {R + h} }} - {{RF} \over {2\mu (R + h)\sqrt {R + h} }}} & { - 1 + {R \over {\mu \sqrt {R + h} }}} } } \right).
The stability analysis is contained in the following propositions.
Proposition 2
Consider the system(4), with 0 < h < 1. The following statements hold.
The equilibrium point E1 = (0,0) is a saddle.
The equilibrium point E2 = (1,0) is :
a stable node if\mu > {1 \over {\sqrt {h + 1} }}
a saddle if\mu < {1 \over {\sqrt {h + 1} }}
Proof
The basic idea of the proof is to analyze the eigenvalues of Jacobian matrix at equilibria E1 and E2. Then, in the vicinity of those points, we can apply Grobman-Hartman Theorem [16] obtaining the topology of equilibria E1 and E2.
The Jacobian matrix (7) evaluated at E1 has the eigenvalues ξ1 = λ and ξ2 = −1. Thus, as ξ1 and ξ2 are real and have opposing signs, by the Grobman-Hartman Theorem, E1 is a saddle.
Similarly, the Jacobian matrix (7) evaluated at E2 has the eigenvalues
{\xi _1} = - 1 + {1 \over {\mu \sqrt {h + 1} }},\,\,\,\,\,{\xi _2} = - \lambda .
Thus, if (8) holds both ξ1 and ξ2 are real and negative. Therefore, by the Grobman-Hartman Theorem, E2 is a stable node. Now, if (9) holds ξ1 and ξ2 are real and have opposing signs, therefore E2 is a saddle.
Note that in the particular case where ξ1 = 0, that is,
\mu = {\mu ^\dagger }: = {1 \over {\sqrt {h + 1} }}
we cannot use the the Grobman-Hartman Theorem, since the hypotheses are not fulfilled. This degenerate case will be contemplated in the Proposition 4(a), which presents a transcritical bifurcation when μ passes though this specific value.
Proposition 3
Consider the system(4), with 0 < h < 1 and λ, μ > 0. Taking{\mu ^*} = {{1 - 2h} \over {\sqrt {3(1 + h)} }}and{\mu ^\dag } = {1 \over {\sqrt {h + 1} }}, the following statements hold.
The Jacobian matrix at E3 can be simplified when we consider the isoclines that define E3\lambda (1 - {R_3}) - {{{F_3}} \over {\mu \sqrt {{R_3} + h} }} = 0,\,\,\,\,\,\,\mu = {{{R_3}} \over {\sqrt {{R_3} + h} }}.
to obtain
\eqalign{ & {{{\left( {{J_{{E_3}}}} \right)}_{11}} = {{\lambda {R_3}} \over {2({R_3} + h)}}( - 3{R_3} + 1 - 2h),\,\,\,\,{{\left( {{J_{{E_3}}}} \right)}_{12}} = - 1,} \cr & {{{\left( {{J_{{E_3}}}} \right)}_{21}} = {{\lambda (1 - {R_3})({R_3} + 2h)} \over {2({R_3} + h)}},\,\,\,\,{{\left( {{J_{{E_3}}}} \right)}_{22}} = 0.} }
Thus
{\rm{det}}({J_{{E_3}}}) = {{\lambda (1 - {R_3})({R_3} + 2h)} \over {2({R_3} + h)}}.
Note that there is a sign change in det(JE3) when R3 = 1, i.e., for μ = μ†. More precisely, det(JE3) is positive for 0 < μ < μ† and negative for μ† < μ. Thus, for μ > μ† the eigenvalues of JE3 are real numbers with opposed signs and E3 is a saddle, therefore unstable, proving (a).
Statements (b) and (c) are proven simultaneously. Indeed, if 0 < μ < μ†, follows that det(JE3) > 0, then for the stability analysis of E3 it is sufficient to study the sign of
{\rm{tr}}\left( {{J_{{E_3}}}} \right) = {{\lambda {R_3}} \over {2({R_3} + h)}}( - 3{R_3} + 1 - 2h).
Let P(μ) := −3R3 + 1 − 2h. Substituting the explicit form of R3 from (6) we have
P(\mu ) = - 3{\mu \over 2}(\mu + \sqrt {{\mu ^2} + 4h} ) + 1 - 2h.
Note that for
h > {1 \over 2}
we have that P is negative.
Now, take
0 \le h < {1 \over 2}
. Note that P is decreasing with respect to μ, so the sign of P is determined by its root μ*, (11). Thus, tr(JE3) is positive for μ < μ* and negative conversely. Therefore, E3 is stable if (11) holds and unstable if (12) is satisfied.
In Table 1 we summarize the behaviour of the equilibrium points of (4).
Behaviour and conditions of feasibility and stability of equilibria for the model (4), with 0 < h < 1, μ > 0, and μ* and μ† respectively defined in (11), (10). Note that if μ > μ†, E3 is unfeasible, having negative coordinates, although without direct biological meaning, this equilibrium collides with E2 when it becomes feasible, going through a transcritical bifurcation.
Equilibria
Feasibility
Stability
E1
always
unstable (saddle)
E2
always
unstable (saddle) if μ < μ†
stable if μ > μ†
E3
0 ≤ μ ≤ μ†
unstable (saddle) if μ > μ†
stable (node/focus) if max{0, μ*} < μ < μ†
unstable (focus) if 0 < μ < μ*
Bifurcations
Using Sotomayor’s Theorem [16] we prove the analytical conditions required for a transcritical bifurcation between E2 and E3 to occur. Also, using well known conditions (see, for example [8, 17]) we have shown that a Hopf bifurcation occcurs in E3 at the threshold value μ*, giving rise to sustained population oscillations [14, 15, 20]. Point E2 cannot give rise to oscillations, since both eigenvalues are always real. Instead, point E3 goes through a supercritical Hopf bifurcation, creating a stable limit cycle.
Proposition 4
In the system(4)we take μ as the bifurcation parameter and obtain:
A transcritical bifurcation between equilibria E2 and E3 when μ passes through the critical value μ†, (10).
A stable limit cycle that bifurcates at μ* (11) in a supercritical Hopf bifurcation.
Proof
(a) The equilibrium point E2 has the same coordinates as equilibrium E3 at the parametric threshold μ†. At this threshold the Jacobian of the system, evaluated at E2 = E3 becones
{J_{{E_2}}}\left( {{\mu ^\dagger }} \right) = \left( {\matrix{ { - \lambda } \hfill & { - 1} \hfill \cr \hfill 0 \hfill & 0 } } \right),
the zero eigenvalue is associated with the corresponding right and left eigenvectors: V = φ(1,−λ)Tand Q = ψ(0,1)T, where φ and ψ are arbitrary nonzero real numbers. Differentiating partially the right-hand sides of the system (4) with respect to μ we find
{f_\mu } = \left( {\matrix{ {{{RF} \over {{\mu ^2}\sqrt {R + h} }}} \cr { - {{RF} \over {{\mu ^2}\sqrt {R + h} }}} } } \right).
Now, calculating the Jacobian matrix of (4) in which the elements of this matrix are differentiated with respect to μ and then evaluated at E2 and μ† we get
D{f_\mu } = \left( {\matrix{ 0 \hfill & \hfill {\sqrt h } \hfill \cr 0 \hfill & { - \sqrt h } } } \right).
To test if Sotomayor’s Theorem conditions (see [16]) are satisfied and thus the existence of a transcritical bifurcation is guaranteed, we calculate D2f ((R,F); μ)(V,V) defined by
\left( {\matrix{ {{{{\partial ^2}{f_1}} \over {\partial {R^2}}}\xi _1^2 + 2{{{\partial ^2}{f_1}} \over {\partial R\partial F}}{\xi _1}{\xi _2} + {{{\partial ^2}{f_1}} \over {\partial {F^2}}}\xi _2^2} \hfill \cr {{{{\partial ^2}{f_2}} \over {\partial {R^2}}}\xi _1^2 + 2{{{\partial ^2}{f_2}} \over {\partial R\partial F}}{\xi _1}{\xi _2} + {{{\partial ^2}{f_2}} \over {\partial {F^2}}}\xi _2^2} } } \right),
where
{f_1} = {{dR} \over {dt}}
,
{f_2} = {{dF} \over {dt}}
of (4) and ξ1, ξ2 are the components of eigenvector V.
Calculating D2f and we obtain that:
\eqalign{ & {{Q^T}{f_\mu }({E_2};{\mu ^\dagger }) = 0,\,\,\,\,\,\,{Q^T}D{f_\mu }({E_2};{\mu ^\dagger })V = \varphi \psi \lambda \ne 0,} \hfill \cr & {{Q^T}{D^2}({E_2};{\mu ^\dagger })(V,V) = - {{2{\varphi ^2}\psi \lambda (2h + 1)} \over {2h + 2}} \ne 0.} }
Thus, all the conditions for the existence of a transcritical bifurcation are satisfied.
(b) The Jacobian matrix (7) evaluated at E3 is
{J_{{E_3}}} = \left( {\eqalign{& {{{\lambda \left( {{\mu ^2}( - 3{\mu ^2} - 8h + 1) + \mu \sqrt {{\mu ^2} + 4h} ( - 3{\mu ^2} - 2h + 1)} \right)} \over {2\mu \sqrt {{\mu ^2} + 4h} + 2{\mu ^2} + 4h}} \,\,\,- 1} \hfill \cr& {{{\lambda (\mu \sqrt {{\mu ^2} + 4h} (2 - 2{\mu ^2} - 4h) + {\mu ^2}(2 - 2{\mu ^2} - 8h) + 8h} \over {4(\mu \sqrt {{\mu ^2} + 4h} + {\mu ^2} + 2h)}}\,\,\,\,\,\,\,0} } } \right).
which yields the characteristic equation: Λ2 − tr(JE3)Λ + det(JE3) = 0 resulting in the pair of imaginary eigenvalues:
\matrix{ {{\Lambda _ \pm } = {{\tau \pm \sqrt {Q(\mu )} } \over 4},\,\,\,\tau = {\rm{tr}}({J_{{E_3}}}),} \hfill \cr {Q(\mu ) = {\tau ^2} - 4D,\,\,\,D = {\rm{det}}({J_{{E_3}}}) = {{\lambda (1 - {R_3})({R_3} + 2h)} \over {2({R_3} + h)}}.} }
Given that τ(μ*) = tr(JE3(μ*)) = 0, we have that Q(μ*) < 0, because
D = {\rm{det}}({J_{{E_3}}}({\mu ^*})) = {{\lambda (4h + 1)} \over 3}
is positive. Then, since Q(μ) is continuous, there is an interval T = (μ*− ε, μ* + ε) around μ*, such that, Q(x) < 0 whenever x ∈ T. The transversality condition
({d \over {d\mu }} \tau\mu))| {_{\mu = {\mu ^*}} \ne 0}
is satisfied because
{d \over {d\mu }}\tau ({\mu ^*}) = {{3\lambda (h + 1)(2h - 1)} \over {(4h + 1)\sqrt {3(h + 1)} }} > 0,
since
h > {1 \over 2}
or negative if
h < {1 \over 2}
.
Moreover, the Hopf bifurcation at equilibrium E3 is supercritical. Letting x = R − R3 and y = F − F3, the third-order Taylor series expansion for (4) evaluated at E3 and at the critical value μ* is given by
P(x,y) = S_R^*x + S_F^*y + {1 \over {2!}}(S_{RR}^*{x^2} + 2S_{RF}^*xy) + {1 \over {3!}}\left( {S_{RRR}^*{x^3} + 3S_{RRF}^*{x^2}y} \right)
and
Q(x,y) = H_R^*x + H_F^*y + {1 \over {2!}}(H_{RR}^*{x^2} + 2H_{RF}^*xy) + {1 \over {3!}}(H_{RRR}^*{x^3} + 3H_{RRF}^*{x^2}y)
where the asterisks denote the derivatives of
S(R,F) = \lambda R(1 - R) - {{RF} \over {\mu \sqrt {R + h} }}
and
H(R,F) = - F + {{RF} \over {\mu \sqrt {R + h} }}
evaluated at E3 and the critical value μ*.
Using the values for the Jacobian (13), we obtain from (14) and (15)P(x,y) = x′ = −y + p(x,y), Q(x,y) = y′ = ω2x + q(x,y), where
\omega = \sqrt {{J_{21}}} = \sqrt {{{\lambda (4h + 1)} \over 3}}
and
p(x,y) = {1 \over {2!}}(S_{RR}^*{x^2} + 2S_{RF}^*xy) + {1 \over {3!}}(S_{RRR}^*{x^3} + 3S_{RRF}^*{x^2}y)q(x,y) = {1 \over {2!}}(H_{RR}^*{x^2} + 2H_{RF}^*xy) + {1 \over {3!}}(H_{RRR}^*{x^3} + 3H_{RRF}^*{x^2}y).
In addition using a suitable changes of variables u = ωx and v = y, multiplying
\dot x
by ω, we obtain a “normal” form for the Hopf bifurcation [8]:
\dot u = - \omega v + \psi (u,v),\,\,\,\,\,\,\dot v = \omega u + \rho (u,v),
where the frequency of the limit cycle is given approximately by ω and
\psi (u,v) = {{S_{RR}^*} \over {2\omega }}{u^2} + S_{RF}^*uv + {{S_{RRR}^*} \over {6{\omega ^2}}}{u^3} + {{S_{RRF}^*} \over {2\omega }}{u^2}v\rho (u,v) = {{H_{RR}^*} \over {2{\omega ^2}}}{u^2} + {{H_{RF}^*} \over \omega }uv + {{H_{RRR}^*} \over {6{\omega ^3}}}{u^3} + {{H_{RRF}^*} \over {2{\omega ^2}}}{u^2}v.
Finally, according to the analytic criterion provided by [8], we have a supercritical Hopf bifurcation because φ < 0 where
\eqalign{ & {\phi = {1 \over {16}}({\psi _{uuu}} + {\psi _{uvv}} + {\rho _{uuv}} + {\rho _{vvv}}) + {1 \over {16\omega }}\left[ {{\psi _{uv}}({\psi _{uu}} + {\psi _{vv}}) - ({\rho _{uv}}({\rho _{uu}} + {\rho _{vv}})} \right.} \cr &\,\,\,\,\, {\left. { - {\psi _{uu}}{\rho _{uu}} + {\psi _{vv}}{\rho _{vv}})} \right] = - {1 \over {16}}{{756{h^2} + 540h + 27} \over {64{h^4} + 160{h^3} + 132{h^2} + 40h + 4}} < 0} }
if h > 0, the subscripts denote partial derivatives evaluated at (0,0).
Considering the parameters h and μ, the transcritical and Hopf bifurcation curves are shown in the Figure 2-a). These curves divide the square (1,0) × (1,0) in three specific regions determined by the stability of equilibrium points E2 and E3. A natural question is what kind of relationship there is between the Hopf bifurcation and the herd behavior effect, it means, if the Hopf bifurcation occurs in the presence or in the absence of the group defense effect; or in both cases. Effectively, the Figure 2-b) answers this question by a geometric analysis. Note that the curve
C = \{ {(h,\mu ),{\mu \over 2} ( {\mu + \sqrt {{\mu ^2} + 4h} } ) = h} \}
determines two regions in the parameter plane (h, μ): with and without group defense, (4) and (5) respectively. Since the curves H and C are transversal, it follow that the Hopf bifurcation occurs both in the presence or in the absence of the group defense effect.
Fig. 2
(a) The transcritical and Hopf bifurcation curves, respectively,
T = \{ {(h,\mu ),\mu = {\mu ^\dagger } = {1 \over {\sqrt {h + 1} }}} \}
and
H = \{ {(h,\mu ),\mu = {\mu ^*} = {{1 - 2h} \over {\sqrt {3(1 + h} )}}} \}
, determine three regions in the parameter plane (h, μ). The equilibrium point E3 is unstable in (1) “saddle” and (3) “focus or node”, and stable in (2) “focus or node”. The equilibrium point E2 is stable in (1) “node” and unstable in (2) and (3) “saddle”. (b) The curve C = {R3 = h} determines two regions in the parameter plane (h, μ). Namely, (4) region of group defense and (5) region without group defense.
Even though the curve C represents a threshold for effective defense, the transition between the regions (4) and (5) is smooth, since the response function f given by (1) is analytic. In [2] the authors perform a study of herd behaviour in which the response function is continuous but not smooth at critical threshold of group size for effective defense.
Numerical simulations
The numerical simulations here reported illustrate graphically the bifurcations obtained in the previous section. Since the empashis is more on the qualitative than on the quantitative results, the free software Geogebra is used to integrate the differential system and to draw the trajectories. Figure 3, illustrates the transcritical bifurcation between E2 and E3 when μ crosses the critical threshold μ†.
Fig. 3
The transcritical bifurcation between the equilibrium points E2 and E3 as function of the bifurcation parameter μ. a) For μ = 1.05, λ = 0.7 and h = 0.15, E2 is a stable node and E3 is a saddle. b) E3 = E2 is a nonhyperbolic equilibrium for μ = μ† ≈ 0.932, λ = 0.7 and h = 0.15. c) For μ = 0.88, λ = 0.7 and h = 0.15, E2 is a saddle and E3 is a stable node.
If we further decrease the parameter past the transcritical bifurcation at μ = μ†, E3 changes from a stable node to a stable focus, then undergoing a Hopf bifurcation at μ = μ* where it loses its stability, giving rise to a stable limit cycle. Figure 4 illustrates the bifurcation at
{\mu _2} = {\sqrt {33} ^{ - 1}}
for the chosen parameters values.
Fig. 4
Taking μ as a variation parameter a limit cycle arises from equilibrium E3R, indicating the onset of a Hopf bifurcation. a) For μ = 0.65 > μ*, λ = 0.7 and h = 0.15, E3 is a stable focus. b) For μ = μ* ≈ 0.376, λ = 0.7 and h = 0.15, E3 is a weak stable focus. c) For μ < μ*, λ = 0.7 and h = 0.15 E3 is an unstable focus.
Bilogical discussion
In order to create some biological interpretations, let’s first try to establish a meaning for the non-dimensional parameters used.
Parameter h has a quite direct interpretation, since it is simply the fraction of the threshold value
\tilde h
over K the carrying capacity of the prey population. So h << 1 implies that even small groups are beneficial to the prey population and display group defense, while h ≈ 1 means that only large groups display effective group defense. Intermediate values define the threshold, in comparison to the carrying capacity, beyond which group defense is effective.
Parameter λ = r/m can be interpreted in term of a comparison between two time scales: average time for prey reproduction and average length of the life of a predator. Assuming simplifying Poisson processes for the reproduction of prey and mortality of predators, 1/m represents the average lifespan of a predator and 1/r the time it takes for one prey to generate a descendant. So λ = (1/m)/(1/r) is the ratio between those two time scales, λ >> 1 means that predators have a long lifespan when compared with the reproductive cycle of the prey species while λ << 1 represents the opposite situation. It’s interesting that the most important qualitative behaviours do not depend on λ, but on the non-dimensional parameters h and μ.
Of course, such observation must not be taken too lightly, because λ does control the speed in which the dynamics take effect and that, depending on the time scale of analysis can make a great difference. Also, this parameter controls the coordinates of the equilibrium value E3 (co-existence) for the predator population, with high values corresponding to larger populations of predators. Finally, for the same initial condition, changing the value of λ can modify drastically the observed trajectory in a finite time scale. This because both the coordinates of the equilibrium change, and also the speed in which the trajectories tend to an stationary point or a stable limit cycle.
Naturally, for low values of λ and certain initial conditions, it is possible to observe trajectories that come very near the extinction of predators, which is expected, since it would mean that predators are dying too quickly while prey reproduce slowly. This result is not absolute, it depends on initial conditions and also relations between other scales, but, keeping fixing other factors, low values of λ imply in a worse situation for the survival of predators. Finally, while proving that the Hopf bifurcation is supercritical, we wrote where
\omega = \sqrt {{{\lambda (4h + 1)} \over 3}}
, where ω is the frequency of the limit cycle near the bifurcation. This means that λ controls the period of each cycle and, thus, the speed of the dynamics. Low values of λ tend to lead to slow cycles, while high values increase the speed of the dynamics as expected from our previous observations.
Analogously, parameter μ can be interpreted also as a ratio between two time scales. The first is LP = 1/m, which we already mentioned, and can be seen as an average lifespan of predators. The second, given by
{R_P} = 1/(ae\sqrt K )
can be interpreted as the average time that a predator would take to generate an offspring, if the prey population was fixed in the carrying capacity. Then, μ = RP/LP, so it is reasonable to expect extinction of predators if μ > 1. Indeed, the threshold for the stability of E2 = (1,0) is
{\mu ^\dag } = 1/\sqrt {h + 1} < 1
(stable if μ > μ†), meaning that E2 is always stable if μ > 1 and the predators tend to extinction.
It is interesting to note the effect of the group defense on the predator population. If we consider a fixed value for μ, which then could be interpreted of how favorable the environment is for the predators, the threshold μ† is the limit for the survival of the predators, since it is a decreasing function of h, the stronger the group defense effect (smaller h) the better it is for the survival of predators, which may come as a surprise. Higher values of h (meaning that only very large groups of prey flock together efficiently), lead to lower values of μ†, making it harder for the survival of the predator population.
One possible interpretation of this initially counter-intuitive result is that the group defense effect avoids over exploitation of the prey population, which would eventually lead to predator extinction. The transcritical bifurcation observed in the model is simply the change of stability in E2, when the environment becomes favorable enough for the establishment of the predator population.
As Figure 2 indicates, the Hopf bifurcation requires both relative low values for both h and μ, as defined by the curve
{\mu ^*} = (1 - 2h)/\sqrt {3(h + 1)}
. We could interpret the necessary conditions as that the predators must have relatively long lifespans and even small groups of prey can display group defense. In Figure 2-b) we also portray the curve
\mu = \sqrt {h/2}
which displays the set where R3, the prey coordinate for the coexistence equilibrium, is exactly equal to h the threshold population for group defense. We observe that the Hopf bifurcation occurs on both sides of the curve meaning that both when the group defense is fully operative or when it is not we may observe sustained oscillations. In fact, for certain combinations of parameters (μ = 0.05, h = 0.4 and λ = 2, for instance) the oscillations remain confined in the region where the prey population is always smaller than h.
On the other side, the necessary condition that h < 1/2 and μ < μ* do indicate that even though weak, the effect generated by the group defense is important in the establishment of the sustained oscillations. So it is seem reasonable to state that the oscillations may be interpreted as a tug-of-war between a situation less favorable for the prey (where population is too small to form group defense) and an alternative scenario (where population levels are closer to the ones necessary to establish an efficient group defense). The slow decaying predator population helps to create the necessary conditions for transitions between those regions. In the classic Lotka-Volterra model, the addition of intra-specific competition in prey is enough to dampen the oscillations, while in this case, with the group defense effect maintain the oscillations.
Conclusion
Our findings indicate that the proposed model shows a behaviour similar to the one found in [1]. In particular it gives rise to stable populations limit cycles. Ecologically this means that the suggested response function may be adequate if we want to model the prey herd behaviour that takes place only for a sizeable population, namely when the population level settles above a certain threshold.
It is also worth to note that the features obtained by using f(R), (1), as the response function are more realistic than those exhibited by the simpler
a\sqrt R F
. On the other hand, the increased mathematical complexity in the formulation of f(R) does not require a much more complicated analysis or shows more complex behavior than those of the system investigated in [1]. On the contrary, the new response function provides a beautiful example of a prey-predator system with fairly simple bifurcations. However, while the dynamics proposed in [1] in suitable circumstances induces total ecossytem collapse in finite time, which is a rather peculiar if not unique feature in continuous dynamics, the new proposed model displays a more regular behaviour, in which the populations in the exctinction scenarios possibly vanish exponentially fast, but not in finite time.
