Fear-driven extinction and (de)stabilization in a predator-prey model incorporating prey herd behavior and mutual interference

1.   Introduction

Interactions among predator and prey population species are modeled by systems of differential equations, and the functional response (number of prey consumed per predator per unit of time) of predators toward the prey is one of the important ecological components, which provides a bedrock for predator-prey dynamics. Mathematical models incorporating a functional response originated from the investigation of chemical reactions and biological interactions ^[1,2]. A large body of scholarly literature has shown that the functional response of predator can have profound impacts on the dynamics in natural predator-prey communities ^{[3,4,5,6,7,8,9,10,11]}. Thus, to make mathematical models more realistic, an appropriate choice of functional response is needed.

Herd behavior refers to the phenomenon in which individuals in a group act collectively for a given period without coordination by a central authority ^[12,13,14]. This behavioral phenomenon has been widely researched across several disciplines. For example, in early economics, Veblen studied herd behavior in sudden shifts in consumer behavior, such as fads and fashions ^[15]. Again, such phenomenon is seen in the Cobb-Douglas production function in the econometrics literature ^[16]. Its effect on population dynamics can be modeled using a functional response. Prey herd behavior is a form of anti-predator behavior and provides protection for the prey species against predators ^[17,18]. One way of modeling prey heard behavior is by using the functional response, $\varphi(u) = cu^p$, with $0 < p < 1$ introduced by Rosenzweig ^[19], where $u = u(t)$ is the population density of the herd and $c$ is the predation rate. Symbiotic, competition and predator-prey models in which the interaction terms use square root of the density of one population was considered by Ajraldi et al. ^[17]. Furthermore, Braza ^[20] studied a predator-prey model with the square root functional response $\varphi(u) = cu^{1/2},$ proposed by Gauss ^[3], implying a strong herd behavior where the predator interacts with the prey along the outskirts of the herd. An ecoepidemic predator-prey model with feeding satiation showing prey herd behavior and abandoned infected prey was investigated by Kooi and Venturino ^[21]. For other ways of modeling prey herd behavior, see ^[6,10] and references therein.

In addition, finite time extinction (FTE) of species exists in ecosystems and it is a significant issue for the management of natural resources ^[19,22]. The functional response $\varphi(u) = cu^p$, for $p < 1$ is non-smooth for $u = 0$ and therefore possess interesting complex dynamics. This response function allows for the extinction of the prey species in finite time, after which the predator population species exponentially decay to zero in infinite time ^{[18,20,23,24,25]}. Non-smooth functional responses (or power incidence functions) have been analyzed in susceptible-infective models, where the host species can potentially go extinct in finite time ^[26].

The need to consider intra-specific behavioral interactions among predators when searching for prey is a vital question for ecologists and conservationist trying to ascertain the dynamics that inform ecosystem balance. These behavioral effects, also known as mutual/predator interference impede the predators' searching efficiency as the density of the predators increases ^{[27,28,29,30]}. Several studies have concluded that mutual interference has a stabilizing effect on population dynamics, see ^[31] and references therein. Freedman and Wolkowicz ^[32] investigated the survival or extinction of predators in a deterministic predator-prey system exhibiting prey group defence. In this study, they determined that extinction due to group defence combined with enrichment can be prevented by introducing mutual interference of predators. For further discussions of mutual interference with other types of functional response, see ^{[33,34,35,36,37,38]}.

Recently, the non-consumptive effects of predation due to fear of predators has become the subject of interest for ecologists and mathematical biologists. Experiments on terrestrial vertebrates showed that the presence of a predator may play an important role by changing the behavior of the prey demography ^[39]. Zanette et al. ^[39] manipulated predation risk of song sparrows for the duration of an entire breeding season. This experiment was conducted to ascertain whether perceived predation risk alone could have an impact on the reproduction of the song sparrows. Suraci et al. ^[40] observed from experimentation that the effect of the manipulation of fear of large carnivores causes a tropic cascade. Hua et al. ^[41] studied how increased perception of predation risk to adults and offspring alters reproductive strategy and performance. In Wirsing and Ripple ^[42], a comparison of shark and wolf research revealed similar behavioral responses by prey. Bauman et al. ^[43] investigated how the effects of fear associated with predator presence and habitat structure interact to change the removal of macroalgal biomass (i.e herbivory) on coral reefs. They observed that the effects of fear due to the presence of predators were highest at low macroalgal density, but lost at higher densities due to increased background risk.

Motivated by these ecological and biological findings, Wang et al. ^[44] introduced a mathematical model incorporating fear. The authors demonstrated that strong fear responses can have a stabilizing effect on a predator-prey model with Holling type Ⅱ response by excluding periodic solutions to the system, resulting in a locally stable point of coexistence between the predator and prey populations ^[44]. Subsequent studies investigated the dynamics of fear in models incorporating hunting cooperation, prey refuge, Leslie-Gower type and a variety of functional responses, such as Beddington-DeAngelis, Holling type Ⅰ, Ⅱ, Ⅲ and Ⅳ ^{[45,46,47,48,49,50]}. A recent work considering the effect of mutual interference and fear on a predator-prey model with a Holling type Ⅰ functional response established that the inclusion of mutual interference promotes system stability ^[51].

The qualitative effect of predator harvesting on the stability of the ecosystem has been investigated extensively ^{[52,53,54,55,56,57]}. Chakraborty et al. ^[58] explored a mathematical study with biological ramifications of a predator-prey model with predator effort harvesting. Their result suggests that harvesting of predator may be one of several ways to observe coexistence of prey and predator population species in the laboratory study and possibly nature. The ratio-dependent predator-prey model where the predator population is harvested at catch-per-unit-effort hypothesis is investigated by Gao et al. ^[59]. Therein, they studied the temporal, spatial and spatiotemporal rich dynamics due to the non-smoothness of the origin.

Our primary contributions in the present manuscript are:

1. We formulate a mathematical model (i.e. model (2.2)) incorporating the combined effects of fear of predator, prey herd behavior and mutual interference.

2. We study the effect of fear of predators on the dynamics of the model (2.2). We note that when there is no fear (i.e. $k = 0$), there is some initial data that converges uniformly to a stable coexistence equilibrium point. With the introduction of fear of predator (i.e. $k > 0$), that same initial data will converge to the predator axis in finite time. This phenomenon is shown analytically via Theorem 5.1 and presented numerically in Figure 6.

3. We analyze the impact of predator harvesting on the dynamics of the modified Lotka-Volterra model with fear effect. Our mathematical conjecture (see Conjecture 4) and numerical simulation (see Figure 8) reveal that, harvesting of predators can prevent finite time extinction of the prey species.

This paper is arranged as follows: In Section 2, we propose a mathematical model of systems of differential equations to incorporate the combined effects of fear of predators, prey herd behavior and mutual interference. Guidelines to dynamical analysis are presented in Section 3, where we investigated the possible existence of biologically feasible equilibrium points and the stability of the coexistence equilibrium. In Section 4, we derive conditions for the existence of local and global bifurcations including saddle-node, Hopf, Bautin, and homoclinic bifurcations. Finite time extinction of the prey species driven by fear of predators were analyzed in Section 5. In Section 6, we investigate the effect of effort harvesting of predators. To illustrate the feasibility of our mathematical analysis and conjectures, extensive numerical solutions are presented in Section 7. The paper ends with concluding remarks in Section 8.

2.   The mathematical model

Consider a modified Lotka-Volterra predator-prey model with predation intensity and mutual interference of predators. Let $u(t)$ and $v(t)$, respectively, denote the prey and predator population densities at any time $t$. The model is given by the following systems of equations

where $0 < m, p\leq 1$. Let $m$ denote mutual interference parameter introduced by Hassell ^[28], $1/p$ is the intensity of predation and $p$ determines the slope of the functional response at the origin, $a$ denotes the birth rate of prey, $d$ denotes the death rate of predator, $b$ reflects the intraspecific competition of the prey, $c$ denotes the rate of predation, and $e$ measures efficiency of biomass conversion from prey to predator. When $p = m = 1$, the model (2.1) degenerates to the classical Lotka-Volterra model ^[1,2]. All parameters are assumed to be positive. The underlying assumptions of the model (2.1) are as follows:

(ⅰ) The first equation in model (2.1) describes the change in prey population with respect to time, and it is separated into three parts, namely birth rate, effect of the density of one species on the rate of growth of the other and functional response of the predator towards the prey.

(ⅱ) The second equation in model (2.1) describes the change in predator population with respect to time and it is separated into two parts, namely death rate, $d$, the predators die out in the absence of its only food source, prey and biomass conversion from prey to predator with rate $e$.

(ⅲ) The term $v^m$ models the intra-specific behavioral interactions among predators when searching for prey. For $m < 1$, there is predator interference, where larger predator densities leads to less consumption per capita. Furthermore, this leads to a nonvertical predator nullcline.

(ⅳ) The predator is consuming the prey with the functional response $\varphi(u) = u^p$, for $0 < p < 1$ (see ^[10] for assumptions of $\varphi$).

Model (2.1) has been well investigated in infectious disease modeling, where $u$ represents the density of susceptible populations, $v$ represents the density of infective populations and the term $u^pv^m$ (i.e. modified Lotka-Volterra interaction term) represents power incidence function ^[26]. As we have seen, the model considered in ^[26] reveals some significant and interesting results due to the power incidence function. A natural question that arises is: how does the effect of fear of predators affect the dynamical behaviors of the model (2.1)? Does it stabilize, destabilize or have no influence?

Now, based on experimental evidence ^[39], we assume fear of predators decreases the birth rate of the prey species. To account for the decrease in the prey population due to fear of predators, the birth rate of the prey is multiplied by the term $\phi(k, v) = \frac{1}{1+kv}$, introduced by Wang et al. ^[44], which is monotonically decreasing in both $k$ and $v$. Here $k$ denotes the strength of fear of predator. Biologically, it is appropriate to assume the following:

To the best of our knowledge, there does not exist any scholarly literature that investigates the combined influence of fear of predators on predator-prey interactions with prey herd behavior and mutual interference. This motivates us to formulate the following model

When $p = m = 1$, we recover the results from Wang et al. ^[44]. Moreover, for $p = 1$ and $0 < m\leq 1$, we recover results from Xiao and Li ^[51]. Recently, Fakhry and Naji ^[60] investigated the model (2.2), where the fear function was multiplied to the logistic growth term i.e. $(au-bu^2)\phi(k, v)$ with square root functional response (i.e. $\; p = 0.5$) and no predator interference (i.e. $m = 1$). Huang and Li ^[61] disproved and also provided an alternative proof for some of the results obtained by Fakhry and Naji. Our model provides a generalization of the models mentioned above, and we will focus on the case where $0 < m, p < 1$.

3.   Dynamical analysis

The dynamical analysis of the model (2.2) is investigated in this section.

Lemma 3.1. Consider the first quadrant $\mathbb{R}_+^2 = \{(u, v):u\geq 0, v\geq 0\}$, then the solutions $(u(t), v(t))$ of model the (2.2) which initiate in $\mathbb{R}_{++}^2$ are nonnegative for all $t\geq 0$. Here, $\mathbb{R}_{++}^2 = \{(u, v):u > 0, v > 0\}$.

Proof. The right hand side of model (2.2) is continuous and locally non-smooth in $\mathbb{R}_{+}^2$. Also, the solution $(u(t), v(t))$ which initiate in $\mathbb{R}_{++}^2$ of model (2.2) exists and is non-unique. From model (2.2), we obtain

However, $v$ stays positive for all $t > 0$.

In theoretical biology and ecology, nonnegativity of the model (2.2) implies survival of the populations over some temporal domain.

Lemma 3.2. The solutions $(u(t), v(t))$ of model the (2.2) which initiate in $\mathbb{R}_{++}^2$ are uniformly bounded and dissipative.

The proof of Lemma 3.2 is standard and therefore omitted in this work.

The boundedness of a system limits total population growth of the interacting species, ensuring that neither population experiences exponential growth over a long time interval. As a condition of this property, total population values will not reach impracticable quantities in a period of time. Also, in a dissipative model, the population of each species is bounded from above for all time. This guarantees that the individual populations of the predator or the prey do not exceed a finite upper limit.

3.1.   Existence of equilibria

We present analytic guidelines in this subsection to analyze the model (2.2) and to investigate its equilibria. Consider the solutions to the steady state equations:

The above equations contain the following non-negative equilibrium points.

(ⅰ) The trivial (extinction) equilibrium point $E_0(0, 0)$, always exist. Here, neither the prey nor predator population survives in the ecosystem.

(ⅱ) An axial (predator free) equilibrium point $E_1(a/b, 0)$, always exist. Here, the predator population goes into extinction and only the prey population survives.

(ⅲ) The coexistence equilibrium point(s) $E_2(u^*, v^*)$. Ecologically, this equilibrium point is important since both prey and predator populations coexist.

The existence of a coexistence equilibrium point is ascertained by finding the intersection(s) of the prey and predator nullclines.

First, the prey nullcline is determined by the equation

If $u = 0$, we obtain

But $v = 0$ implies $\sigma(v) = a\neq 0$, since $a$ is a positive constant. Furthermore,

and thus, there exist some $v > 0$ such that $\sigma(v) = 0$. Also, if $v = 0$ in equation (3.1), then $u = \frac{a}{b} > 0$.

Moreover, we assume that

and observe that $u > 0, v > 0$ and

Now the graph of the prey nullcline is concave and intersects the prey axis at $u = 0$ and $u = \frac{a}{b}$.

Additionally, the predator nullcline is determined by the equation

Solving for $v$ in equation (3.3) yields

Clearly, the point $(0, 0)$ lies on the predator nullcline. By computing the first and second derivatives with respect to $u$, we obtain

From equation (3.5), the sign of the second derivative depends on $p+m-1$.

Hence, the predator nullcline goes through $(0, 0)$, and as $u$ increases, the predator nullcline increases monotonically. Now by the intermediate value theorem, the prey and predator nullclines will intersect in $\mathbb{R}_{++}^2$ to produce a unique (i.e. $E_2(u^*, v^*)$) or two (i.e. $E_2^i(u_i^*, v_i^*)$, for $i = 1, 2$ and $0 < u_1^* < u_2^* < \frac{a}{b}$) coexistence equilibrium points.

Remark 1. Since $0 < p, m < 1$, there is singularity in the Jacobian at $E_0$ and $E_1$. Hence we cannot analyze the stability of $E_0$ and $E_1$ by linearizing the model (2.2).

3.2.   Stability analysis at a coexistence equilibrium point

We discuss in this subsection the local stability at any coexistence equilibrium point.

Theorem 3.3. Consider the model given by (2.2).

(a) For $p+m\geq 1$, there exists a unique coexistence equilibrium point $E_2(u^*, v^*)$ which is locally asymptotically stable (LAS) by the Routh-Hurwitz criterion.

(b) For $p+m < 1$, either there exist two coexistence equilibrium points or none. However, if there exist two coexistence equilibrium points, then $E^1_2(u_1^*, v_1^*)$ is a saddle and $E^2_2(u_2^*, v_2^*)$ is LAS.

Proof. The linearized model (2.2) at any coexistence equilibrium point $(u^*, v^*)$ is given by the Jacobian Matrix ${\bf{J}}$

where

The characteristic equation at the coexistence equilibrium is

where

and

By using implicit function theorem as used in ^[62], we obtain

where $\frac{dv^{*(g_2)}}{du^*}$ and $\frac{dv^{*(g_1)}}{du^*}$ are the slopes of the tangents of the predator and prey nullclines at $E_2$ respectively.

(a) Now we assume $p+m\geq 1$. The two possibilities for the sign of the Jacobian matrix at $E_2$ are:

or

If $b > c(1-p)u^{*p-2}v^{*m}$, then $\det \, ({{\bf{J}}}) > 0$ and $\operatorname{tr} \, ({{\bf{J}}}) < 0$. Thus in (3.7), $E_2$ is LAS by Routh-Hurwitz criterion. In (3.8), the $\det \, ({{\bf{J}}}) > 0$ since $\frac{dv^{*(g_2)}}{du^*} > \frac{dv^{*(g_1)}}{du^*}$. This is clearly seen in Figure 1[(a) and (b)]. Thus, $E_2$ is LAS if the $\operatorname{tr} \, ({{\bf{J}}}) < 0$.

(b) Assume $p+m < 1$ and there exist two coexistence equilibrium points. The sign of the Jacobian matrix at $E^1_2$ is given by

and $E^2_2$ is

or

In (3.9), $b < c(1-p)u^{*p-2}v^{*m}$ and $\frac{dv^{*(g_2)}}{du^*} < \frac{dv^{*(g_1)}}{du^*}$, hence $E^1_2$ is a saddle point since the $\det \, ({{\bf{J}}}) < 0$. This is evident in Figure 1(c). In (3.10), $b > c(1-p)u^{*p-2}v^{*m}$, thus $E^2_2$ is LAS since $\det \, ({{\bf{J}}}) > 0$ and $\operatorname{tr} \, ({{\bf{J}}}) < 0$. In (3.11), when $b < c(1-p)u^{*p-2}v^{*m}$ then $\det \, ({{\bf{J}}}) > 0$ if $\frac{dv^{*(g_2)}}{du^*} > \frac{dv^{*(g_1)}}{du^*}$. Therefore, $E^2_2$ is LAS if $\operatorname{tr} \, ({{\bf{J}}}) < 0$.

4.   Bifurcation analysis

4.1.   Local bifurcation

In this subsection, we investigate the qualitative changes in the dynamical behavior of model (2.2) under the effect of varying the strength of the fear of predator $k$. The conditions and restrictions for the occurrence of saddle-node and Hopf bifurcations are derived analytically and their classification is of co-dimension $1$ bifurcations. Additionally, we present numerically the two-parameter projection of the Hopf-bifurcation diagrams of the model (2.2).

4.1.1.   Saddle-node bifurcation

Saddle-node bifurcation occurs when shifting a parameter value causes two equilibria of contrasting stability to collide and mutually disappear, forming an instantaneous saddle-node at the point of their collision. This is a one parameter bifurcation and hence of codimension 1. In the next theorem, we show that using the strength of fear as a bifurcation parameter, the model (2.2) satisfies the conditions for saddle-node bifurcation.

Theorem 4.1. Model (2.2) admits a saddle-node bifurcation around $E_2$ at $k_s$ when the model parameter values satisfy the conditions $\frac{dv^{*(g_2)}}{du^*} = \frac{dv^{*(g_1)}}{du^*}$ and $\operatorname{tr} \, ({{\bf{J}}}) < 0$.

Proof. In order to verify the conditions for the existence of saddle-node bifurcation, we employ Sotomayor's theorem ^[63] at $k = k_s$. At $k = k_s$, we obtain that $\frac{dv^{*(g_2)}}{du^*} = \frac{dv^{*(g_1)}}{du^*}$ and $\operatorname{tr} \, ({{\bf{J}}}) < 0$, which shows that the Jacobian $\, ({{\bf{J}}})$ has a zero eigenvalue. Let $W$ and $Z$ be the eigenvectors corresponding to the zero eigenvalue of the matrix $\, {{\bf{J}}}$ and $\, {{\bf{J}}}^T$ respectively. Here, $W = (w_1, w_2)^T$ and $Z = (z_1, z_2)^T$, where $w_1 = -\frac{J_{12}w_2}{J_{11}}$, $z_1 = -\frac{J_{21}z_2}{J_{11}}$ and $w_2, z_2 \in \mathbb{R} \setminus \{0\}$.

Furthermore, let $H = (F, G)^T$ and $\tilde{M} = (u^*, v^*)^T$, where $F, G$ are defined in (2.2). Thus

and

Therefore model (2.2) admits a saddle-node bifurcation when $k = k_s$.

Remark 2. Additionally, the model (2.2) undergoes saddle-node bifurcation around $E_2$ with respect to the following parameters, $d, b, a, c$ and $e$. See Figure 9 in the appendix for numerical verification.

4.1.2.   Hopf-bifurcation

Similar to a saddle-node bifurcation, a Hopf-bifurcation describes a local change in the stability of an interior equilibrium point due to an alteration of a parameter. However, for a Hopf-bifurcation, varying the bifurcation parameter does not annihilate or create new equilibrium points. Rather, at the point where system stability shifts (i.e. Hopf point) – a stable or unstable periodic orbit develops. This is a one parameter bifurcation and hence of codimension 1. The conditions for the existence of Hopf-bifurcation of the model (2.2) is derived in the theorem below.

Theorem 4.2. Model (2.2) experiences Hopf-bifurcation around the coexistence equilibrium point $E_2$ at $k = k_h$, where

when the following conditions are satisfied:

$S(k) = 0$, $M(k) > 0$, and $\dfrac{d}{dk}\left. \mathit{\text{Re}} [\eta_i (k) \right]|_{k = k_h}\neq 0$ for $i = 1, 2$.

Proof. Using the strength of fear as a bifurcation parameter, consider the Jacobian matrix (3.6) around the coexistence equilibrium $E_2$. The characteristic equation at $E_2$ is given by

where $S = \operatorname{tr} \, ({{\bf{J}}}) = J_{11} + J_{22}$ and $M = \det \, ({{\bf{J}}}) = J_{11} J_{22}- J_{12} J_{21}.$ The zeros of the equation (4.1) are

At $k = k_h$, $S(k) = 0$ implies

The characteristic equation (4.1) becomes

at $k = k_h$. Solving for the zeros of equation (4.3) yields $\eta_{1, 2} = \pm i\sqrt{M}$. Thus, a pair of purely imaginary eigenvalues. Furthermore, we substantiate the transversality condition. For any $k$ in the neighborhood of $k_h$ in (4.2), let $\xi(k) = \text{Re} \left[\eta_i (k) \right] = \frac{1}{2}S(k)$ and $\mu(k) = \sqrt{M(k)-\frac{[S(k)]^2}{4}}$. Thus,

Thus the tranversality condition is satisfied if $\dfrac{d}{dk}\left. \text{Re} [\eta_i (k) \right]\neq 0$ at $k = k_h$. Therefore, by the Hopf-bifurcation Theorem ^[64], the model (2.2) experiences a Hopf-bifurcation around $E_2$ at $k = k_h$.

4.1.3.   Direction of Hopf-bifurcation

We investigate the stability and direction of the periodic cycles emitted via Hopf-bifurcation around the coexistence equilibrium point by computing the first Lyapunov coefficient ^[63]. We first translate the coexistence equilibrium $E_2$ of the model (2.2) to the origin by using the transformation $x = u-u^*$ and $y = v-v^*$. Now, the model (2.2) becomes

Applying Taylor series expansion at $(x, y) = (0, 0)$ up to third order, we obtain the following planar analytic model

where

Since $a_{10}, a_{01}, b_{10}$ and $b_{01}$ are the components of the Jacobian matrix ${\bf{J}}$ evaluated at the coexistence equilibrium point $E_2$ at $k = k_h$, then $S = a_{10}+b_{01} = 0$ and $M = a_{10}b_{01}-a_{01}b_{10} > 0$.

The first Lyapunov coefficient $L$ ^[63] is computed by the formula

Now, if $L < 0$, then the Hopf-bifurcation is supercritical and subcritical if $L > 0$.

4.1.4.   Generalized Hopf-bifurcation

We note here that from the two-dimensional projection of the Hopf-bifurcation diagrams of the model (2.2) in Figure 4, a generalized Hopf-bifurcation or Bautin bifurcation is observed. The generalized Hopf-bifurcation is a local bifurcation of co-dimension 2 and this happens when the first Lyapunov coefficient is zero, and the coexistence equilibrium point has a pair of purely imaginary eigenvalues. The generalized Hopf-bifurcation point separates branches of subcritical and supercritical Hopf-bifurcation in the parameter plane.

Now, we present a conjecture that pertains to generalized Hopf-bifurcation.

Conjecture 1 (Existence of Generalized Hopf-bifurcation). Assume that the model (2.2) admits a Hopf-bifurcation as in Theorem 4.2 for a given set of parameters. If the first Lyapunov coefficient becomes zero and the coexistence equilibrium point has a pair of purely imaginary eigenvalues, then the model (2.2) undergoes a Bautin or generalized Hopf-bifurcation.

4.2.   Global bifurcation

Now, there exist a unique value of parameter $k$ for which $W^s(E_0)$ and $W^u(E_0)$ coincide, and that implies existence of a homoclinic loop in model (2.2). In particular, we state a conjecture concerning the effect of fear of predators on the global dynamics of model (2.2).

Conjecture 2 (Existence of Homoclinic Bifurcation). Consider the model (2.2) where all parameters are fixed except $k\geq 0$. There exists $k^* > 0$ such that a homoclinic loop occurs when $k = k^*$.

5.   Fear-Driven Finite Time Extinction (FDFTE)

We seek to investigate the effect of introducing fear of predator in the predator-prey model – that is, is it possible for the fear effect to drive a stable coexistence equilibrium point to extinction in finite time?

Thus, we state our result concerning finite time extinction driven by fear of predators.

Theorem 5.1 (FDFTE). Consider the predator-prey model given by (2.2), and a certain parameter set, and certain initial data $(u(0), v(0))$ that converges uniformly to a stable coexistence equilibrium point $(u^*, v^*)$ for $k = 0$. Then there exists $k > 0$ such that all trajectories initiating from the same initial data $(u(0), v(0))$ will lead to finite time extinction of $u$, followed by $v$ going extinct asymptotically.

Proof. We argue by contradiction. We begin by assuming not. Thus for a certain parameter set, with $k = 0$ and certain initial data $(u^{*}(0), v^{*}(0))$ that converges uniformly to a stable coexistence equilibrium point $(u^*, v^*)$, there exists a $k > 0$ s.t for trajectories emanating from the same initial data, and parameters, we will have

on $[0, T^{*}], \ \forall T^{*}$. Now,

This implies,

Note, via (5.2), the upper bound on $u$, and positivity of solutions, we have

However, we see that the solution to $u$,

will go extinct in finite time if

Thus, given an initial data $(u^{*}(0), v^{*}(0))$, and $T^{*} > 0$ we can choose $k > > 1$, s.t.

so that we obtain,

on $[0, T^{*}]$, with $u$ being driven to extinction in finite time, from which the asymptotic extinction of $v$ follows. Since $T^{*}$ is arbitrary, we have derived a contradiction, and so the theorem is proved.

Remark 3. We see from Theorem 5.1 that $k$ chosen sufficiently large will lead to prey extinction in finite time; however in practice, $k$ need not be large, see Figure 6. Thus, a necessary condition on the size of $k$ to cause finite time prey extinction remains unproven.

6.   Impact of effort harvesting of predators

In this section, we consider the impact of an external effort dedicated to harvesting of predators in the model (2.2). A natural question that arises is: how do the external effort of predator harvesting affect the FTE dynamic of the prey population species? Here, the harvesting is proportional to the density of harvested predator population species. The corresponding differential equations can be represented as:

Let the parameter $q > 0$ represents the external effort dedicated to predator harvesting. Here $0 < r\leq 1$. We recover the model (2.2) when $q = 0$. Nonnegativity of solutions of model (6.1) follows from Lemma 3.1. The model (6.1) contains a trivial equilibrium point $E_0(0, 0)$, an axial equilibrium point, $E_1(a/b, 0)$, and coexistence equilibrium point(s) $E_2$ (or $E^i_2$, for $i = 1, 2$).

6.1.   Dynamical Guidelines

The linearized model (6.1) at any coexistence equilibrium point $(u^*, v^*)$ is given by the Jacobian Matrix ${\bf{J^*}}$

where

$c_{11} = J_{11}, c_{12} = J_{12},$ and $c_{21} = J_{21}$. Here $J_{11}, J_{12}$ and $J_{21}$ are given by (3.6). Now

Theorem 6.1. Consider the model given by (6.1) and assume $r\geq m$.

(a) There exists a unique coexistence equilibrium point $E_2(u^*, v^*)$ which is LAS.

(b) There exist two coexistence equilibrium points such that $E^1_2(u_1^*, v_1^*)$ is a saddle and $E^2_2(u_2^*, v_2^*)$ is LAS.

Proof. The proof of Theorem 6.1 is similar to proof in Theorem 3.3 and therefore omitted.

Theorem 6.2. Model (6.1) experiences Hopf-bifurcation around the coexistence equilibrium point $E_2$ at $q = q_h$, where

provided $S(q) = 0$, $M(q) > 0$, and $\dfrac{d}{dq}\left. \mathit{\text{Re}} [\lambda (q) \right]|_{q = q_h}\neq 0$ for $i = 1, 2$.

Example 1. To validate Theorem 6.2, we consider the following parameter values $m = 0.6, \; a = 3, \; k = 0.08, \; b = 0.2, \; p = 0.5, \; d = 1, \; c = 2, \; e = 1.1, \; q = 1, \; r = 1$ (see Figure 7(a)). Hopf-bifurcation is obtained at $q = q_h^* = 0.27068$ around the coexistence equilibrium point $E_2(2.54542, 2.24175)$. Furthermore, at $q = q_h$, $S(q) = \operatorname{tr} \, ({{\bf{J^*}}}) = 0$, $M(q) = \det \, ({{\bf{J^*}}}) = 0.29264 > 0$ and $\frac{d}{dq}\left. \text{Re} [\lambda (q) \right]|_{q = q_h}\neq 0$. Hence, all necessary and sufficient condition for Hopf-bifurcation to occur are satisfied.

We now state two conjectures concerning the effect of effort harvesting of predators on the dynamics of model (6.1).

Conjecture 3 (Existence of Homoclinic Bifurcation). Consider the model (6.1) where all parameters are fixed except $q\geq 0$. There exists $q^* > 0$ such that a homoclinic loop occurs when $q = q^*$.

Conjecture 4 (Harvesting-Induced Recovery). Consider the predator-prey model given by (6.1). Then there exists a harvesting effort, $0 < c_{1} < q^* < c_{2}$ such that the solution to the prey equation does not go extinct in finite time, and in particular for these levels of predator harvesting the solution can be driven to a coexistence state, if initiated from certain initial conditions.

7.   Numerical simulations

In this section, we shall visualize the role of fear, herd behavior, mutual interference, and harvesting in the models (2.2) and (6.1). Indeed in model (2.2), the coexistence equilibria are not analytically accessible. Numerical simulations are provided as guidelines in Figure 1, to show the existence of a unique coexistence equilibrium point for $p+m\geq1$ and two coexistence equilibria for $p+m < 1$.

We obtained saddle-node bifurcations of the model (2.2) as we increase the strength of fear of predator with some appropriate choice of parameters, see Figure 2.

From the numerical simulations in Figure 3, we show that the fear of predator has an effect in altering the stability of the coexistence equilibrium solution $E_2(u^*, v^*)$ of the model (2.2) via Hopf-bifurcation. In Figure 3(a), the coexistence equilibrium solution changes from an unstable zone to a stable zone around $E_2(0.43392, 0.14327)$ as the strength of fear of predator $k$ crosses the Hopf point at $k_h^* = 15.093353$. We used MATCONT ^[65] to generate the bifurcation diagrams and obtained $L = 1.83219e^{-01} > 0$, hence subcritical Hopf-bifurcation. Furthermore, in Figure 3(b), the coexistence equilibrium solution changes from a stable zone to an unstable zone around $E_2(2.26530, 3.16304)$ as the parameter $k$ crosses the Hopf point at $k_h^* = 0.061382$ with $L = 1.76388e^{-02} > 0$, hence subcritical Hopf-bifurcation. In summary, we conclude here that with appropriate parameters the effect of fear of predator can have a stabilizing and destabilizing effect on the coexistence equilibrium solutions of the model (2.2).

From the two-dimensional projections of Hopf-bifurcation curves, we observed a generalized Hopf-bifurcations or Bautin bifurcation which is local and of co-dimension 2 (see Figures 4$(a)-(d)$). Next, we explain Conjecture 2 via numerical simulations in Figure 5. In Figure 5(a), $E_2$ is unstable and the unstable manifold ($W^u(E_0)$) surrounds the stable manifold ($W^s(E_0)$). In Figure 5(b), $E_2$ is stable and a homoclinic loop is formed as the strength of fear is increased.

Hopf-bifurcation occurs at $q = q_h^* = 0.27068$ around $(2.54542, 2.24175)$ with first Lyapunov coefficient $L_H = 8.38051e^{-03}$ in the model (6.1). The periodic orbits emitted at the Hopf point is subcritical. A two-parameter bifurcation diagram in $q-k$ parametric plane is shown in Figure 7(b). Clearly, increasing the harvesting effort $q$, will cause a lowering of the predator nullcline, bringing down the coexistence equilibrium, see Figure 8. In Figure 8(a), $E_2$ is unstable and all the points are attracted to the predator axis when there is no harvesting effort. A homoclinic loop is formed in Figure 8(b) when $W^s(E_0)$ coincide with $W^u(E_0)$ for $0.321 < q < 0.323$. Here $E_2$ is stable. Also, in Figure 8(c), for $q = 1$, $E_2$ is stable and all initial data below $W^s(E_0)$ converges to $E_2$ whilst those above are attracted to the predator axis.

For large values of $q$, this could be brought as close to the predator free equilibrium as desired, see Figure 7.

8.   Conclusions

In this paper, we have proposed and investigated the rich dynamical behavior of a predator-prey model (2.2), incorporating the effect of fear, prey herd behavior and mutual interference of predators. The prey herd behavior is governed by a modified Holling type-Ⅰ functional response that allows for finite time extinction of the prey species due to the non-smoothness at the trivial equilibrium point ^[19,20]. The stable manifold ($W^s(E_0)$) of the trivial equilibrium $E_0$ splits the phase plane into two, where solutions with initial conditions above the stable manifold are attracted towards the predator axis in finite time. This posses a problem of non-uniqueness of solutions in backward time.

The fear of predation risk excited by predators can drive a stable prey population species to extinction in finite time, and consequently, to the extinction of the predator population species. To this end please see Theorem 5.1. We provide numerical justification in Figure 6.

Taking the strength of fear of predator $k$ as a bifurcation parameter, we have shown analytically and numerically various local and global bifurcations. We observed from our investigation that fear of predator has the tendency to stabilize and destabilize a coexistence equilibrium point by producing limit cycles via subcritical Hopf-bifucation. Biologically, a strong strength of fear can stabilize an unstable coexistence equilibrium of interacting species, see Figure 3(a). Also, with weak strength of fear of predator and certain parameter sets, the stable coexistence between a predator and prey can be destabilized, see Figure 3(b). In Figure 3, the effects of fear reduces both predator and prey population densities. There exist a critical strength of fear where the stable and unstable manifold of the trivial equilibrium meet (i.e. homoclinic bifurcation), see Figure 5(b). This is conjectured in Conjecture 2. Here we observe that all solutions with initial conditions inside the loop goes to the stable coexistence equilibrium whilst those outside goes to prey extinction in finite time. Hence, these obtained results are interesting and provide further justification that fear of predation risk plays a crucial role in ecosystem balance ^[66].

We conjecture via Conjecture 4 that harvesting of predators can prevent finite time extinction of the prey species and yield persistence of the predator and prey population species. Proving this conjecture would make interesting future work. See Figure 8. Additionally, from our numerical simulations in Figure 8, when the unstable manifold of $E_0$ is above the stable manifold of $E_0$ coupled with an unstable coexistence equilibrium point, we observed a homoclinic loop by introducing an external effort dedicated to predator harvesting. Biologically, we are able to stabilize the predator and prey population species that initiated inside the loop with low harvesting effort. Also, when the external effort dedicated to predator harvesting rate is very high, the prey population species approach its carrying capacity and the predator population species get close to extinction, see Figure 7(a). Thus, the predator population species may not survive at a very high harvesting effort. Note, harvesting finds large scale applications in current bio-control applications ^[67]. A full dynamical analysis of (6.1) would make an interesting future endeavor, as the FTE dynamic in predator could counteract with the FTE dynamic in prey, to generate rich dynamical behavior.

Another bio-control application is in pest management where the fear of natural enemy is introduced to drive an invasive pest into extinction. This study should, therefore, prove to be a useful tool in resource management and control. A further interest, which the authors are currently investigating, is the interplay of fear of predator between aggregating prey species and predator interference models, such as those considered in ^[10,68,69].

Acknowledgments

We are grateful to the handling editor and the anonymous reviewers for their insightful comments and suggestions which led to the improvement of the quality of this work. KAF would like to express his gratitude for an internal Faculty Development Grant (FDG084) from Samford University.

Conflict of interest

The authors declare there is no conflict of interest in this paper.

Appendix

We provide numerical simulations to visualize the saddle-node bifurcation of parameters $b,c$ and $d$.