Impact of Cooperation and Intra-Speciﬁc Competition of Prey on the Stability of Prey–Predator Models with Refuge

: The main objective of this study is to ﬁnd out the inﬂuences of cooperation and intra-speciﬁc competition in the prey population on escaping predation through refuge and the effect of the two intra-speciﬁc interactions on the dynamics of prey–predator systems. For this purpose, two mathematical models with Holling type II functional response functions were proposed and analyzed. The ﬁrst model includes cooperation among prey populations, whereas the second one incorporates intra-speciﬁc competition. The existence conditions and stability of different equilibrium points for both models were analyzed to determine the qualitative behaviors of the systems. Refuge through intra-speciﬁc competition has a stabilizing role, whereas cooperation has a destabilizing role on the system dynamics. Periodic oscillations were observed in both systems through Hopf bifurcation. From the analytical and numerical ﬁndings, we conclude that intra-speciﬁc competition affects the prey population and continuously controls the refuge class under a critical value, and thus, it never becomes too large to cause predator extinction due to food scarcity. Conversely, cooperation leads the maximal number of individuals to escape predation through the refuge so that predators suffer from low predation success.


Introduction
Prey-predator interaction has immense importance in the existence of life.In predatorprey systems, the prey population uses different defense mechanisms to escape from predators [1,2].Among them, refuge is the most common in the literature [3,4], and it has become a significant issue in theoretical ecology [5][6][7].One famous proverb about the Fox-Rabbit chase says, "The fox is only running for his dinner, but the rabbit is running for his life [8]".If the rabbit can successfully save its life, it is called the rabbit escaping predation from the fox through refuge.Thus, the refuge can be defined as any strategy that decreases the rate of predation [9,10].
Refuge can be attained in two ways-(i) constant numerical refuge, where a definite number of individuals in the population escape from predation, and (ii) proportional refuge, where a fixed proportion of individuals in the population escape from predation.In both cases, the whole prey population cannot 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 individuals in the population can enjoy the facility of refuge due to the availability of a fixed number of safe refuge places, two habitats are available for the prey population-either they live in the refuge habitat or they roam around in the open habitat where they are available to predators [11].
In [12], the authors considered a prey-predator model with a constant numerical prey refuge and showed that there are two major factors that control this type of refuge in the population: the population density of the prey (x) and the available refuge (M) (which is the maximum number of prey individuals that can be supported in the refuge habitat).When x ≤ M, competition among prey for refuge is negligible, since prey fitness depends only on predation [11].Under such conditions, cooperation among individuals in the prey population must be seen.But, this is not the case at all times, because this type of refuge can include a definite number of individuals, and thus, they cannot move freely between habitats when the population size exceeds the refuge number (x > M), and obviously, intraspecific competition to take refuge will increase proportionally as the density increases.The conclusions regarding this type of refuge are that it has a strong stabilizing effect [13] and it can prevent the extinction of prey, but when the population is too large for all members to stay all in the refuge, coexistence with predators is expected, as the predators will obtain continuous resources [11].So, in that case, the mortality rate of the prey increases when the population density goes beyond the refuge capacity, and this will cause a negative feedback stabilizing effect on the system [14].
In case of a proportional refuge, there are also two habitats available for the prey population-either they live in the refuge habitat or they roam around in the open habitat, where they are available to predators [11].But, under such conditions, movement of prey individuals always occurs between the open and refuge habitats [15].In [16], Kar proposed the Lotka-Volterra type prey-predator system, including proportional prey refuges, where m is denoted as the fractional refuge of prey from predators (which is different from the previously discussed refuge, M), and it ranges from 0 < m < 1.Now, if we consider the whole prey population as x, then the mx portion of the prey is in the refuge habitat, and the (1 − m)x portion of the prey is in open habitats and is thus available to the predators.For the proportional refuge, the whole population never has the opportunity to attain refuge, and thus, the process of finding a refuge place is either cooperation-or competitionbased.This type of refuge can reduce the mortality rate of prey by reducing the predation pressure and can affect positively affect the growth of prey and negatively affect that of predators.Thus, in this way, refuge can act as a stabilizing effect in prey-predator dynamics [17].Extreme refuge may destabilize prey-predator dynamics, as the predator population becomes extinct by suffering from lower predation success due to the large refuge size [18].
Mathematically, prey-predator interaction was first described by Alfred James Lotka in 1925 [19] and Vito Volterra in 1927 [20].In 1934, Georgy Frantsevich Gause [21] first found refuge in his experiment with the protozoan predator Didinium feeding on Paramecium.One major drawback of the Lotka-Volterra model is that the predator population never becomes saturated.Crawford Stanley Holling solved this problem [22][23][24] by suggesting three kinds of functional responses (Holling types I, II, and III) for different species of predators.In 1963, Rosenzweig and MacArthur [25] first replaced the linear functional response of the Lotka-Volterra model with the Holling type II response function and first documented that prey-predator coexistence is not limited to a stable equilibrium; a limit cycle appears when the stable equilibrium undergoes Hopf bifurcation [26].Since then, various modeling approaches for the prey-predator system, including refuge, have been conducted by researchers [16,[27][28][29].
The stabilizing role of refuge in prey-predator systems was shown many years ago, and since then, various aspects of prey-predator systems with refuge have been analyzed using Holling type II nonlinearities [30,31].Mapunda et al. [32] and Kumar et al. [30] showed the effect of intra-specific competition among predators on the prey-predator system in their respective studies.The effect of cooperation among predators on the preypredator system through cooperative hunting strategies was shown by Saha and Samanta in their study [33].The effect of inter-specific competition among predator populations in a three-species system was analyzed by Sarwadi et al. [29,34].The impact of cooperation among two prey species in a three-species prey predator system was analyzed by Rihan et al. [35].The present study fills the research gap by dealing with the impacts of cooperation and intra-specific competition among the prey population when escaping from predators through refuge on the prey-predator system, and it is found that when the prey individuals become competitive with each other for shelter in the refuge habitat, the system dynamics achieve stability faster.On the other hand, any cooperation among prey populations regarding refuge turns the stable equilibrium condition into a limit cycle condition.
In previous articles ([29,30,32,34] and the references therein), authors have considered competition among predator populations.To our knowledge, this is the first article to consider the impacts of cooperation and intra-specific competition within the prey population on attaining refuge.The key question of this study is "how are the dynamics of the prey-predator system impacted when the decision about refuge is made through cooperation and the same for intra-specific competition?"To obtain an answer to this problem, we developed two mathematical models.The first model includes cooperation among prey populations, whereas the second model contains intra-specific competition.The dynamics of the models were analyzed with qualitative theories.
This study is presented sequentially in the latter sections as follows: Section 2 deals with the mathematical model formation, including both the prey refuge and the cooperation among prey.The local stability and bifurcation analysis of the equilibrium points is provided.The global stability of the co-existing equilibrium is analyzed in Section 2.4.Another mathematical model considering the prey refuge and the intra-specific competition among prey is formed.The existence of Hopf bifurcation and a stability analysis are conducted in Section 3. The direction of Hopf bifurcation is analyzed in Section 3.4.Detail numerical simulations using Matlab are delineated in Section 4. The final discussion and biological interpretation of the results are included in the conclusions in Section 5.

Mathematical Model for Prey Refuge and Cooperation
Here, we propose the Lotka-Volterra type prey-predator system, including proportional prey refuge and cooperation, which makes the following assumptions: There are two populations in the model.x(t) is the prey population size and y(t) 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.γ is the death rate of the predator in the absence of prey, β is the predation rate, and c is the conversion factor that denotes the efficiency of predators for each captured prey.Also, β a is the maximum number of prey that the predator can eat per unit time; 1 a is the half-saturation constant.
The term βx 1+ax denotes the predator's response to prey.This type of predator response function is known as a Holling type II functional response.
Let m be the fractional refuge of prey from the predators, where 0 < m < 1 and is constant.So, the mx(t) portion of the prey is in the refuge habitat, and the (1 − m)x(t) portion of the prey is in open habitats and is thus available to the predator, and θ is the cooperation coefficient among the prey population trying to escape predation pressure through refuge.
Using the assumptions discussed above, we obtain the following mathematical model: with the initial biological condition It is assumed that all system parameters are positive constants.A 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 the biological feasibility of the model (1).  1) for the Lotka-Volterra type prey-predator system, including the proportional prey refuge and cooperation among the prey population trying to escape predation pressure through refuge.

Boundedness of Solutions
Biologically, all populations in the model should be finite.For this, a boundedness analysis is significant.We derived the following theorem for the boundedness of the solutions to the model (1).

Theorem 1.
All solutions to the model (1) with the given initial conditions (2) are non-negative and bounded for all time periods t > 0.
Proof.We rewrite the first equation of (1), where Since x(0) > 0, we can conclude that x(t) is always positive.Similarly, the nonnegativity of y(t) can be established.
For the boundedness analysis, we define W = x + y.Now, adding the two equations of (1), we can obtain By definition, c ≤ 1.Using this fact, we can write 2(r−kθ) , and the maximum value is M θ = k(r+γ) 2 4(r−kθ) .Thus, from (4), we obtain Using the theory of functional differential equation, we have The set D is defined by where M θ = k(r+γ) 2 4(r−kθ) , is positively invariant and each solution that starts in D remains in D.

Existence of Equilibria
In this subsection, we analyze the existence and stability of the positive solution to the system (1).
Using the fact that 0 < m < 1 and from the above discussion, we have the following result for the existence of the equilibria of the system (1).

Proposition 1.
(i) The predator-free equilibrium of the system (1) is feasible if θ < r k , and (ii) the coexisting equilibrium Ē * is feasible if cβ a > γ and x * > rk r−θk .
Remark 1.When the cooperation for attaining refuge (θ) is less than the overall intra-specific competition in the logistic growth of prey ( r k ), the predator-free equilibrium Ē1 of the system (1) is feasible; otherwise, it will not exist.On the other hand, when the utilization efficiency of predators after successful predation ( cβ a ) is greater than the mortality rate of predators in the absence of prey (γ), then the coexisting equilibrium Ē * will exist.

Stability Analysis and Hopf Bifurcation
Here, we examine the stability of each biologically feasible equilibrium.We perturb the ODE system near to any equilibrium point Ē(x, y) and obtain the system where X = [x, y] T and J is the Jacobian matrix.
The Jacobian J for the system (1) is A characteristic equation of the Jacobian matrix is given by An equilibrium point E(x, y) is stable if both the roots of ( 8) are negative or have negative real parts.The conditions for this are 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 the bifurcation parameter).Suppose η is a bifurcation parameter.Then, for bifurcation to occur, the following condition must be satisfied: (ii) 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 Ē0 (0, 0) has a positive eigenvalue r; thus, it is always unstable.
At the predator-free equilibrium Ē1 ( x1 , 0), the Jacobian matrix is where x1 = rk r−θk .Ē1 is always feasible as r k > θ is true for the model (1).The Jacobian matrix has eigenvalues ξ 1 = r x(1 − 2 x k ) + 2θ x and is negative for r k > θ, directly following on from the feasibility of Ē1 .Another eigenvalue is The characteristic equation at the coexisting equilibrium Ē * is The eigenvalues are negative or have negative real parts if l * 1 > 0 as l * 2 > 0 always holds.
From the above discussion, we have the following theorem for the stability of the three equilibriums of system (1).
(ii) Ē1 is always feasible and is stable if where θ * is the critical value of the bifurcation parameter θ.

Global Stability of E *
The following theorem establishes the global stability of the interior equilibrium E * .
Theorem 3. The interior equilibrium E * is globally asymptotically stable if the conditions mentioned below are satisfied: (ii) Proof.We construct the suitable Lyapunov function V as It is easy to verify that V is positive definite for all (x, y) ∈ R 2 + \ (x * , y * ).Now, we take the derivative of V along the solutions of the system (1) and obtain the following: By rearranging equation ( 14) we obtain the following: Thus, if the conditions given in (12) hold, then the right-hand side of ( 14) is negative definite.Thus, dV dt < 0 along all trajectories in the first quadrant except (x * , y * ).Also, we have dV dt | E * = 0.The proof follows on from the suitable Lyapunov function V and Lyapunov-LaSalle's invariant principle (cf.Hale [36]).Hence, the interior equilibrium point E * is globally asymptotically stable if condition (12) is satisfied.

The Mathematical Model for Prey Refuge and Intra-Specific Competition
The flow chart of this situation is given in Figure 2. From this, we propose the Lotka-Volterra type prey-predator system, including the proportional prey refuge and competition, as follows: where α is the competition coefficient among the prey population trying to escape predation pressure through refuge, and all other parameters were already defined in the previous model (1).From the analysis presented in the previous section (using Theorem 1) we derived the following theorem.

Boundedness of Solutions
Theorem 4. Each solution to model (14) with positive initial conditions remains non-negative and bounded in S for all time points t > 0, where where M α = k(r+γ) 2 4(r+kα) .
Remark 3. In model (1), the existence of the predator-free equilibrium Ē1 depends on θ, but in model (14), the equilibrium always exists.One major explanation for the cooperation for attaining refuge (θ) may cause behavioral contradictions 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 the equilibrium.

Stability and Hopf Bifurcation
The Jacobian at any steady point Ẽ(x, y) of system ( 14) 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 in the system (14).

Direction of Hopf Bifurcation
In this section, we analyze 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) into the following system, where a ij and b ij are determined by using the following relations: where a 1 = ∂F 1 ∂x , a 2 = ∂F 2 ∂y , etc., and We obtain the following relations from Hassard et al. [37], where Thus, We calculate the values of g 20 , g 02 , g 11 and g 21 by comparing the coefficients of z 2 , z2 , z z and z 2 z.Hence, where The stability and direction can be determined from the following theorem [29,38].

Remark 4.
A similar theorem can be obtained for the parameter θ from Theorem 6.

Numerical Simulations
In this section, we provide numerical examples for visualizing the dynamical regimes that can be seen in our model systems analyzed in the previous sections.Parameter values with detailed descriptions and the units used for the simulation are given in Table 1.The value of the carrying capacity is taken as k = 4, assuming that the maximum number of four prey can be supported within a 3 m × 3 m square area.Although the conversion efficiency (c) varies within a range of 0-1, we considered c = 0.75 for the simulations.The results obtained from simulations, using the ode45 solver in Matlab are discussed in the following subsections.

Numerical Simulation of System (1)
In this subsection, a simulation of the system (1) is presented in Figures 3-8 based on the calculations made in Section 2. In Figure 3, we plot the transcritical forward bifurcation between Ē1 and E * by varying the parameter θ.For some large values of θ, Ē * becomes feasible (according to Theorem 2).In region A, Ē * is not feasible, but Ē1 is stable there; in region B, Ē * is stable; in Region C, Ē * is unstable.In (a), m = 0.3 and in (b), β = 0.18, the other parameters are as given in Table 1.
In Figure 4, we can see the limit cycle for m = 0 (other parameters are given in Table 1), i.e., the limit cycle can occur if refuge is neglected.But, if we take m = 0.3 (i.e., 30% of the total prey population is in the refuge habitat and the other 70% is vulnerable to predation), the limit cycle disappears, and the coexisting equilibrium E * becomes a stable focus.
In Figure 5, we use m = 0.3 and θ = 0 (i.e., there is no cooperation among the prey population for finding a refuge habitat), and the system is stable.Then, we increase the value of θ (i.e., increase cooperation among the prey individuals for finding a safe refuge place) and set the limit cycle to θ = 0.005 (other parameters are given in Table 1).From this, we can conclude that the impact of θ is destabilization of the stable, steady state.
Figure 6 plots three different initial values of the model variables.This shows the nonlinear stability of the coexisting equilibrium E * .In Figure 7, we plot the Hopf bifurcation diagram, taking θ as the main parameter.A stable limit cycle appears, as the parameter θ crosses the critical value (Theorem 2).
In Figure 8, we plot the stable regions of the coexisting equilibrium in the θ-β and θ-m parameter planes.Figure 8a shows how parameters θ and β critically affect the stability of the system.We can observe that for the very low value of θ, Ē * is not feasible but Ē1 is stable, and as θ crosses its critical value, Ē * becomes stable and the stability of Ē * depends on the combined value of θ and β, as higher values can make the stable system unstable via Hopf bifurcation.In Figure 8b, a direct relationship between the parameters θ and m and the stability of the system is shown, where with an increased value of θ, only a very high value of m can make the system stable.
From the simulations, we observe that θ is one of the critical parameters for the model's system (1), and the stability of different equilibria depends on it.For very low values of θ, Ē * is not feasible, but Ē1 is stable, and as θ crosses its critical value, Ē * becomes stable via transcritical bifurcation, and finally, further higher values can make Ē * unstable via Hopf bifurcation (Theorem 2).

Numerical Simulation of System (14)
Simulations are presented for the dynamics of the system ( 14) in Figures 9-13 based on the analytical analysis presented in Section 3.
In Figure 9, we plot the transcritical forward bifurcation of Ẽ1 and the existence of Ẽ * with variations in the parameter α (Theorem 5, (ii)).A Hopf bifurcation diagram is plotted in Figure 10 ,taking α as the main parameter (using (Theorem 5 (iii)).The limit cycle disappears as the parameter α crosses the critical value α * .1.The dotted curve represents the maximum and minimum values of the periodic solutions and the solid blue line is the stable coexisting equilibrium Ẽ * .
In Figure 12, we show m = 0.3, α = 0.001 (i.e., there is very little competition among the prey individuals for finding a refuge habitat), and the system is stable.Now, if we increase the value of α 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 α has a stabilizing effect on the system.In Figure 13, we plot the stable regions of the coexisting equilibrium in the α-β and α-m parameter planes.The area of the stable region increases as both m and α are increased.But, β can make the stable system unstable via Hopf bifurcation (Theorem 5 (iii)).The critical value of β i.e., β * depends on the other model parameters, such as α.Here, we show that the critical value of β significantly depends on α.In (b), we use a higher value of β, and then we vary α and m.We observe that Ẽ * changes its stability level to become stable, except for at higher values of m and α.The stable condition is only shown in E 1 for higher values of m and α.From the numerical simulations, we observe that for model system ( 14), the stability of different equilibria depends on the parameter α.Here, for very low values of α, Ẽ * is unstable, and it becomes stable as α crosses the critical value via Hopf bifurcation, but for very high values of α, Ẽ * is not feasible, whereas Ẽ1 becomes stable via transcritical bifurcation (Theorem 5, (ii)).

Discussion and Conclusions
A prey-predator system including prey refuge has been a sought-after topic among ecological and mathematical modelers for almost the last five decades.Previous studies have assured that refuge has a stabilizing effect on prey-predator dynamics.Several models have already been established to represent various conditions and situations.Still, another finding is 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 cannot 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 intra-specific interactions must affect prey-predator dynamics.This study identified the impacts of cooperation and intra-specific competition among prey populations trying to escape predation through refuge on the stability of prey-predator dynamics.
The ecological motivation behind the formation of the mathematical models comes from the Hoogly-Matla estuarine mangrove ecosystem, where several islands are situated and are crisscrossed by numerous creeks that 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 such fish was described by Sarwardi et al. [34].In their study, Sarwardi et al. [34] described one important detritivorous fish from these creeks, namely Liza parsia, which lives in the creeks, and its predator Otolithoides pama visits the creeks during high tide only to capture food.The mangrove plants in this region extend from the supra-littoral zone to the lower-littoral 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.In many cases, due to an inadequate refuge habitat, all individuals of the prey population can not escape predation.In this situation, the prey population finds refuge either through cooperation or competition and avoids predation pressure.The present study focuses on the effects of these two important intra-specific interactions on the dynamics of prey-predator interaction with the prey refuge.
In a basic predator-prey model, the population dynamics of the predators and prey are influenced by factors such as predation, reproduction, and competition for resources.However, if we introduce cooperation and competition among the prey species themselves, the model becomes more complex.For the two situations mentioned above, in this paper, we considered a prey-predator system incorporating a prey refuge and assumed that the predator response function ia Holling type II.Cooperation (model ( 1)) and intra-specific competition (model ( 14)) among the prey population for escaping predation pressure through refuge are incorporated into the system to find out the process (either cooperation or intra-specific competition) that leads the population towards better stabilization.
Firstly, the stabilizing effect of refuge on the system dynamics was shown (Figure 4).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.Over-exploitation of the prey population by predators moves 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 about 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 the maximum number of individuals can escape predator easily, which is favorable for them but on the other hand, predators 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 the equilibrium state faster when the refuge facility is achieved through intra-specific competition among the prey population (Figure 12).Refuge during intra-specific competition has a stabilizing effect on the system dynamics.As intra-specific competition has a negative impact on individuals in the prey population, the refuge class is constantly controlled and never becomes too large to cause prey scarcity for a predator.A controlled amount of the prey population will always escape predators through refuge, which improves the survival of the prey population, and the predator population will not become extinct in any situation, as it never faces a food crisis.Due to this balanced interaction, the system dynamics reach an equilibrium faster.
In reality, the decision-making process to escape predation pressure by refuge either through competition or cooperation within a prey population is complex and does not follow any hard and fast rules.Cooperation or competition in various prey organisms depends on the spatial heterogeneity of the habitat, as when the shelter place (refuge size) is large enough to protect all prey individuals, only cooperation acts as the deciding factor, but any scarcity in the 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 by living in a group and cooperating conspecifically.Group living in a 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 that are attuned to predation threats [39,40].So, cooperation plays a major role in attaining refuge within a living prey population group, but this is not necessarily always the case.Underlying intraspecific competition may be found among the group members, which can be explained by the selfish-herd effect [41,42].Though each member enjoys the benefit of decreasing the domain of danger when living in a group than as a solitary forager, those individuals who are present outside the group still have a larger domain of danger compared to those at the core, and thus, they may compete with each other to take shelter at the core of the group [39,43].This study fills the research gap by dealing with the impacts of cooperation and intra-specific competition among the prey population for escaping from predators through refuge on the prey-predator system, and it was found that when the prey individuals become competitive with each other for shelter in the refuge habitat, the system dynamics achieve stability faster.On the other hand, any cooperation among the prey population regarding refuge turns the stable equilibrium into a limit cycle condition.
Interactions in the living world are very complex, and it is difficult to conclude whether competition or cooperation is helpful for refuges.A high rate of parental care has been observed in many bird species, including the protection of offspring from predators 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 successful refuge.In reality, when prey populations are small, they need cooperation rather than competition, but the reverse situation is noticed in case of large prey populations.Therefore, both cooperation and competition are inevitable processes during the refuge of the prey population, and the selection of either cooperation or competition depends on the conditions within the prey population.

Figure 1 .
Figure 1.Schematic diagram of the proposed model (1) for the Lotka-Volterra type prey-predator system, including the proportional prey refuge and cooperation among the prey population trying to escape predation pressure through refuge.

Figure 2 .Remark 2 .
Figure 2. Schematic diagram of the proposed model (14) for the Lotka-Volterra type prey-predator system including the proportional prey refuge and intra-specific competition among the prey population trying to escape predation pressure through refuge.Remark 2. From a mathematical point of view, model (1) and model(14) differ only in one term: −αx 2 in place of θx 2 .Thus, the analytical results (dynamics) of model (14) can be derived from the analysis of model (1).

Figure 4 .Figure 5 .Figure 6 .
Figure 4.The effect of refuge on the system's populations is shown: (a) prey population, (b) predator population, and (c) phase-portrait in x-y plane.The dashed line represents m = 0 (i.e., without refuge), and a solid line is used for m = 0.3.Other parameters are the same as those in Figure 3.

Figure 7 .Figure 8 .
Figure 7. (a,b): Bifurcation diagram is plotted taking the rate of cooperation (i.e., θ) as the bifurcation parameter.The set of parameters is β = 0.18, m = 0.3.The dotted curve represents the maximum and minimum values of the periodic solutions, and the solid black line is the stable coexisting equilibrium Ē * .

Figure 9 .Figure 10 .
Figure 9. Transcritical bifurcation (forward) is plotted with variations in the parameter α.The set of parameter values is taken as β = 0.3, m = 0.1, and the other parameters are taken from Table1.

Figure 11 .Figure 12 .
Figure 11.Phase trajectories are plotted in the x-y plane by taking three different initial conditions.The parameter values are same as those presented in Figure 5.

Figure 13 .
Figure 13.The region of stability of the equilibrium points of system (14) is shown in (a) α-β, (b) α-m parameter planes.In region A, Ẽ * is unstable; in region B, Ẽ * is stable; and in region C, Ẽ * is not feasible, but Ẽ1 is stable.In (a), m = 0.1 and in (b), β = 0.3, the other parameters are as given in Table1and Figure4.

Figure 11
Figure 11 also plots two different initial values of the model variables.This figure shows that the bifurcating periodic solutions are stable and the Hopf bifurcation is of the supercritical type (Theorem 6).From the numerical simulations, we observe that for model system (14), the stability of different equilibria depends on the parameter α.Here, for very low values of α, Ẽ * is unstable, and it becomes stable as α crosses the critical value via Hopf bifurcation, but for very high values of α, Ẽ * is not feasible, whereas Ẽ1 becomes stable via transcritical bifurcation (Theorem 5, (ii)).

Author Contributions:
Conceptualization, S.P. and S.R.; methodology, S.P. and F.A.B.; software, S.P. and F.A.B.; validation, S.P. and F.A.B.; formal analysis, S.P. and F.A.B.; investigation, S.P. and F.A.B.; writing-original draft preparation, S.P., S.R. and F.A.B. writing-review and editing, S.P. and F.A.B.; visualization, S.P. and F.A.B.; supervision, F.A.B.All authors have read and agreed to the published version of the manuscript.All authors contributed equally to this work.Funding: This research received no external funding.

Table 1 .
List of parameters used in numerical simulations and their references.
The existence and stability of the coexisting equilibrium are also shown.The set of parameter values is c = 0.75, m = 0.3, and the rest of the parameters are taken from Table1.