1. Introduction
Preypredator interaction has immense importance in the existence of life. In predatorprey systems, the prey population shows a variety of defense to escape from predation [
1,
2]. Among them, the refuge is most common in literature [
3,
4], and it has become a significant issue in theoretical ecology [
5]. One famous proverb about the FoxRabbit chase says: “The fox is only running for his dinner, but the rabbit is running for his life [
6]”. If the rabbit can successfully save its life, it’s called the rabbit refuge of the fox. So, the refuge can be defined as any strategy that decreases the rate of predation [
7].
Refuge can be done in two ways– (i) constant numerical refuge, where an definite number of individuals in the population refuge predation, and (ii) proportional refuge, where a fixed proportion of individuals in the population refuge predation. In both cases, the whole prey population can’t attain the facility of refuge, and thus, a particular number or fraction of the population is always exposed to the predator. If not, the predator population will no longer be able to persist, and the system will be destabilized.
If we consider a situation where a constant number of safe refuge places is available and thus only a definite number of individuals in the population can enjoy the facility of refuge. In such condition, two habitats are available for the prey population– either lives in the refuge or roam around in the open habitat where they are available to predator [
22]. In [
23], authors considered a preypredator model with constant numerical prey refuge and showed that, there are two major factors that control this type of refuge in the population: (i) population density of prey (
x) and available refuge (
m). When
$x\le m$, competition among prey for refuge is negligible since prey fitness depends only on predation [
22]. In such condition cooperation among individuals of the prey population must be seen. But this is not the case for all time because this type of refuge can include a definite number of individuals and thus they cannot move freely between habitats when population size exceeds refuge number (
$x>m$) and obviously intraspecific competition to take refuge will increase proportionally as density increases. The conclusion regarding this type of refuge is that it has a strong stabilizing effect [
24] it can prevent the extinction of prey, but when the population is too large to stay all in the refuge then coexistence with predators is expected as the predator will get continuous resources [
22]. So, in that case, mortality rate of the prey increases when population density goes beyond the refuge capacity and this will cause a negative feedback stabilizing effect to the system [
25].
Another situation is where a proportion of individuals in the population refuge the predation. In such condition, there are also two habitats available for the prey population– either lives in the refuge or roam around in the open habitat where they are available to predator [
22]. But in that case a continuous movement of prey individuals occur between open habitat and refuge [
26]. In [
27], Kar proposed the LotkaVolterra type preypredator system including proportional prey refuge where,
m is denoted as the fractional refuge of prey for predator and it ranges from
$0<m<1$. Now, if we consider the whole prey population as
x then
$mx$ portion of the prey is in refuge habitat, and
$(1m)x$ of the prey is in open habitat and thus available to the predator. In case of proportional refuge, the whole population never get the opportunity to attain refuge and thus, the process of finding refuge place is either cooperation or competition based. This type of refuge can reduce mortality rate of prey by reducing predation pressure and can effect positively the growth of prey and negatively that of predator. Thus, in this way refuge can act as a stabilizing effect in preypredator dynamics [
32]. On other extreme refuge may destabilize preypredator dynamics as the predator population become extinct by suffering from lower predation success due to large refuge size [
33].
Mathematically preypredator interaction was first described by Alfred James Lotka in 1925 [
8] and Vito Volterra in 1927 [
9]. Georgy Frantsevich Gause in 1934 [
10] first found refuge in his experiment with protozoan predator
Didinium feeding on
Paramecium. One major drawback of the LotkaVolterra model was that the predator never becomes saturated. Crawford Stanley Holling solved this problem [
11,
12,
13] by suggesting three kinds of functional responses (Holling type I, II, and III) for different species of predator. Rosenzweig and MacArthur in 1963 [
14] first replaced the linear functional response of Lotka–Volterra model with the Holling typeII response function and first documented that preypredator coexistence is not limited to a stable equilibrium; a limit cycle appears when the stable equilibrium undergoes the Hopf bifurcation [
15]. Hassel in 1978 [
20] showed that adding a large number of refuges to a model in which the absence of a refuge exhibits divergent oscillation replaces the oscillatory behavior with a stable equilibrium. Since then, various modeling approaches for the preypredator system, including refuge, have been done by researchers [
16,
17,
18,
19,
27,
28].
In the previous articles ([
28,
29,
31] and the references therein), authors have considered competition among predator populations. In our knowledge, this article is the first article to consider the impact of cooperation and intraspecific competition within the prey population for attaining refuge. The key question of this study is– “how does the preypredator system dynamics being impacted when the decision of refuge is made through cooperation and same for the intraspecific competition?” To get an answer to this problem, we have developed two mathematical models. The first model includes the cooperation among prey populations, whereas the second model contains the intraspecific competition. The dynamics of the models have been analyzed with qualitative theories.
This paper has been carried out sequentially in the latter sections as follows.
Section 2 deals with the mathematical model formation with prey refuge and cooperation and local stability and Hopf bifurcation also shown. Another mathematical model with prey refuge and intraspecific competition is formed and the existence of Hopf bifurcation and stability analysis is done in
Section 3. Global stability of coexisting equilibrium are analyzed in
Section 4.
Section 5 deals with findings of direction of Hopf bifurcation. Numerical simulation results are delineated in
Section 6. The final discussion and biological interpretation of the results of the present study are included in the concluding (
Section 7).
2. Mathematical Model for Prey Refuge and Cooperation
Here, we propose the LotkaVolterra type preypredator system including proportional prey refuge and cooperation making the following assumptions:
There are two populations in the model namely $x\left(t\right)$ is the prey population size, $y\left(t\right)$ is the population size of the predator at any time t. r is the intrinsic growth rate, and k is the carrying capacity of the prey, $\gamma $ is the death rate of the predator in the absence of prey, $\beta $ is the predation rate, and c is the conversion factor denoting the efficiency of predators for each captured prey.
Also, $\frac{\beta}{a}$ is the maximum number of prey that the predator per unit time can eat; $\frac{1}{a}$ being the halfsaturation constant.
The term $\frac{\beta x}{1+ax}$ denotes the predator’s response to prey. This type of predator response function is known as Holling typeII functional response.
Let m be the fractional refuge of prey for predator, $0<m<1$ and constant. So, $mx\left(t\right)$ portion of the prey is in refuge habitat and $(1m)x\left(t\right)$ of the prey is in open habitat and thus available to the predator, and $\theta $ is the cooperation coefficient among prey population trying to escape predation pressure through refuge.
Using the assumptions discussed above, we get the following mathematical model:
with the initial biological condition
It is assumed that all the system parameters are positive constants. The schematic diagram of the model (
1) is shown in
Figure 1.
In the next subsection, we determine all the possible equilibrium points and existence conditions for their biological feasibility of the model (
1).
2.1. Boundedness of Solutions
Biologically, all population of the model should be finite. For this boundedness analysis is significant. We have derived the following theorem for boundedness of the solutions of the model (
1).
Theorem 1. All solutions of the model (1) with the given initial conditions (2) are nonnegative and bounded for all time $t>0$.
Proof. We rewrite the first equation of (
1),
where,
$h(x,y)=r\left(\right)open="["\; close="]">1\frac{x}{k}$. Therefore
Since $x\left(0\right)>0$, we can conclude that $x\left(t\right)$ is always positive. Similarly, the nonnegativity of $y\left(t\right)$ can be established.
For boundedness analysis, we define
$W=x+y$. Now, adding the two equations of (
1), we can obtain
As by definition,
$c\le 1$. Using this fact, we can write
Now,
$g\left(x\right)=rx\left(\right)open="["\; close="]">1\frac{x}{k}$ is a function of
x only. It has a maximum value at
$x=\frac{k(r+\gamma )}{2(rk\theta )}$ and the maximum value is
${M}_{\theta}=\frac{k{(r+\gamma )}^{2}}{4(rk\theta )}$. Thus, from (
4), we get
Using the theory of functional differential equation, we have
□
The set
$\mathcal{D}$ defined by
where
${M}_{\theta}=\frac{k{(r+\gamma )}^{2}}{4(rk\theta )}$, is positively invariant and each solution that starts in
$\mathcal{D}$ remains in
$\mathcal{D}$.
2.2. Existence of Equilibria
In this subsection, we analyze the existence and stability of positive solution of the system (
1).
Model (
1) has three equilibrium points viz (i) the trivial equilibrium
${\overline{E}}_{0}(0,0)$, (ii) the predatorfree equilibrium
${\overline{E}}_{1}\left(\right)open="("\; close=")">\frac{rk}{r\theta k},0$, and (iii) the coexisting equilibrium
${\overline{E}}^{*}({\overline{x}}^{*},{\overline{y}}^{*})$, where
Now ${y}^{*}$ is positive if $\theta {\overline{x}}^{*}+\frac{r(k{\overline{x}}^{*})}{k}>0$ which after simple calculations gives ${\overline{x}}^{*}>\frac{rk}{r\theta k}$.
Using the fact that
$0<m<1$ and from above discussion, we have the following result for the existence of the equilibria of the system (
1).
Proposition 1.

(i)
Predatorfree equilibrium of the system (1) is feasible if $\theta <\frac{r}{k}$, and

(ii)
the coexisting equilibrium ${\overline{E}}^{*}$ is feasible if $\frac{c\beta}{a}>\gamma $ and ${\overline{x}}^{*}>\frac{rk}{r\theta k}$.
Remark 1. When the cooperation for attaining refuge (θ) is less than the overall intraspecific competition on logistic growth of prey ($\frac{r}{k}$), the predatorfree equilibrium ${\overline{E}}_{1}$ of the system (1) is feasible, otherwise it will not exist. On the other hand, when the utilization efficiency of predator after successful predation ($\frac{c\beta}{a}$) is greater than the mortality rate of predator in absence of prey (γ), then the coexisting equilibrium ${\overline{E}}^{*}$ will exist.
2.3. Stability Analysis and Hopf Bifurcation
Here, we examine the stability of each biological feasible equilibrium. We perturb the ODE system nearby any equilibrium point
$\overline{E}(x,y)$, and get the system
where
$X={[x,y]}^{T}$ and
J is the Jacobian matrix.
The Jacobian
J for the system (
1) is:
Characteristic equation of the Jacobian matrix is given by
An equilibrium point
$E(x,y)$ is stable if both the roots of (
8) is negative or having negative real parts. The conditions for this is
where
The equilibrium $E(x,y)$ is unstable if any of the above conditions are violated.
Hopf bifurcation will occur if both the roots of (
8) are purely imaginary for a particular value of the model parameter (called bifurcation parameter). Suppose
$\eta $ a bifurcation parameter, then for bifurcation to occur, the following condition must be satisfied:
Condition (i) will give the critical value of the bifurcation parameter and condition (ii) is the transversality condition.
The Jacobian matrix J at the trivial equilibrium ${\overline{E}}_{0}(0,0)$ has a positive eigenvalue r, thus it is always unstable.
At the predatorfree equilibrium
${\overline{E}}_{1}(\overline{{x}_{1}},0)$, the Jacobian matrix is
where
$\overline{{x}_{1}}=\frac{rk}{r\theta k}$.
${\overline{E}}_{1}$ is feasible when $\frac{r}{k}>\theta $. Jacobian matrix has eigenvalues ${\xi}_{1}=r\overline{x}(1\frac{2\overline{x}}{k})+2\theta \overline{x}$ which is negative for $\frac{r}{k}>\theta $, directly follows from the feasibility of ${\overline{E}}_{1}$. Another eigenvalue is ${\xi}_{2}=\gamma +\frac{c\beta (1m)\overline{x}}{(1+a(1m)\overline{x})}$ which is negative if $\theta >\frac{r\gamma rk\left[c\beta \right(1m)a(1m\left)\gamma \right]}{\gamma k}$.
The characteristic equation at the coexisting equilibrium
${\overline{E}}^{*}$ is
where
The eigenvalues are negative or having negative real parts if ${\overline{l}}_{1}^{*}>0$ as ${\overline{l}}_{2}^{*}>0$ holds always.
From the above discussion we have the following theorem for the stability of the three equilibrium of system (
1).
Theorem 2.

(i)
The trivial equilibrium ${\overline{E}}_{0}(0,0)$ is always unstable.

(ii)
${\overline{E}}_{1}$ is feasible if $\frac{r}{k}>\theta $ and is stable if $\theta <\frac{r\gamma rk(1m)[c\beta a\gamma ]}{\gamma k}$, consequently transcritical bifurcation occurs if $\theta =\frac{r\gamma rk(1m)[c\beta a\gamma ]}{\gamma k}$.

(iii)
${\overline{E}}^{*}$ is stable if ${\overline{l}}_{1}^{*}>0$ and unstable for ${\overline{l}}_{1}^{*}<0$. Hopf bifurcation occurs if ${\overline{l}}_{1}^{*}=0$ and $\frac{\partial {l}_{1}^{*}}{\partial \theta}{}_{\theta ={\theta}^{*}}\ne 0$, where ${\theta}^{*}$ is the critical value of the bifurcation parameter θ.
3. The Mathematical Model for Prey Refuge and Intraspecific Competition
The flow chart of this situation is given in
Figure 2. From this we propose the LotkaVolterra type preypredator system including proportional prey refuge and competition as follows:
Where,
$\alpha $ is the competition coefficient among prey population trying to escape predation pressure through refuge and all other parameters are already defined in the previous model (
1).
Remark 2. From the mathematical point of view, model (1) and model (12) is differ only in one term $\alpha {x}^{2}$ in place of $\theta {x}^{2}$. Thus in the analytical results (dynamics) of model (12) can be derived from the analysis of model (1).
From the analysis of the previous section (using Theorem 1) we have derived the following theorem.
3.1. Boundedness of Solutions
Theorem 3.
Each solution of model (12) with positive initial conditions remains nonnegative and bounded in $\mathcal{S}$ for all time $t>0$, where,
where ${M}_{\alpha}=\frac{k{(r+\gamma )}^{2}}{4(r+k\alpha )}.$
3.2. Existence of Equilibria
The model (
12) has three equilibrium points viz (i) the trivial equilibrium
${\tilde{E}}_{0}(0,0)$, (ii) predator free equilibrium
${\tilde{E}}_{1}(\frac{rk}{r+\alpha k},0)$, and (iii) the coexisting equilibrium
${\tilde{E}}^{*}({\tilde{x}}^{*},{\tilde{y}}^{*})$, where
In model (
12),
${\tilde{E}}_{1}$ always exists and
${\tilde{E}}^{*}$ exists for
$\frac{c\beta}{a}>\gamma $ and
$\alpha {\tilde{x}}^{*}+\frac{r(k{\tilde{x}}^{*})}{k}>0$. Last condition gives
${\tilde{x}}^{*}<\frac{rk}{k\alpha +r}$. Thus we have the following result.
Proposition 2.

(i)
Predatorfree equilibrium of the system (12) always exists.

(ii)
The coexisting equilibrium ${\tilde{E}}^{*}$ is feasible if $\frac{c\beta}{a}>\gamma $ and ${\tilde{x}}^{*}<\frac{rk}{r+\alpha k}$.
Remark 3. In model (1), existence of predatorfree equilibrium ${\overline{E}}_{1}$ depends on θ but in model (12) the equilibrium always exists. One major explanation behind this may be cooperation for attaining refuge (θ) may cause behavioral contradiction among the prey population as they compete in almost all aspects except finding refuge, and for that reason, θ becomes a critical factor for the existence of equilibrium.
3.3. Stability and Hopf Bifurcation
The Jacobian is at any steady point
$\tilde{E}(x,y)$ of system (
12) is given below:
Using the stability analysis of the previous section (using Theorem 2, we can establish the following theorem for the stability of the equilibrium points of the system (
12).
Theorem 4.

(i)
The trivial equilibrium ${\tilde{E}}_{0}(0,0)$ is always unstable.

(ii)
${\tilde{E}}_{1}$ is stable if $\alpha <\frac{rk(1m)[c\beta a\gamma ]r\gamma}{\gamma k}$, and unstable for $\alpha >\frac{rk(1m)[c\beta a\gamma ]r\gamma}{\gamma k}$ consequently transcritical bifurcation occurs at $\alpha =\frac{rk(1m)[c\beta a\gamma ]r\gamma}{\gamma k}$.

(iii)

${\tilde{E}}^{*}$ is stable if ${\tilde{l}}_{1}^{*}>0$ and unstable for ${\tilde{l}}_{1}^{*}<0$. Hopf bifurcation occurs if ${\tilde{l}}_{1}^{*}=0$ and $\frac{\partial {\tilde{l}}_{1}^{*}}{\partial \alpha}{}_{\alpha ={\alpha}^{*}}\ne 0$, where ${\alpha}^{*}$ denotes the critical value of the bifurcation parameter α. In this case,
4. Global Stability of ${\mathbf{E}}^{*}$
The following theorem establish the global stability of the interior equilibrium ${E}^{*}$.
Theorem 5.
The interior equilibrium ${E}^{*}$ is globally asymptotically stable if the condition mentioned below are satisfied:
Proof. We construct the suitable Lyapunov function
V as,
It is easy to verify that
V is positive definite for all
$(x,y)\in {R}_{+}^{2}\setminus ({x}^{*},{y}^{*})$. Now we take the derivative of
V along the solutions of the system (
1) and obtain the following,
Rearranging the equation (18) we get, following,
Thus, if the conditions given in (
16) holds, then the right hand side of (18) is negative definite. Thus
$\frac{dV}{dt}<0$, along all the trajectories in the first quadrant except
$({x}^{*},{y}^{*})$. Also we have
$\frac{dV}{dt}{}_{{E}^{*}}=0$. The proof follows from the suitable Lyapunov function
V and LyapunovLaSalle’s invariant principle (cf. Hale [
39]). Hence, the interior equilibrium point
${E}^{*}$ is globally asymptotically stable if the condition (
16) is satisfied. □
5. Direction of Hopf Bifurcation
In this section, we analyse the stability nature and direction of the Hopf bifurcating periodic solutions. For this, we consider the transformation,
${u}_{1}=x{x}_{*}$,
${u}_{2}=y{y}_{*}$, which transform the system (
1) to the following system,
where
${a}_{ij}$ and
${b}_{ij}$, are determined by using the following relations,
where,
${a}_{1}=\frac{\partial {F}_{1}}{\partial x}$,
${a}_{2}=\frac{\partial {F}_{2}}{\partial y}$, etc., and
We get the following relations from Hassard et al. [
21],
where
We calculate the values of
${g}_{20}$,
${g}_{02}$,
${g}_{11}$ and
${g}_{21}$ by comparing the coefficients of
${z}^{2}$,
${\overline{z}}^{2}$,
$z\overline{z}$ and
${z}^{2}\overline{z}$. Hence
where,
The stability and direction can be determined from the following theorem [
28,
38].
Theorem 6. The nature of the bifurcating periodic solutions from the values of ${\beta}_{2},\phantom{\rule{3.33333pt}{0ex}}{\mu}_{2}$ and ${\tau}_{2}$ in the center manifold at the threshold value ${\alpha}^{*}$.
(i) If ${\mu}_{2}$ is positive (negative), then the system (1) emits a stable (unstable) limit cycle for $\alpha >{\alpha}^{*}$ ($\alpha <{\alpha}^{*}$).
(ii) ${\beta}_{2}$ determines the stability of the bifurcating periodic solutions: the bifurcating periodic solutions are stable (unstable) if ${\beta}_{2}<0$$({\beta}_{2}>0)$ and
(iii) ${\tau}_{2}$ determines the period of periodic solutions: the period increases (decreases) if ${\tau}_{2}>0$$({\tau}_{2}<0)$.
A similar theorem can be obtained for the parameter $\theta $ from Theorem 6.
6. Numerical Simulations
In this section, we provide numerical examples for visualizing the dynamical regimes that can be seen in our model systems analysed in the previous sections. Parameter values with detailed descriptions and units used for the simulation are given in
Table 1. The value of carrying capacity is taken as
$k=4$ assuming that the maximum number of 4 prey can be supported within a
$3m\times 3m$ square area. Though conversion efficiency (
c) varies within a range of
$01$, we have considered
$c=0.75$ for the simulations. The results obtained from simulations are discussed in the following subsections.
6.1. Numerical Simulation of System (1)
Here, simulation of system (
1) is discussed. In
Figure 3, we have plotted the transcritical forward bifurcation of
${\overline{E}}_{1}$ and existence of
${\overline{E}}^{*}$ by varying parameter
$\theta $. For some large value of theta,
${\overline{E}}^{*}$ becomes unstable (Theorem 2).
In
Figure 4, we have seen limit cycle for
$m=0$, (other parameters are given in
Table 1) i.e. limit cycle can occur if refuge is neglected. But if we take
$m=0.3$ (i.e.,
$30\%$ of the total prey population are in refuge habitat and rest
$70\%$ are vulnerable to predation) the limit cycle disappears and coexisting equilibrium becomes a stable focus.
In
Figure 5, we have taken
$m=0.3,\phantom{\rule{3.33333pt}{0ex}}\theta =0$ (i.e., there is no cooperation among the prey population for finding refuge habitat) the system is stable. Then we increase the value of
$\theta $ (i.e., increasing cooperation among the prey individuals for finding safe refuge place) and seen limit cycle for
$\theta =0.005$ (other parameters are given in
Table 1). From this we can conclude that the impact of
$\theta $ is destabilising the stable steady state.
Figure 6 is plotted for three different initial values of the model variables. This shows the nonlinear stability of coexisting equilibrium
${E}^{*}$.
Figure 7 is also plotted for two different initial values of the model variables. This figure shows that the bifurcating periodic solutions are stable and Hopf bifurcation is of supercritical type. In
Figure 8, we have plotted the Hopf bifurcation diagram taking
$\theta $ as the main parameter. Limit cycle appear as the parameter
$\theta $ crosses the critical value (Theorem 2).
In
Figure 9, we have plotted the stability regions of coexisting equilibrium in
$\theta \beta $ and
$\theta m$ parameter planes. In
Figure 9a, it is shown how parameters
$\theta $ and
$\beta $ critically affect the stability of the system. We have observed that for the very low value of
$\theta $,
${\overline{E}}^{*}$ is not feasible but
${\overline{E}}_{1}$ is stable and as
$\theta $ crosses its critical value
${\overline{E}}^{*}$ becomes stable and the stability of
${\overline{E}}^{*}$ depends on the combined value of
$\theta $ and
$\beta $ as higher values can make the stable system unstable via Hopf bifurcation. In
Figure 9b direct relationship between parameter
$\theta $ and
m to the stability of the system has shown, where with the increased value of
$\theta $ only a very high value of
m can make the system stable.
From the simulations, we observe that
$\theta $ is one of the critical parameters for the model system (
1), stability of different equilibria depends on it. For very low value of
$\theta $,
${\overline{E}}^{*}$ is not feasible but
${\overline{E}}_{1}$ is stable and as
$\theta $ crosses its critical value
${\overline{E}}^{*}$ becomes stable via transcritical bifurcation and finally further higher values can make the
${\overline{E}}^{*}$ unstable via Hopf bifurcation.
6.2. Numerical Simulation of System (12)
In
Figure 10, we have plotted the transcritical forward bifurcation of
${\tilde{E}}_{1}$ and existence of
${\tilde{E}}^{*}$ varying the parameter
$\alpha $ (Theorem 4).
Hopf bifurcation diagram is plotted in
Figure 11 taking
$\alpha $ as the main parameter. The Limit cycle disappears as the parameter
$\alpha $ crosses the critical value
${\alpha}^{*}$
In
Figure 12, we have taken
$m=0.3$,
$\alpha =0.001$ (i.e. there is very less competition among the prey individuals for finding refuge habitat) the system is stable. Now, if we increase the value of
$\alpha $ to 0.01 (i.e., the prey individuals become ten times more competitive in the process of finding a safe refuge place), the system stabilizes faster (other parameters are given in
Table 1). From this, we can conclude that
$\alpha $ has stabilizing effect on the system.
In
Figure 13, we have plotted the stability regions of coexisting equilibrium in
$\alpha \beta $ and
$\alpha m$ parameter planes. The area of stability regions increased as both the value of
m and
$\alpha $ are increased. But
$\beta $ can make the stable system unstable via Hopf bifurcation. The critical value of
$\beta $, i.e.,
${\beta}^{*}$ depends on the other model parameters such as
$\alpha $. Here we have shown that the critical value of
$\beta $ significantly depends on
$\alpha $. In (b) we have taken higher value of
$\beta $, then we vary
$\alpha $ and
m. We have observed that
${\tilde{E}}^{*}$ changes its stability to become stable, then is not feasible for higher values of
m and
$\alpha $, but
${\tilde{E}}_{1}$ is stable there.
Thus from the numerical simulations, we observe that for the model system (
12), the stability of different equilibria depends on the parameter
$\alpha $. Here, for very low value of
$\alpha $,
${\tilde{E}}^{*}$ is unstable and it becomes stable as
$\alpha $ crosses the critical value via Hopf bifurcation but for very high value of
$\alpha $${\tilde{E}}^{*}$ is not feasible where
${\tilde{E}}_{1}$ becomes stable via transcritical bifurcation.
7. Discussion and Conclusion
In HooglyMatla estuarine mangrove ecosystem, there are several islands and those are crisscrossed by numerous creeks which originate from the main rivers. These creeks have rich detritus loads which come from the adjacent mangrove litter. Therefore, these creeks support many detritivorous fishes and one of such was described by Sarwardi et al. [
29]. In their study (Sarwardi et al., 2013) described one important detritivorous fish of these creeks, namely
Liza parsia which lives in the creeks and its predator
Otolithoides pama visits the creeks only during high tide for capturing food. The mangrove plants in this region are extended from the supralittoral zone to the lowerlittoral zone up to the creek bed. During high tide as the water level comes up the lower parts of the plants submerge under the water and the bushy part of the submerged mangrove plants creates shelter for
Liza parsia to avoid its predator. So, this spatial refuge is possible due to enhanced environmental heterogeneity in this ecosystem. But there is not enough refuge habitats to support all the individuals of the prey population and if it’s the case then the refuge decision among the prey population is made either though cooperation or by competition and our study is focused to find out the impact of these two intraspecific interaction in the system dynamics.
Preypredator system with prey refuge has been a soughtafter topic among ecological and mathematical modelers for almost the last five decades. Previous studies have assured that refuge has a stabilizing effect on preypredator dynamics. Several models are already established with various conditions and situations. Still, another finding was that if the refuge class or fraction becomes too large, the predator faces food scarcity. The system becomes destabilized due to the extinction of the predator. But in any natural system, all members of a prey population can’t avoid predation at any given time; thus, continuous decision pressure works among them. The refuge among the prey population can happen in two ways – cooperatively or competitively, and these two inevitable intraspecific interactions must affect preypredator dynamics. This study finds the impact of cooperation and intraspecific competition among prey populations trying to escape predation through the refuge on the stability of preypredator dynamics.
For the two situations mentioned above, in this paper, we have considered a preypredator system incorporating a prey refuge and assumed that the predator response function is of Holling type II. Cooperation (model (
1)) and intraspecific competition (model (
12)) among the prey population for escaping predation pressure through refuge are incorporated into the system to find out the process (either cooperation or intraspecific competition) for leading the population towards better stabilization.
Firstly, the stabilizing effect of refuge on the system dynamics has been shown (
Figure 4), and after the addition of refuge, the system dynamics settle from limit cycle behavior to a stable equilibrium. The refuge protects the prey population from predators, and it reduces predation pressure. Overexploitation of prey population by predator turns the prey population into a divergent oscillation state. But this vulnerable state is transformed into a stable state with the help of the prey population’s refuge strategies. The stable state of the prey population also makes a stable predator population.
When the decision of refuge is made cooperatively among the prey population, the stable equilibrium again goes back to a limit cycle condition which clearly denotes that cooperation has a destabilizing effect on the system dynamics during refuge (
Figure 5). From the prey’s point of view, if the refuge is attained cooperatively, then maximum individuals can escape predation easily, which must be favorable for them but on the other hand predator will suffer from low predation success, and in an extreme situation, they will be eliminated from the system if the system is linear as a predator has no alternative option to predate.
Our study shows that the system reaches faster in the equilibrium state when the refuge facility is achieved through intraspecific competition among the prey population (
Figure 12). Refugees during intraspecific competition have a stabilizing effect on the system dynamics. As intraspecific competition has a negative impact on individuals of the prey population, the refuge class is constantly checked and never too large to cause prey scarcity for a predator. A controlled amount of the prey population will always escape predation through the refuge, which helps the prey population for better survival and the predator population will not go extinct in any situation as they never face a crisis for food. Due to this balanced interaction, the system dynamics reach equilibrium faster.
In reality, the decisionmaking process to escape predation pressure by refuge either through competition or cooperation within a prey population is complex and doesn’t follow any hard and fast rules. Cooperation or competition in various prey organisms depends on the spatial heterogeneity of their habitat as when the shelter place (refuge size) is large enough to protect all prey individuals, then only cooperation acts as the deciding factor, but any scarcity in protected habitat turns the cooperative behavior into a competitive one.
It is more difficult for a prey population to escape predation through refuge individually than living in a group and cooperating conspecifically. Group living in prey population has many facilities from the perspective of refuge as it reduces the chance of detection by the predator due to the presence of many eyes and ears tuned to predation threats [
34]. So, cooperation plays a major role in attaining refuge within a group living prey population, but not necessarily always it is the case. An underline intraspecific competition may be found among the group members, which can be explained by the selfishherd effect [
35]. Though each member enjoys the benefit of decreasing the domain of danger when living in a group than as a solitary forager, still those individuals who are present outside the group have a larger domain of danger compared to those at the core, and thus, they may compete each other to take shelter at the core of the group [
34].
Interactions of the living world are very complex, and difficult to conclude whether competition or cooperation is helpful for refugees. A high rate of parental care has been observed in many bird species, including the protection of offspring from a predator by forming nests. In such cases, the parents cooperate with their offspring for successful refuge, but at the same time, the offspring may compete to achieve a successful refuge. In reality, when prey populations are small, they need cooperation rather than competition, but the reverse situation is noticed in the case of the large prey population. Therefore, both cooperation and competition are inevitable processes during the refuge of the prey population, but the selection of either cooperation or competition depends on the situation that nature wants.
Conflicts of Interest
Authors declare that there is no conflict of interest.
References
 Baldwin, I. T. Inducible defenses and population biology. Tree 1996, 11, 104–105. [Google Scholar] [CrossRef] [PubMed]
 Chakraborty, S.; Bhattacharya, S.; Feudel, U.; Chattopadhyay, J. The role of avoidance by zooplankton for survival and dominance of toxic phytoplankton. Ecol Compl 2012, 11, 144–153. [Google Scholar] [CrossRef]
 GonzalezOlivares, E.; RamosJiliberto, R. Dynamic consequences of prey refuges in a simple model system: more prey, fewer predators and enhanced stability. Ecol Model 2003, 166, 135–146. [Google Scholar] [CrossRef]
 Abrams, P. A. Measuring the impact of dynamic antipredator traits on predatorpreyresource interactions. Ecology 2008, 89, 1640–1649. [Google Scholar] [CrossRef]
 Zhuang, K.; Wen, Z. Dynamical behaviors in a discrete predatorprey model with a prey refuge. Int J Comput Math Sci 2011, 5, 196–197. [Google Scholar]
 Cronin, W. T. The Visual Ecology of PredatorPrey Interaction. Ecology of PredatorPrey Interaction; Oxford University Press, 2005; pp. 105–138. [Google Scholar]
 Anderson, T. W. Predator responses prey refuges and densitydependent mortality of a marine fish. Ecology 2001, 82, 245–257. [Google Scholar] [CrossRef]
 Lotka, A.J. Elements of Physical Biology; Williams and Wilkins: Baltimore, 1925. [Google Scholar]
 Volterra, V. Fluctuations in the abundance of a species considered mathematically. Nature 1926, 118, 558–560. [Google Scholar] [CrossRef]
 Gause, G. F. The Struggle for Existence; Williams and Wilkins: Baltimore, 1934. [Google Scholar]
 Holling C., S. The components of predation as revealed by a study of small mammal predation of the European pine sawfly. Can. Entomol. 1959, 91, 293–320. [Google Scholar] [CrossRef]
 Holling, C.S. Some characteristics of simple types of predation and parasitism. Can. Entomol. 1959, 91, 385–395. [Google Scholar] [CrossRef]
 Holling, C.S. The functional response of predators to prey density and its role in mimicry and population regulation. Mem Entomol Soc Ca 1965, 45, 3–60. [Google Scholar] [CrossRef]
 Rosenzweig, M. L.; MacArthur, R. H. Graphical representation and stability conditions of predatorprey interactions. Am. Nat. 1963, 97, 209–223. [Google Scholar] [CrossRef]
 Gkana, A.; Zachilas, L. Incorporating prey refuge in a prey–predator model with a Holling type I functional response: random dynamics and population outbreaks. J Biol. Phys. 2013, 39, 587–606. [Google Scholar] [CrossRef] [PubMed]
 Das, A.; Samanta, G.P. Stochastic preypredator model with additional food for predator. Physica A: Statistical Mechanics and its Applications 2018, 512, 121–141. [Google Scholar] [CrossRef]
 Das, A.; Samanta, G.P. Modelling the fear effect on a stochastic preypredator system with additional food for predator. Journal of Physics A: Mathematical and Theoretical 2018, 51, 465601. [Google Scholar] [CrossRef]
 Das, A.; Samanta, G.P. A preypredator model with refuge for prey and additional food for predator in a fluctuating environment. Physica A: Statistical Mechanics and its Applications 2020, 538, 122844. [Google Scholar] [CrossRef]
 Das, A.; Samanta, G.P. Influence of environmental noises on a preypredator species with predatordependent carrying capacity in alpine meadow ecosystem. Mathematics and Computers in Simulation 2021, 190, 1294–1316. [Google Scholar] [CrossRef]
 Hassell, M.P. The dynamics of arthropod predatorprey systems. Princeton University Press, 1978. [Google Scholar]
 Hassard, B.; Kazarinof, D.; Wan, Y. Theory and Applications of Hopf Bifurcation; Cambridge University Press: Cambridge, 1981. [Google Scholar]
 Cressman, R.; Garay, J. A predatorprey refuge system: Evolutionary stability in ecological systems. Theoretical Population Biology 2009, 76, 248–257. [Google Scholar] [CrossRef]
 Chen, L.; Chen, F.; Chen, L. Qualitative analysis of predatorprey model with Holling type II functional response incorporating a constant prey refuge. Nonlinear Analysis: Real World Applications 2010, 11, 246–252. [Google Scholar] [CrossRef]
 Maynard Smith, J. Models in ecology; Cambridge Univ. Press, 1974. [Google Scholar]
 Mukherjee, D. The effect of refuge and immigration in a predatorprey system in the presence of a competitor for the prey. Nonlinear Analysis: Real World Applications 2016, 31, 277–287. [Google Scholar] [CrossRef]
 Eisenberg, J.N.S.; Washburn, J.O.; Schreiber, S.J. Generalist feeding behaviour of Aedes sierrensis larvae and their effects on protozoan populations. Ecology 2000, 81, 921–935. [Google Scholar] [CrossRef]
 Kar, T. K. Stability analysis of a preypredator model incorporating a prey refuge. Communication in Nonlinear Science and Numerical Simulatio 2005, 10, 681–691. [Google Scholar] [CrossRef]
 Sarwardi, S.; Mandal, P.K.; Ray, S. Analysis of a competitive prey–predator system with a prey refuge. Biosystems 2012, 110, 133–148. [Google Scholar] [CrossRef] [PubMed]
 Sarwardi, S.; Mandal, P.K.; Ray, S. Dynamical behaviour of a twopredator model with prey refuge. J. Biol Phys 2013, 39, 701–722. [Google Scholar] [CrossRef] [PubMed]
 Roy, U.; Sarwardi, S.; Majee, N. C.; Ray, S. Effect of salinity and fish predation on zooplankton dynamics in Hooghly–Matla estuarine system, India. Ecological informatics 2016, 35, 19–28. [Google Scholar] [CrossRef]
 Kumar, G.R.; Ramesh, K.; Lakshminarayan, K.; Rao, K.K. Hopf Bifurcation and Stochastic Stability of a PreyPredator Model Including Prey Refuge and Intraspecific Competition Between Predators. Int. J. Appl. Comput. Math 2022, 8, 209. [Google Scholar] [CrossRef]
 Wang, Y.; Wang, J. Influence of prey refuge on predatorprey dynamics. Nonlinear Dyn 2012, 67, 191–201. [Google Scholar] [CrossRef]
 Johnson, W. D. Predation, habitat complexity and variation in density dependent mortality of temperate reef fishes. Ecology 2006, 87, 1179–1189. [Google Scholar] [CrossRef]
 Beauchamp, G. Group Living. In Encyclopedia of animal behavior; Breed, M.D., Moore, J., Eds.; Academic Press: Oxford, 2010; pp. 21–24. [Google Scholar]
 Hamilton, W.D. Geometry for the selfish herd. Journal of theoretical Biology 1971, 31, 295–311. [Google Scholar] [CrossRef]
 Zhang, Y.; Huang, J. Bifurcation analysis of a predator–prey model with Beddington–DeAngelis functional response and predator competition. Mathematical Methods in the Applied Sciences. 2022. [Google Scholar] [CrossRef]
 Haque, M.; Sarwardi, S. Dynamics of a harvested prey–predator model with prey refuge dependent on both species. International Journal of Bifurcation and Chaos 2018, 28, 1830040. [Google Scholar] [CrossRef]
 Kuznetsov, Y.A. Elements of Applied Bifurcation Theory; Springer: New York, NY, USA, 2004. [Google Scholar]
 Hale, J. Ordinary differential equation; Klieger Publishing Company: Malabar, 1989. [Google Scholar]

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content. 
© 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).