Fear effect in prey and hunting cooperation among predators in a Leslie-Gower model

1.   Biological background and formulation of the model

Dynamical complexities of interacting predator-prey model are extensively studied by several researchers to understand the long-time behaviour of the species. Most of the predator-prey models are based upon classical Lotka-Volterra formalism, where the prey consumption rate of predators is the growth rate of predators with a conversion factor. However, in this paper we confine ourselves with Leslie-Gower model, whose main feature is logistic type growth function for predators ^[1]. Leslie ^[2] introduced that the environmental carrying capacity ($K_P$) of the predators is proportional to the number of prey abundance (N) i.e. $K_P = kN$, for some positive constant $k$, which is the conversion factor of prey into predators. Also, a positive constant term added with $K_P$ to avoid a mathematical singularity when prey population becomes zero. From a biological point of view, for the case of severe scarcity of prey, predators can switch over other populations (alternative food), but its growth will be restricted due to non availability of its favorite food. Also, adding this type of positive constant introduces a maximum decrease rate, which stands for environment protection.

Thus we consider the following modified Leslie-Gower model within homogeneous environment with Holling functional response of type Ⅱ corresponds to the prey consumption rate of predators:

where all variables and parameters are positive, and they are defined in Table 1.

Many authors studied the Leslie-Gower type predator-prey model and its modified version. Aziz-Alaoui and Okiye ^[3] showed the global stability of the positive interior equilibrium point in a modified Leslie-Gower model. Nindjin et al. ^[4] analyzed the modified Leslie-Gower model with time delay. Gupta and Chandra ^[5] studied the modified Leslie-Gower model with Michaelis-Menten type prey harvesting and show rich complex-dynamical behaviours. Zhu and Wang ^[6] showed the global attractive behaviour of positive periodic solutions for the modified Leslie-Gower model, where predator-prey interaction follows Holling type Ⅱ scheme.

Now, we consider two different aspects to describe the predator-prey interactions:

1. Predators cooperate during hunting and as a result, success of capturing a prey increases with predator density.

2. Due to fear of predation risk, the birth rate of prey population reduces.

The role of mutual beneficial interactions (e.g. cooperation) among social animals is a rapidly growing research field in population dynamics and conservation biology. Cooperative hunting is the most widely distributed form of cooperative behaviour in animals. Cooperative hunting is widespread among vertebrates (e.g., lions ^[7], African wild dogs ^[8], chimpanzees ^[9], wolves ^[10]), birds ^[11], primates ^[12] etc.

Different animals adopt different strategies for making a group. For instance, wolves have the strategies to enclose the prey and catch it, for the acquisition of larger and faster prey ^[13]. Group hunting in lionesses generally involves a formation whereby some lionesses circle prey while others wait for prey to move towards them ^[7]. For wild chimpanzees, when a large prey is subdued, all the chimpanzees scream with excitement and trying to reach it ^[9]. Painted dog forms a rally, rousing one another and getting each other for group hunting ^[14]. Bailey et al. ^[15] discussed different strategies of cooperative hunting, mutual benefit for living in a social group and how predators form a cooperative group with active individuals.

In general, animals of same species compete with each other for resources, but because of some specific mechanism is at work, they help each other for mutual benefit and show cooperative behaviour. Two main constraints that lead to cooperative behaviour of the individuals are need for assistance in hunting and killing large prey, and the need for defence against attack by other predators ^[16,17]. There are several benefits of group hunting, such as the probability of capturing large prey increases ^[18], rate of hunting success increases with the number of adults ^[8], chasing distance decreases ^[8], locate food more quickly with increasing group size ^[19], helps to prevent the carcass from being stolen by other predators ^[20,21] etc. Packer and Ruttan ^[22] enlisted $61$ species of mammals, birds, vertebrates and invertebrates, their hunting strategy, which prey or types of prey to predate, single or multiple preys to predate at a time and percentage of their hunting success, from different sources.

Predators can impact the ecology and evolution of their prey directly by eating them, but also indirectly by influencing the behaviour and physiology of survivors ^[23]. This types of behavioural and physiological changes due to the fear of predation are more powerful and longer-lasting evolutionary consequence compared to direct predation ^[23,24]. For example, birds respond to the sounds of predators with anti-predator defences and flee from their nests at the first sign of danger ^[24]. Mule deer spends less time for foraging due to the predation risk of mountain lions ^[25]. Creel et al. ^[26] observed that reproductive physiology of elk changes due to the predation risk of wolves.

Many experimental studies show that manipulation of fear is strong enough to affects the population dynamics of ecological systems. Zanette et al. ^[27] manipulated fear of song sparrows during the entire breeding season by eliminating direct predation and observed that perception of predation risk alone reduces the number of offspring produced per year by 40%. Hua et al. ^[28] observed that perception of predation risk negatively affect blue bird's reproduction by manipulating the vocal cues of their predators. Suraci et al. ^[29] experimentally showed that manipulation of fear of large carnivore's vocalization in free living mesocarnivores reduces their foraging, which is beneficial for mesocarnivore's prey. This type of anti-predator behaviour although beneficial for increasing the probability of survival but driving long-term costs on reproduction that may affect population numbers ^[24]. However, very recently, Peers et al. ^[30] argued that predator cue experiments may not be realistically comparable to those of the prey's natural sensory environment, and the understanding of fear effects will continue to rely on untested assumptions in the wild. How prey population sustains potentially risky activities such as feeding or mating against anti-predator demands, has been a major issue in modern behavioural ecology ^[31,32,33]. In the presence of predators, prey change their behaviour due to predation risk. These behavioural responses to predation risk include reduce foraging time ^[34], increase vigilance ^[35], reduce movement ^[36], retreat to low-risk areas or refuge ^[37] and change in group size ^[38].

In mathematical modelling approach, many authors investigated the impact of hunting cooperation ^{[39,40,41,42,43]} and the impact of fear ^{[44,45,46,47,48]} in the predator-prey systems separately. In the context of hunting cooperation, Duarte et al. ^[39] investigated population dynamics of cooperative hunting in a McCann-Yodzis model by considering a fraction of predators cooperates during hunting. Berec ^[40] studied predators foraging facilitation or hunting cooperation in a Rosenzweig-MacArthur model and proposed a family of functional responses by considering attack rate and handling time of predators varies with predator density, in a Holling type Ⅱ functional response. Teixeira Alves and Hilker ^[41] showed that cooperative hunting in predators enhances the stability scenario of the coexistence steady state and it produces oscillatory dynamics in a predator-prey model, which is not possible in the absence of hunting cooperation. Recently, Pal et al. ^[43] investigated a discrete-time predator-prey model with hunting cooperation and observed that hunting cooperation has the potential to modify the well-known period-doubling route to chaos by reverse period-halving bifurcations and makes the system stable, but very high hunting cooperation can be detrimental for the system since it ultimately leads to extinction. On the other hand, recently, few researchers have investigated the impact of fear effect in a predator-prey system with the help of mathematical modeling. In this context, Wang et al. ^[44] first proposed a mathematical model by incorporating fear of the predator on prey, where the cost of fear reduces in the birth rate of prey. Their mathematical analysis suggests that high levels of fear can stabilize the predator-prey dynamics by excluding the existence of periodic solution. Wang and Zhou ^[45] investigated a predator-prey model with maturation delay for prey under the cost of fear and adaptive avoidance of predators. Their numerical findings showed that both strong adaptation of adult prey and the large cost of fear destabilize the predator-prey dynamics but higher population density stabilize the system. Panday et al. ^[46] studied a three-species food chain model, by considering the growth rate of prey and middle predator reduced due to the cost of fear of middle predator and top predator, respectively. Their numerical findings showed that fear has the potential to stabilize a chaotic system. Pal et al. ^[48] studied the impact of fear in a predator-prey model with Beddington-DeAngelis functional response and observed that the model exhibits multiple Hopf-bifurcations. Sasmal ^[47] studied ecological and epidemiological models with strong Allee effect and cost of fear in prey reproduction. Recently, Sha et al. ^[49] investigated an eco-epidemiological system with fear effects and observed that fear can produce backward bifurcation and chaotic oscillations in a simple eco-epidemic model.

In the present paper, we are interested to observe the dynamics of the such system where predators show cooperative behaviour during hunting and also create fear upon their prey ^[50,51,52]. For example, wolves cooperate during hunting, and when wolves are present, elk used anti-predator strategies and avoided areas frequented by the wolves ^[53]. Elk responds in the presence of wolves by altering patterns of agammaegation, habitat selection, vigilance, foraging, and sensitivity to environmental conditions ^{[38,54,55,56]}. Lioness shows cooperative behaviours during hunting ^[7]. Due to the fear of predation risk, zebras reached areas where the encounter with lionesses less frequent ^[51]. Thus by both killing and frightening, predators could have a dual impact upon their prey. The above observations motivate us to study the dynamics of predator-prey system including cooperative hunting by predators and a reduction of prey birth rate due to fear of predation risk. To the best of our knowledge, the effect of hunting cooperation and fear in the predator-prey system simultaneously has not yet been studied, and is the first attempt in this direction. The main focus of our present study is to explore the following ecological issues:

1. How does cooperative hunting impact (stabilize or destabilize) on the dynamics of a modified Leslie-Gower system?

2. Can fear effect stabilize the predator-prey system?

3. Can cooperative hunting and fear effect simultaneously promote species persistence and enhance ecosystem stability?

To incorporate the cooperative behaviour of predators in the system (1.1), we briefly outline Berec's approach for cooperative behaviours among predators ^[40]. Berec discussed cooperative phenomena in a Holling type Ⅱ functional response with two different approaches.

Firstly, when encounter rate of predators ($\rho$) increases with increasing predator density and handling time ($h$) is constant i.e.

then type Ⅱ functional response with this $\rho$ and $h$ is called encounter-driven functional response.

Secondly, if the encounter rate of predators $\rho$ is constant and handling time $h$ decreases with increasing predator density i.e.

then type Ⅱ functional response with this $\rho$ and $h$ is called handling-driven functional response.

In this paper, we confine ourselves with Berec's encounter-driven functional response. In particular, we choose $w = -1$ and so $\rho(P) = \rho_0(e+P) = \mu+\nu P$, where $\mu$ corresponds to $\rho_0e$ and $\nu$ corresponds to $\rho_0$. Also for the Leslie-Gower model (1.1), positive constant $k$ which is stated earlier as the conversion factor of prey into predators, affected with the predator density. Thus we consider $k = k_0(\mu+\nu P)$.

In Leslie-Gower or modified Leslie-Gower model, the carrying capacity of the predator population is variable and grows linearly with prey population density. However, in the present paper, we assume that the environment carrying capacity of the predator grows with both prey and predators abundances. In terms of foraging theory, resource supply rates likely increase due to cooperative foraging, such as reduced search costs and increased encounter rates ^[57]. So, due to cooperation, predators can consume more prey, which can increase the carrying capacity for the predator population. Zhang ^[58] theoretically demonstrated that cooperation or mutualism can promote the competitive ability of competitors through increasing their carrying capacities. Hamilton et al. ^[59] suggested a new mechanism for the biogeographic expansion of modern humans by a modified density-dependent model with varying intraspecific competition for finite resources. They showed that cooperative behaviors reduce intraspecific competition and increase population carrying capacity. They also suggested that cooperation and sociality in modern humans led to increased regional carrying capacities of humans.

Following the above observations, while formulating the following mathematical model, we use the formalism of increasing the cooperative behavior of predator increases its carrying capacity. Therefore, the model (1.1) with the Berec's encounter-driven functional response becomes:

It is to be noted that in the absence of hunting cooperation ($\nu = 0$), the above model (1.2) will be exactly same as the modified Leslie-Gower model (1.1).

Our second assumption is that birth rate of prey population is reduced due to fear induced by predators. To incorporate the fear phenomenon, we multiply the reproduction term i.e., birth rate ($b_0$) of prey population with a decreasing function of the predator population size, $\phi(\tau, P) = \frac{1}{(1 + \tau P)}$. Similar fear function was considered by many authors ^{[44,46,47,60]}. The parameter $\tau(>0)$ represents the level of fear that reduces the growth rate of prey. Expression of $\phi(\tau, P)$ is appropriate since it satisfies the following biological assumptions:

1. $\phi(0, P) = 1:$ if there is no anti-predator behaviours of prey, then the birth rate of prey remains unchanged.

2. $\phi(\tau, 0) = 1:$ if the predator population becomes zero, then no reduction in prey reproduction due to anti-predator behaviours.

3. $\lim_{\tau\to\infty} \phi(\tau, P) = 0:$ when anti-predator behaviour is very large and predator density $P > 0$, then prey reproduction decreases and ultimately becomes zero.

4. $\lim_{P\to\infty} \phi(\tau, P) = 0:$ when predator population is very large and strength of fear $\tau > 0$, then prey reproduction decreases and ultimately becomes zero, due to large anti-predator behaviours.

5. $\frac{\partial \phi(\tau, P)}{\partial \tau} < 0:$ the reproduction of prey decreases with the increasing of anti-predator behaviours.

6. $\frac{\partial \phi(\tau, P)}{\partial P} < 0:$ the reproduction of prey decreases with the increasing of predator populations.

Thus the model (1.2) with fear phenomena becomes

To reduce the number of parameters, we rescale the variables and parameters and the system (1.3) gets transformed to

This is the model we will analyze in the following. The new variables and parameters are given as follows: $x = \frac{a_0}{c_0}N$, $y = \frac{a_0}{hc_0^2}P$, $t = c_0T$, $b = \frac{b_0}{c_0}$, $d = \frac{d_0}{c_0}$, $\alpha = \frac{ \nu hc_0^2}{\mu a_0}$, $\beta = \frac{hc_0^2\tau}{a_0}$, $p = \frac{a_0}{\mu hc_0}$, $q = \frac{h}{\mu k_0}$, $r = \frac{m_0a_0}{\mu c_0k_0}$.

2.   Mathematical preliminaries

In this section, we present some basic results, such as positivity, boundedness and permanence of the solutions of system (1.4).

2.1.   Positivity and boundedness of the solutions

To prove positivity of the model (1.4), we follow the theorem given in H.R. Thieme's book ^[61].

From the standpoint of ecological meaning, we consider the following region $\mathbb{R}^2_+ = \{(x, y): x\geq0, y\geq0\}$. Now the function $F_1(x, y) = xf(x, y)$ and $F_2(x, y) = yg(x, y)$ of the system (1.4) are continuously differentiable and locally Lipschitz in $\mathbb{R}^2_+ = \{(x, y): x\geq0, y\geq0\}$. Therefore, Theorem A.4, page 423 in H. R. Thieme's book ^[61] implies that the solutions of the initial value problem with non-negative initial conditions exist on the interval $[0, S)$ and unique, where $S$ is a sufficiently large number.

In theoretical ecology, boundedness of the system (1.4) implies that the system is well behaved. Boundedness of the solutions entails that none of the interacting species grow abruptly or exponentially for a long-time interval. The number/abundance of each species is bounded due to limited resource. To prove the boundedness of the system (1.4), we state the following lemma whose proof is given in ^[62].

Lemma 2.1. Assume that $\xi, \eta > 0$ with $u(0) > 0$. Then for the differential inequation $\frac{du}{dt}\leq u(t)(\xi-\eta u(t))$,

and also for the differential inequation $\frac{du}{dt}\geq u(t)(\xi-\eta u(t))$,

Proposition 2.2. Assume that $b > d$ and $q > \alpha(b-d+\epsilon_1)$, where $\epsilon_1 > 0$ is a arbitrary small number. Then all the solutions $(x(t), y(t))$ of the system (1.4) with initial conditions $x(0) > 0, \; y(0) > 0$ are bounded.

Proof. From the first equation of the system (1.4)

Assume that $b > d$. Now applying the Lemma 2.1, we have

Thus for arbitrary small $\epsilon_1 > 0$, there exists a positive real number $T_1 > 0$ such that

Also, the second equation of the system (1.4) can be written as

Thus, for all $t > T_1$, we have

Assume that $q > \alpha (b-d+\epsilon_1)$. Then again applying the Lemma 2.1, we have

Thus, for arbitrary small $\epsilon_2 > 0$, there exist a positive real number $T_2 > T_1$ such that

Therefore, the solutions of the system (1.4) are bounded.

2.2.   Permanence:

Geometrically, uniform permanence means the existence of a region in the phase plane at a non-zero distance from the boundary in which population species enter and must lie ultimately that ensures the long time survival of species in biological sense. Analytically the system (1.4) is said to be permanent if there exist positive constants $m$ and $M$ such that each positive solution $(x(t, x_0, y_0), y(t, x_0, y_0))$ of the system (1.4) with initial condition $(x_0, y_0)\in int (\mathbb{R}^2_+)$ satisfies

Proposition 2.3. The system (1.4) with initial conditions $x(0) > 0$, $y(0) > 0$ is permanent if $\frac{b}{1 +\beta M_2}-d-\frac{ (1+\alpha M_2)M_2}{p} > 0.$

Proof. From the prey equation of the system (1.4), we have

Assume that $w_0 > 0$. Then by applying the Lemma 2.1, we have

Also, from the predator equation of system (1.4), we have

Again, by applying the Lemma 2.1, we get the following result

Choose $m = \min\{w_0, \frac{r}{q}\}$ and from (2.1) and (2.2), choose $M = \max\{b-d, \frac{r+b-d+\epsilon_1}{q-\alpha (b-d+\epsilon_1)}\}$. Then the condition of permanence of the system (1.4) follows.

3.   Equilibria & stability analysis

The equilibria of the system (1.4) occur at the intersection of prey nullclines and predator nullclines governed by the following equations:

The non-trivial prey nullcline $f(x, y) = 0$ represents a cubic equation of $x$ and $y$ and non-trivial predator nullcline $g(x, y) = 0$ represents a hyperbolic curve of the form

The predator-prey model (1.4) has following trivial and boundary equilibria:

(ⅰ) Trivial equilibrium point $E_0 = (0, 0)$.

(ⅱ) Prey only equilibrium point $E_1 = (b-d, 0)$.

(ⅲ) Predator only equilibrium point $E_2 = (0, \frac{r}{q})$.

Now, the positive interior equilibrium points of the system (1.4) are the solution of the following system of equations:

Let $(\bar x, \bar y)$ be a positive solution of the system of equations (3.1). Then we have

and $\bar y$ is a positive root of the equation

where

Here $\sigma_0$ and $\sigma_1$ are positive, so the above equation (3.3) can have atmost three positive roots depending on the sign of $\sigma_2$, $\sigma_3$ and $\sigma_4$.

Now using the linearized stability analysis, we can obtain the local stability criteria of the above mentioned equilibrium points. The community matrix at $E_0 = (0, 0)$, $E_1 = (b-d, 0)$ and $E_2 = (0, \frac{r}{q})$ are respectively given by

Thus the equilibrium points $E_0 = (0, 0)$ and $E_1 = (b-d, 0)$ both are unstable. The prey free equilibrium point $E_2 = (0, \frac{r}{q})$ is locally asymptotically stable if the cooperation strength $\alpha$ greater than the threshold value $\frac{bpq^3}{r^2(\beta r+q)}-\frac{dpq^2}{r^2}-\frac{q}{r}$, or equivalently, the fear strength $\beta$ greater than the threshold value $\frac{bpq^3}{r^2(\alpha r+q)+dpq^2r}-\frac{q}{r}$. Also the community matrix at $\bar E = (\bar x, \bar y)$ is given by

The trace and determinant value of the above community matrix are given by

Both of $T$ and $D$ can be positive or negative. So the stability of interior equilibrium point depends on both the trace and determinant value of the community matrix. Therefore, by using Routh-Hurwitz criterion, the interior equilibrium $\bar E(\bar x, \bar y)$ is locally asymptotically stable if $T(\bar x, \bar y) < 0$ and $D(\bar x, \bar y) > 0$.

4.   Hopf-bifurcation and existence of multiple limit cycles

Nonlinear mathematical models of interacting populations show rich and complex dynamical behaviours even when system complexity is low (two or three species). Oscillating behaviour is a most frequent dynamical behaviours in population dynamics. Oscillating behaviour or the existence of limit cycle lead to the Hopf-bifurcation of the system. Hopf-bifurcation is defined as the appearance or disappearance of a periodic orbit through a local change in the stability properties of an equilibrium point. Here we explore the possibility of occurrence of Hopf-bifurcation and the direction of Hopf-bifurcation around the interior equilibrium point $\tilde E(\tilde x, \tilde y)$ with respect to bifurcating parameter $\alpha$.

4.1.   Hopf-bifurcation

The following theorem gives the necessary and sufficient condition for Hopf-bifurcation of the system (1.4).

Theorem 4.1. The necessary and sufficient conditions for the system (1.4) undergoes Hopf-bifurcation at $\alpha = \alpha_h$ around the interior equilibrium point $\tilde E$ are that $D(\tilde x(\alpha), \tilde y(\alpha)) > 0$ and $\frac{d}{d\alpha}T(\tilde x(\alpha), \tilde y(\alpha))\neq0$ at $\alpha = \alpha_h$.

Proof. Consider the community matrix (3.4) evaluated at $\tilde E$. For Hopf-bifurcation, the matrix $J_{\tilde E}$ must have a pair of purely imaginary eigenvalues. At a critical value $\alpha = \alpha_h$, $T(\tilde x(\alpha), \tilde y(\alpha)) = 0$. So the characteristic equation of the community matrix becomes

The above equation must have a pair of purely imaginary roots $\lambda_{1, 2} = \pm i\theta_0$, where $\theta_0 = \sqrt{D(\tilde x(\alpha), \tilde y(\alpha))|_{\alpha = \alpha_h}}$, if $D(\tilde x(\alpha), \tilde y(\alpha)) > 0$ at $\alpha = \alpha_h$.

Now we check the transversality condition which confirms that the eigenvalues crosses the imaginary axis of the complex plane with non-zero speed. For this, let at any neighbouring point $\alpha$ of $\alpha_h$, the eigenvalues of the community matrix are $\lambda_{1, 2} = \chi(\alpha)\pm i\theta(\alpha)$, where $\chi(\alpha) = \frac{T(\tilde x(\alpha), \tilde y(\alpha))}{2}$ and $\theta(\alpha) = \sqrt{D(\tilde x(\alpha), \tilde y(\alpha))-\frac{T^2(\tilde x(\alpha), \tilde y(\alpha))}{4}}$. Now,

Thus, the transversality condition satisfied if $\frac{d}{d\alpha}T(\tilde x(\alpha), \tilde y(\alpha))\neq0$ at $\alpha = \alpha_h$.

Therefore, the system (1.4) undergoes Hopf-bifurcation at $\alpha = \alpha_h$.

4.2.   Direction of Hopf-bifurcation

Here, we discuss the direction and stability properties of the bifurcating periodic solutions originating from the coexisting equilibrium point $\tilde E(\tilde x, \tilde y)$ via Hopf-bifurcation. To investigate the stability and direction of Hopf-bifurcation, we calculate the $1^{st}$ Lyapunov coefficient by using the results given in ^[63].

First, we transform the equilibrium point $\tilde E(\tilde x, \tilde y)$ of the system (1.4) into the origin by letting $z_1 = x-\tilde x, z_2 = y-\tilde y$. Then the system (1.4) becomes

Expanding Taylor's series of the above system at $(z_1, z_2) = (0, 0)$ up to terms of order $3$ produces the following system:

where

By neglecting the higher order terms of degree $4$ and above, the system (4.1) can be written as

where

The eigenvector $\tilde v$ of community matrix $J_{\tilde E}$ corresponding to the eigenvalue $i\theta_0$ at $\alpha = \alpha_h$ is $\tilde v = (p_{01}, i\theta_0-p_{10})^T$. Now define

Let $Z = QU$ or $U = Q^{-1}Z$, where $U = (u, v)^T.$ Under this transformation, system (4.2) becomes

This can be written as,

where $S^1$ and $S^2$ are non-linear in $u$ and $v$ given by,

with

Now, we calculate the $1^{st}$ Lyapunov coefficient, based on the normal form (4.3), which determines the stability and direction of periodic solution. The $1^{st}$ Lyapunov coefficient is

where all the partial derivatives are calculated at the bifurcation point, i.e., $(u, v; \alpha) = (0, 0;\alpha_h)$.

By applying the result given in ^[63] (see equation (20.2.14), page 385), Hopf-bifurcation is supercritical if $l_1 < 0$ and it is subcritical if $l_1 > 0$. It is to be noted that when $l_1 = 0$, the system (1.4) exhibits generalized Hopf-bifurcation (or, Bautin bifurcation) at which the interior equilibrium has a pair of purely imaginary eigenvalues. The generalized Hopf-bifurcation point separates branches of subcritical and supercritical Hopf-bifurcations in the parameter plane.

5.   Bogdanov-Takens bifurcation

In this section, we discuss Bogdanov-Takens (BT) bifurcation of the system (1.4) in a small neighbourhood of the interior equilibrium point $\tilde E(\tilde x, \tilde y)$. BT bifurcation is a bifurcation of an equilibrium point in a two-parameter family of autonomous system and occurs for a system when the community matrix evaluated at the equilibrium has a zero eigenvalue of algebraic multiplicity $2$. In such a situation the community matrix is similar to the Jordan block of the form $\begin{pmatrix}0 & 1\\0 & 0\end{pmatrix}$. We follow the techniques and steps of Kuznetsov's book ^[64] to derive a normal form by using a series of transformations. Now we consider the system (1.4) and choose $\alpha$ and $p$ as an arbitrary family of parameters and keeping all the other parameters are fixed. Without loss of generality, we define $\tilde x(\alpha, p)\equiv \tilde x$ and $\tilde y(\alpha, p)\equiv \tilde y$. It is to be noted that BT bifurcation may occur for other parameters also. Since the all the parameters are positive and $y$ is non-negative, so the functions $F_1(x, y) = xf(x, y)$ and $F_2(x, y) = yg(x, y)$ in the system of equations (1.4) are smooth function of $x$ and $y$. Assume that at $\alpha = \alpha_{bt}$ and $p = p_{bt}$, $T(\tilde x, \tilde y) = 0, D(\tilde x, \tilde y) = 0$. Then the community matrix in (3.4) evaluated at $\tilde E(\tilde x, \tilde y)$ has a zero eigenvalue of algebraic multiplicity $2$. Now we consider the following parametric region

Mathematically, the surface represented by $\bar\Omega$ is called a BT surface. In the following, we reduce the system (1.4) in the canonical form of a BT bifurcation by employing a series of $C^{\infty}$ change of coordinates ^[64] in a small neighborhood of $(0, 0)$. For these, first we use the transformation $x = u_1+\tilde x$ and $y = u_2+\tilde y$, so that the equilibrium point of the system (1.4) shifted to the origin. Then the model system (1.4) can be written as

Taking the small perturbation of the model (5.1) at $\alpha = \alpha_{bt}$ and $p = p_{bt}$, we get

where $\delta = (\delta_1, \delta_2)$ is a parameter vector vary in a small neighborhood of $(0, 0)$. Expanding the Taylor's series of the above system at $(u_1, u_2) = (0, 0)$ up to the terms of order $2$, produces the following system:

where

Here note that $a_{00}(0) = 0 = b_{00}(0)$ and

Now we introduce affine transformation

which reduce the system (5.2) to

where

Next, under the following $C^{\infty}$ change of coordinates in a small neighbourhood of $(0, 0)$

the system (5.3) becomes

where

It is easily observed that $e_{00}(0) = 0 = e_{10}(0) = e_{01}(0)$, since $d_{00}(0) = 0 = d_{10}(0) = d_{01}(0)$ and

We now use the following parameter-dependent shift of coordinates in the $x_1$-direction to annihilate the $x_2$ term on the RHS of the second equation of (5.4)

Using the above transformation system (5.4) reduces to

where

Assume that

Using standard arguments based on the Implicit Function Theorem gives local existence of a smooth function

annihilating the term proportional to $y_2$ in the $2$nd equation of (5.5). Here we note that $\phi(0) = 0$, thus $f_{00}(0) = 0 = f_{10}(0) = f_{01}(0)$, $f_{11}(0) = e_{11}(0)$ and $f_{20}(0) = e_{20}(0)$. Now we introduce a new time scale, defined by $dt = (1+\Lambda z_1)d\theta$, where $\Lambda = \Lambda(\delta)$ is a smooth function to be defined later. Also the direction of time is preserved near the origin for small $||\delta||$. This transformation reduces the system (5.5) to

Now assuming $z_1 = y_1, \quad z_2 = y_2(1+\Lambda y_1)$, the system (5.6) becomes

To remove the $z_2^2$ term we specify the time reparametrization as $\Lambda(\delta) = -f_{02}(\delta)$. Consequently, we have

where

We now introduce a new time scale and denote it by $t$ again as follows $t = |\frac{N(\delta)}{M(\delta)}|\theta$. Since $N(0) = f_{11}(0) = e_{11}(0)\neq0$ due to (BT.1), the time scaling above will be well defined if we further assume

Simultaneously, we introduce new variables $\chi = \frac{N^2(\delta)}{M(\delta)}z_1$, $\psi = sign\left(\frac{N(\delta)}{M(\delta)}\right)\frac{N^3(\delta)}{M^2(\delta)}z_2$, which takes the system (5.7) into the following form

where

In order to define an invertible smooth change of parameters near the origin, we assume that $rank(\frac{\partial\beta}{\partial\delta})\Big|_{\delta = 0} = 2,$ i.e., at $\delta = (\delta_1, \delta_2) = (0, 0)$

The system (5.8) is called normal form of Bogdanov-Takens bifurcation. Thus under the assumptions of non-degeneracy conditions (BT.1), (BT.2) and (BT.3), the system (1.4) can be reduced to the system (5.8). Also by using Lemma $8.8$, page $324$ in ^[64], the system (5.8) is locally topologically equivalent near the origin to the system

Thus the system (1.4) undergoes BT bifurcation by choosing $\alpha$ and $p$ as a bifurcation parameter when the parameters $\delta_1$ and $\delta_2$ vary in a small neighbourhood of the origin.

6.   Numerical simulations

To justify our analytical findings and to depict the effect of cooperation and fear on the dynamics of the system (1.4), we perform extensive numerical simulations. In the present paper, all the numerical findings are performed in Matlab by using $ode45$ method and Matcont package with suitably chosen small time steps. In order to understand the role of cooperation and fear on the stability and bifurcation of the coexistence equilibrium point, first we consider the system (1.4) without hunting cooperation ($\alpha = 0$) and without fear effect ($\beta = 0$). We choose the parameter values as

with initial condition $(x_0, y_0) = (4, 1)$. First, we observe that the system (1.4) for the above set of parameter choice shows stable dynamics (Figure 1).

6.1.   Multiple Hopf-bifurcations and existence of limit cycles

Now, our aim is to investigate the dynamics of the system with hunting cooperation ($\alpha\neq0$) but without fear effect ($\beta = 0$). We observe that an increase in the strength of hunting cooperation makes the system unstable, and both the predator and prey populations show periodic oscillations (Figure 2). The system (1.4) loses its stability via Hopf-bifurcation at $\alpha = 0.512218 \; (\alpha_l)$. We also observe that further increase of the strength of $\alpha$, the system (1.4) undergoes another Hopf-bifurcation at $\alpha = 6.510087 \; (\alpha_u)$ and becomes stable. However, for very large values of $\alpha\; (\alpha > \alpha_e = 20.5)$ prey population goes to extinction due to overexploitation by the predator population (Figure 3). Therefore, the system (1.4) undergoes multiple Hopf-bifurcations at $\alpha_l = 0.512218$ and $\alpha_u = 6.510087$, and shows periodic oscillation when strength of hunting cooperation is in the range $\alpha_l < \alpha < \alpha_u$. We notice that there is a hump-shaped relationship between predator density at the equilibrium and the parameter $\alpha$ (Figure 3). The predator density increases when the strength of hunting cooperation low, but declines when the strength of hunting cooperation large due to overexploitation of prey by the predators. Such type of relationship between predator density at the equilibrium $(\tilde y)$ and the cooperation strength ($\alpha$) was also observed by Teixeira Alves and Hilker ^[41] in a Lotka-Volterra predator-prey model with hunting cooperation among predators. The biological implications for predator population increases for small values of $\alpha$ because of predators better forage on prey, but for large values of $\alpha$ predator population decreases as decrease in prey density. Also, due to overexploitation of predators, extinction of prey population occurs at the threshold value $\alpha = 20.5$ (Figure 4).

6.2.   Supercritical Hopf-bifurcation

A supercritical Hopf-bifurcation means that the stable periodic solution overlaps an unstable equilibrium point. In the Subsection 4.2, we have obtained the sufficient condition for supercritical Hopf-bifurcation of the system (1.4) theoretically. Now for the set of parameter values (6.1), we calculate the $1^{st}$ Lyapunov coefficient corresponding to the threshold value $\alpha = \alpha_l$. The $1^{st}$ Lyapunov coefficient is $L_l = -0.072533$, which confirms that the limit cycle obtained in our system is stable i.e. the trajectories starting from different initial conditions converge to the same limit cycle. We draw the limit cycle when $\alpha = 1$ and observe that the stable interior equilibrium point becomes unstable and a stable limit cycle is born (Figure 5). Therefore, the system undergoes supercritical Hopf-bifurcation.

6.3.   Subcritical Hopf-bifurcation

A subcritical bifurcation means that the unstable periodic solution overlaps a stable equilibrium point. Now, we calculate the first Lyapunov coefficient for the system (1.4) at the threshold value $\alpha = \alpha_u$. The first Lyapunov coefficient is $L_u = 0.069969$, which confirms that there exists a subcritical Hopf-bifurcation. We draw the solution trajectories when $\alpha = 6.7$ and observe that the unstable equilibrium point becomes stable and there is an unstable limit cycle present around the stable equilibrium point $(0.0635, 1.8501)$ (Figure 6). We also observe that the system shows a stable limit cycle outside of the unstable limit cycle. So, the system (1.4) also exhibits subcritical Hopf-bifurcation and hunting cooperation can change the stability property of the limit cycle. The subcritical Hopf- bifurcation leads to a significant domain of bi-stability, where stable steady state and time-periodic state coexist. Here, the prey and predator population tend to the interior equilibrium point, if the initial populations start inside the unstable limit cycle and if the initial populations start outside the unstable limit cycle, then the prey and predator population converge to the stable limit cycle. So, in the absence of fear effect $(\beta = 0)$ the system (1.4) shows stable dynamics and limit cycle (stable/unstable) oscillations depending on the strength of hunting cooperation $\alpha$.

6.4.   Fear induced stabilization

Now, we are going to investigate the impact of fear effect on the dynamics of the system (1.4). We observe that for a fixed value of $\alpha$, if we increase the strength of fear parameter ($\beta$), then the system (1.4) remains stable or becomes stable. Here, we choose two different values of $\alpha$, where in the absence of fear effect $(\beta = 0)$ the system (1.4) shows stable focus $(\alpha = 0.1)$ and periodic oscillations $(\alpha = 1)$. Now, if we incorporate fear phenomenon into the system with fear strength $\beta = 0.1$, then stable system remains stable and oscillating system becomes stable (Figure 7). We also draw two parametric bifurcation diagram in $\alpha-\beta$ parametric space and observe that fear factor ($\beta$) has a stabilizing effect on the dynamics of the system (Figure 8). It is also observed that the system (1.4) exhibits GH bifurcation (generalized Hopf bifurcation) at the threshold value $\alpha = 4.7255712(\equiv \alpha_{gh})$ and $\beta = 0.062264(\equiv \beta_{gh})$, where the $1^{st}$ Lyapunov coefficient becomes zero. Further, for large values of hunting cooperation ($\alpha$) and fear factor ($\beta$), the prey population goes to extinction from the system whereas predator population converges to a fixed value (Figure 9).

6.5.   Bi-stability

Previously, we have observed the bi-stability behaviour between coexistence equilibrium and periodic oscillations. Here, we also observe bi-stability behaviour between prey-free equilibrium point and coexistence equilibrium point. First, we draw the bifurcation diagram of the system (1.4) with respect to the parameter $\alpha$ (Figure 10), where $p = 0.3$, $\beta = 0.1$ and other parameters are same as in (6.1). Here, we draw black curve for stable branch and magenta curve for unstable branch of the interior equilibrium, respectively. We also draw green line for stable branch and red line for unstable branch of the prey-free equilibrium, respectively. We observe that there exists an interval $(0.1536 < \alpha < 0.2212)$, where both the interior equilibrium point and the prey-free equilibrium point are stable (Figure 10). Further, by considering parameters from the bistable region, we draw the basin of attraction of the prey-free and interior attractors. We observe that for initial conditions in the black region, the trajectories converge to the interior equilibrium point $(1.1371, 2.7662)$ and for initial conditions in green region the trajectories converge to the prey-free equilibrium point $(0, 1)$ (Figure 11(a)). For better visualization, we also draw phase portrait of the system (1.4) by considering different initial conditions from green and black regions (Figure 11(b)).

6.6.   Bogdanov-Takens bifurcation

The Bogdanov-Takens bifurcation (BT bifurcation) is a two-parametric bifurcation of an equilibrium point and at the bifurcation point, the community matrix evaluated at the interior equilibrium point has a zero eigenvalue of algebraic multiplicity two. In the neighbourhood of the bifurcation point, the system has two equilibrium points (one saddle and other non-saddle), which will collide and disappear via a saddle-node bifurcation. The non-saddle equilibrium undergoes Hopf-bifurcation, which degenerates into an orbit homoclinic to the saddle and disappears via a saddle homoclinic bifurcation. We draw two parametric bifurcation diagram in $\alpha-p$ parametric space, when $\beta = 0.1$ and other parameters are same as in (6.1). We observe that the system undergoes BT bifurcation at $\alpha = 0.8993835(\equiv \alpha_{bt})$ and $p = 0.4800637(\equiv p_{bt})$ (Figure 12). It is also observed that when both the eigenvalues are zero at $(\alpha_{bt}, p_{bt})$, there exists interior equilibrium point $(0.0878, 1.1810)$, which is a cusp of codimension 2 (Figure 13).

6.7.   Oscillation to stable

Again, we consider the parameter values as

so that the system (1.4) without hunting cooperation ($\alpha = 0$) and without fear effect ($\beta = 0$) shows periodic oscillations (Figure 14). First, we draw the bifurcation diagram of the system (1.4) w.r.t. cooperation parameter $\alpha$ (with $\beta = 0$) and observe that hunting cooperation alone can stabilize the oscillating system via Hopf-bifurcation (Figure 15). Again, we draw the bifurcation diagram of the system (1.4) w.r.t. fear parameter $\beta$ ($\alpha = 0$) and observe that fear phenomenon alone can also stabilize the oscillating system via Hopf-bifurcation (Figure 16). It is also observed that fear factor ($\beta$) has more stabilizing effect compare to hunting cooperation ($\alpha$).

7.   Discussion

In ecology and evolutionary biology, the mechanism of predator-prey interaction is a central topic. Both prey and predators use different strategies to maximize their biomass. To increase the success rate of catching a prey, many predators cooperate during hunting. On the other hand, due to fear of predation, prey populations show different types of anti-predator behaviours as their survival strategies. Prey population show anti-predator defences such as changes its habitat zone, forage less, increase vigilance, adjust their reproductive strategies etc. Such type of anti-predator behaviour, reduces the birth rate of prey populations.

In the present paper, we have considered a modified Leslie-Gower model with hunting cooperation and fear. Our main assumptions are: (ⅰ) predators can cooperate during hunting which increases their effectual consumption and therefore increases the carrying capacity of the predator population; (ⅱ) fear of predation can reduce the birth rate of prey population. We investigated the model step by step by considering different strengths of hunting cooperation and fear factor. To study the impact of hunting cooperation and fear in the predator-prey system, we considered two situations where the ecosystem showed stable dynamics and unstable (limit cycle oscillations) dynamics. Then we gradually increased the strength of hunting cooperation and the cost of fear and explored the impact of cooperation and fear in the ecological system. We observed that in the absence of fear if we increase the strength of hunting cooperation, then the system becomes unstable around the interior equilibrium and undergoes multiple Hopf-bifurcations, where the first Hopf-bifurcation is supercritical and the second one is a subcritical. So, the coexisting equilibrium point may lose stability via supercritical Hopf-bifurcation with an intermediate value of hunting cooperation and regain stability via subcritical Hopf-bifurcation if the strength of hunting cooperation is large. We have also observed the existence of multiple limit cycles and bi-stability phenomena via subcritical Hopf-bifurcation. Here, we found that the solution trajectories tend to a coexisting equilibrium point or oscillate periodically depending on the initial population size. We also investigated the impact of hunting cooperation and fear on the dynamics of the system when the system shows limit cycle oscillations around the interior equilibrium. We observed that if we increase the strength of hunting cooperation or cost of fear, then the population oscillations will be replaced by stable focus. Further, we have demonstrated the impact of the fear phenomenon that can have in the predator-prey interactions, when predators show cooperative behaviour during hunting. We have observed that fear factor has a stabilizing effect and with fear effect, the stable system remains stable and the oscillating system becomes stable by excluding the existence of periodic solutions. We also found that fear factor has more stabilizing effect compared to hunting cooperation and makes the system more robust. It is interesting to note that with increase in hunting cooperation and fear factor, the density of both prey and predator populations eventually will decrease. From the two parametric bifurcation diagram, we observed that for any cooperation strength, if the fear parameter is increased, then the system becomes stable. However, for low level of fear, if the cooperation level is increased then the system switches it's dynamics from stable to limit cycle oscillation and from limit cycle oscillation to stable focus. It is also to be noted that when the strengths of hunting cooperation and fear factor are very high, prey population goes to extinction due to overexploitation by the predators.

Cooperative hunting is a form of foraging facilitation, whose special feature is that it generates strong Allee effect in predators ^[40,41] in Lotka-Volterra model. Recently, in another manuscript, we have modified the Lotka-Volterra model by assuming that cooperation among predators induces fear in prey population ^[65]. We have observed strong demographic Allee phenomenon among predators, backward bifurcation and many complex dynamical behaviors. In the present paper, since we have considered a modified Leslie-Gower model (logistic-type growth for predators), Allee phenomenon is not possible in our model, although the predators show cooperative behaviour during hunting. We like to point out that in Lotka-Volterra model, hunting cooperation among predators enhances the growth rate of predator whereas in the modified Leslie-Gower model, predator hunting cooperation has a positive effect on the carrying capacity of predator. Allee phenomenon has ecological significance on the perspective of the persistence of the species. In Lotka-Volterra model with predator hunting cooperation Allee phenomenon is common ^[40,41,65] whereas for modified Leslie-Gower model this is not obvious. Thus from a predator persistence point of view modified Leslie-Gower model is more robust compared to Lotka-Volterra model. The system with hunting cooperation and fear phenomena exhibits rich dynamical behaviours including both supercritical and subcritical Hopf-bifurcation, bi-stability between steady state and limit cycle, bi-stability between prey-free equilibrium point and coexisting equilibrium point, generalized Hopf-bifurcation, Bogdanov-Takens bifurcation. So, hunting cooperation and fear effect are important factors and play vital roles in determining the long-term population dynamics of ecosystems stability.

Acknowledgments

Research work of Saheb Pal is supported by the Junior Research Fellowship from the UGC, Government of India. The authors are thankful to the three anonymous reviewers and the Editor for their valuable comments and suggestions, which helped us to improve the paper.

Conflict of interest

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