Abstract
In this work, our aim is to investigate the impact of a non-Kolmogorov predator-prey-subsidy model incorporating nonlinear prey refuge and the effect of fear with Holling type II functional response. The model arises from the study of a biological system involving arctic foxes (predator), lemmings (prey), and seal carcasses (subsidy). The positivity and asymptotically uniform boundedness of the solutions of the system have been derived. Analytically, we have studied the criteria for the feasibility and stability of different equilibrium points. In addition, we have derived sufficient conditions for the existence of local bifurcations of codimension 1 (transcritical and Hopf bifurcation). It is also observed that there is some time lag between the time of perceiving predator signals through vocal cues and the reduction of prey’s birth rate. So, we have analyzed the dynamical behaviour of the delayed predator-prey-subsidy model. Numerical computations have been performed using MATLAB to validate all the analytical findings. Numerically, it has been observed that the predator, prey, and subsidy can always exist at a nonzero subsidy input rate. But, at a high subsidy input rate, the prey population cannot persist and the predator population has a huge growth due to the availability of food sources.
1. Introduction
In the ecological system, the predator-prey interaction is one of the most significant tools which is comparatively easy to observe in the field. But fear of the predator felt by the prey (indirect effect) also plays a vital role since its effect is stronger than direct predation [1, 2]. The cost of fear can reduce the reproduction rate of prey because it affects the physiological condition of prey population. As a result, the prey species may get a long-term loss. In support of this, it is mentioned that, in the Greater Yellowstone Ecosystem, wolves (Canis lupus) affect the reproductive physiology of elk (Cervus elaphus) [3]. When the prey species recognize the predator signal (chemical/vocal), they spend more time to become keenly watchful to detect danger rather than in foraging. So, the birth rate of the scared prey reduces and adopts some survival mechanisms like starvation [1, 2]. For examples, some birds react to the sound of predator with antipredator defenses [1, 2] and they flee from their nests at the first sign of danger [2]. This antipredator behaviour may affect survival and reproduction of the birds [2]. It has been experimentally investigated that, in the absence of direct killing, the reproduction of the offspring of song sparrows (Melospiza melodia) could be reduced by as a result of impact of feeling fear created by the predator [4]. So, this reduction caused by the antipredator behaviour affects the birth rate and survival of offspring. Thus, the cost of fear (apart from direct predation) should be introduced in a predator-prey interaction. Mathematical formulation of the impact of fear on the two species prey-predator system has been initiated by Wang et al. [5] in 2016 introducing fear factor: . It involves a parameter denoting the level of fear to represent the antipredator behaviour of the prey. Some research works have already been done on the ecological system under the influence of predation fear felt by prey species [6–17]. Moreover, the impact of fear in a two-species predator-prey model with prey refuge was analyzed by many researchers [11, 18, 19].
In evolutionary biology, prey refuge is a concept which helps an organism to protect themselves from predation by hiding in an area inaccessible to the predator, for example, in a wolf-ungulate system, ungulates may seek refuge by migrating to areas outside the core territories of wolves. Also, it has many significant roles on the dynamics of predator-prey interactions: prey refuge may decrease the chance of extinction of prey. Researchers have mainly used the dynamic nature of predator-prey model with linear prey refuge (that is, amounts of prey are unavailable to the predator, where is the coefficient of refuge and is the biomass of prey species) with Volterra response [20–22]. Recently, Mondal and Samanta have studied the dynamics of the predator-prey system with prey refuge dependent on both species (that is, amounts of prey are free from the predator risk, where , is prey refuge coefficient, is the biomass of prey population, and is the biomass of predator population) in the presence of additional food (for details, see [23]). In 2020, Mondal and Samanta [11] have also analyzed the dynamics of predator-prey interaction having nonlinear prey refuge function which is the amount of prey that are free from predation, where is half saturation constant and is the biomass of predator population.
Many experimental studies suspect that the introduction of resource subsidies may disrupt otherwise stable food web linkages [24–26]. Such concept is significant for resource management purposes. It is learnt that reintroduced wolves in Yellowstone Park switch to bison when their preferred ungulate prey, namely, elk, are rare in the concerned ecosystem [27]. Mathematically, the influence of resource subsidy on the predator-prey model has been initiated in the work of Nevai and Van Gorder [28]. They have discussed how different subsidy input rate may affect the prey and predator population to persist in the ecosystem.
Generalist predator can consume more than one food source: either multiple prey population or a combination of prey population and resource subsidy. There are many rich theoretical research on ecological systems involving generalist predator [29–31]. Also, there are a variety of real-life applications for such systems [32–34]. From literature surveys, it has been shown that generalist predator can persist in an ecosystem even if one particular prey species is going towards extinction [29–31].
In 2012, Nevai and Van Gorder [28] first extended the Kolmogorov model to a non-Kolmogorov predator-prey-subsidy model. It has been observed that the predator-prey-subsidy model occurs in the arctic foxes (predator), lemmings (prey), and seal carcasses (subsidy). Motivated by the works of Das and Samanta [35], Nevai and Van Gorder [28], and Xu et al. [36], we have analyzed the dynamical behaviour of a mathematical model of non-Kolmogorov form that includes the three components (predator, prey, and subsidy) with the impacts of nonlinear prey refuge function and the fear effect felt by the prey in the presence of the predator. To the best of our knowledge, there does not exist any mathematical model to explore the impact of fear effect incorporating nonlinear prey refuge function in predator-prey-subsidy interaction.
The organization of this work is structured as follows: in Section 2, a mathematical model has been formulated with the influences of nonlinear prey refuge and fear effect. Section 3 shows that the proposed model is well-behaved. In Section 4, feasibility criteria and stability of all the equilibria of the proposed system (in absence of delay) have been studied. The equilibria can change their stability nature through transcritical and Hopf bifurcation which are also analyzed in this section. Generally, the reduction of prey’s birth rate due to the effect of fear will not be an instantaneous biological process but deviated through some time lag, so the study of time-delay is very meaningful to obtain the more realistic dynamics. So, Section 5 deals with the dynamic behaviour of the delayed system for two equilibrium points (subsidy free) and (interior), respectively. Section 6 provides the numerical computations which support the analytical calculations. Section 7 provides a brief conclusion about the system dynamics.
2. Model Formulation
In 2020, Mondal and Samanta [11] analyzed the dynamics of a delayed predator-prey interaction incorporating nonlinear prey refuge function under the influence of fear effect and additional food. Motivated by the work of Mondal and Samanta [11], we have first considered a predator-prey-subsidy model with nonlinear prey refuge function where the prey and subsidy occur in the same habitat and they are both consumed by a single generalist predator according to the following differential equations:where is prey population, denotes the population of the subsidy, and is the generalist predator which exploits both the prey and subsidy. For example, wolves (predator) consume both deer (prey) and salmon carcasses (subsidy) [37].
The term represents the quantity of prey available to the predator, i.e., amounts of prey are free from predation risk where is designated as nonlinear prey refuge function. Also, we have modeled the dynamics of a generalist predator with Holling type II [38–41] response function in the presence of nonlinear prey refuge function.
All parameters are positive (except ) and biologically meaningful. Parameters are described in Table 1.
Apart from direct consumption, feeling of fear among the individuals of the prey species in presence of predator is very common in predator-prey interaction which changes life-history, behavioural responses, and reproduction capability of prey species. In ecology, effect of fear is a common factor, but there does not exist any considerable attention to introduce the impact of fear in the mathematical modeling. Experimental studies indicate that the feeling of fear among the individuals of the prey species in presence of predator reduces the prey’s birth rate. So, birth rate of prey species is multiplied by a monotone decreasing function , where is a level of fear [5]. The fear function satisfies the following conditions:(1): when there is no fear effect on the prey species, the birth rate of the prey is not reduced(2): when there is no predator, the birth rate of the prey species is not reduced in the presence of fear effect(3): when fear effect increases, the birth rate of the prey reduces(4): when predator species increases, prey population reduces
Our main focus is to analyze the dynamic nature of the predator-prey-subsidy model with the influence of nonlinear prey refuge and fear effect. So, system (1) can be modified in the following aspects:with initial conditions:
Throughout the analysis of this work, we have taken which is biologically meaningful.
3. Positivity and Uniform Boundedness
Theorem 1. Every solution of system (2) with (3) uniquely exists and is positive for all .
Proof. Solution () of (2) with (3) exists and is unique on , where [42].
From (2) with (3),Now, we claim that for all . If it does not hold then there exists such that , , and on . From the second equation of (2),a contradiction with . So, .
Hence, solutions of (2) stay positive for all .
Theorem 2. All solutions of system (2) which start in are asymptotically uniformly bounded.
Proof. Case 1: if , from the first equation of (2), Let us take . Differentiating both sides with respect to , we obtain Let Then, Using the Gronwall inequality, we obtain Thus, all solutions of system (2) enter into the region: Case 2: if , from the first equation of (2) we obtain .Now, from the second and third equations of (2), we haveFor large ,LetThen,Using Gronwall inequality, we obtainHence, the theorem.
4. Equilibrium Points and Stability Analysis
4.1. Equilibria
4.1.1. Trivial Equilibrium Point
Extinction: .
4.1.2. Axial Equilibrium Points
(i)Subsidy only: (ii)Prey only: exists if and
4.1.3. Planer Equilibrium Points
(i)Subsidy free: exists if , , and where and can be obtained by solving the equations: and we get where is a positive root of the equation: Here, and exists if and .(ii)Prey free: exists if and .(iii)Predator free: exists if .
4.1.4. Interior (Coexistence) Equilibrium Point
Solving the following system of equations,we can obtain using the software MATHEMATICA with the following existence conditions:(1)(2)(3) (otherwise, predator population goes into extinction)
4.2. Stability Analysis
Now, we will study the stability conditions of all equilibria for the proposed system (2).
The Jacobian matrix at is given by
The eigenvalues of are , , and . Then. we have stated the following theorem.
Theorem 3. Trivial equilibrium point is locally asymptotically stable (LAS) if and unstable if .
The Jacobian matrix at is as follows:
We see thatare the eigenvalues of the matrix . Thus, we have the following theorem.
Theorem 4. Axial equilibrium point (subsidy only) is locally asymptotically stable if and and unstable if either or or .
The Jacobian matrix at is given by
The eigenvalues of are , or . Thus. we have stated the following theorem.
Theorem 5. Axial equilibrium point (prey only) is locally asymptotically stable ifIf with , then is unstable.
The Jacobian matrix at is as follows:where
The characteristic equation corresponding to is expressed aswhere , , and .
Theorem 6. Subsidy-free equilibrium point is locally asymptotically stable if , , and .
The Jacobian matrix corresponding to is given bywhere
The characteristic equation corresponding to is expressed aswhere
Theorem 7. Prey-free equilibrium point is locally asymptotically stable if and .
The Jacobian matrix corresponding to is given bywhere
The characteristic equation corresponding to is expressed aswhere
Theorem 8. Predator free equilibrium point is locally asymptotically stable ifprovided and .
The Jacobian matrix corresponding to is as follows:where
The characteristic equation corresponding to is expressed aswhere
Theorem 9. The coexistence equilibrium is locally asymptotically stable if , , and , where are stated in (41).
4.3. Local Bifurcations of Codimension 1
4.3.1. Transcritical Bifurcation
Theorem 10. System (2) undergoes a transcritical bifurcation around if and ( stands for transcritical bifurcation).
Proof. We apply Sotomayor’s theorem [43] to prove the occurrence of a transcritical bifurcation around with as bifurcation parameter. For applicability of Sotomayor’s theorem, exactly one of the eigenvalues of the Jacobian matrix at must be zero and other eigenvalues must have negative real parts. So, we need to fulfill the condition .
The eigenvectors of and corresponding to the zero eigenvalue of are obtained as and , respectively, where , , , , , , , , , , .
Compute , , and as follows:where and , , and are given byTherefore, by Sotomayor’s theorem [43], system (2) undergoes a transcritical bifurcation at around the axial equilibrium point .
Theorem 11. System (2) exhibits a transcritical bifurcation around if
Proof. Let us apply Sotomayor’s theorem [43] to prove the occurrence of a transcritical bifurcation around with as bifurcation parameter. For applicability of Sotomayor’s theorem, exactly one of the eigenvalues of the Jacobian matrix at must be zero and other eigenvalues must have negative real parts.
The eigenvectors of and corresponding to the zero eigenvalue of are obtained as and , respectively, where , , , , , , , , , , , .
Compute , , and as follows:where and , , and are given byTherefore, by Sotomayor’s theorem [43], system (2) exhibits a transcritical bifurcation at around the axial equilibrium point .
Theorem 12. System (2) undergoes a transcritical bifurcation around if
Proof. Proof is the same as in Theorem 10.
Theorem 13. System (2) undergoes a transcritical bifurcation around if
Proof. Proof is the same as in Theorem 11.
4.3.2. Hopf Bifurcation around
Let us consider as a bifurcation parameter of system (2) where the characteristic equation at is
Then, Hopf bifurcation theorem is stated as follows.
Theorem 14. (Hopf bifurcation theorem [44]). If , , and are the smooth functions of in , , for which the characteristic equation (50) has the following:(i)A pair of imaginary eigenvalues with and so that they become purely complex at and (ii)The other eigenvalue is negative at ; then, a Hopf bifurcation appears around at
Theorem 15. If is locally asymptotically stable, then a Hopf bifurcation is exhibited around when passes through its critical value provided , , and ( is a positive root of equation ).
Proof. At , we can write equation (50) asThe roots of equation (51) are , , and . Also , , and are the smooth functions of . So, the roots of equation (59) have the form , , and where are real functions of in an open neighborhood of for . Next, we verify the transversality condition:Putting in (59), we getDifferentiating both sides with respect to , we haveComparing real and imaginary parts from both sides, we obtainwhereand .
Multiplying (55) by and (56) by and then adding, we getAt , Case 1: , . Then, , , , and . . Case 2: , . Then, , , , and . .Also, .
Hence, this theorem is proved by virtue of Theorem 14.
4.3.3. Hopf Bifurcation around
Let us consider as a bifurcation parameter of system (2) where the characteristic equation of isand then Hopf bifurcation theorem is stated as follows.
Theorem 16. (Hopf bifurcation theorem [44]). If , , and are the smooth functions of in , , for which the characteristic equation (59) has the following:(i)A pair of imaginary eigenvalues with and so that they become purely imaginary at and (ii)The other eigenvalue is negative at ; then a Hopf bifurcation occurs around at
Theorem 17. If is locally asymptotically stable, then a Hopf bifurcation appears around subsidy-free equilibrium when passes through its critical value provided , , and ( is a positive root of equation ).
Proof. Proof is the same as in Theorem 15.
5. Delayed Dynamical System
In biological point of view, many processes, both natural and man-made, include time-delay. The study of delay factor makes our system much more realistic than non-delayed system. Also, a delay differential equation reveals much more complicated dynamics than an ordinary differential equation (for details, see [10–13, 23, 45–49]).
In reality, after sensing the vocal cue, individuals of prey species take some time for assessing the predation risk. So, the effect of fear (felt by prey) of predator does not respond spontaneously on the birth rate of prey population; some time lag must be needed. In view of this fact, the predator-prey-subsidy interactions (2) can be modified as follows:
The initial conditions are assumed as ()
Let us linearize (60) using the following transformations:
It leads towhere ,
The characteristic equation corresponding to (62) iswhere
If , of system (60) is LAS provided equation (64) has no purely imaginary roots and it is also LAS for . Further, it has been shown that stability nature of switches at . Already, it has been derived that is LAS provided , , and for (non-delayed system). Let us discuss if the real part of the roots of equation (64) gradually increases to reach zero and eventually turns to a positive value when increases.
Substituting in equation (64), we have
Equating respective real and complex parts from both sides, we obtain
Now, let us examine whether equation (64) has purely imaginary roots or not. For this purpose, let us take . Then, equations (67) and (68) become
Eliminating from (69) and (70) (squaring and adding), we get
Putting , we have
This is a cubic equation of . It is noticed that . So, equation (72) has exactly one positive real root if , i.e., if .
Let be a positive root of (77); then, .
Lemma 1. (see [50]). Consider the exponential polynomial:where and are constants. As vary, the sum of the orders of zero of in the open half plane can change only if a zero appears on or crosses the imaginary axis.
Now, let us discuss the existence of Hopf bifurcation around with as a bifurcation parameter.
Theorem 18. Suppose exists and is locally asymptotically stable for system (2) when . If , then there exists a critical value such that of system (60) is LAS when and unstable when , whereand (minimum value). Also, system (60) exhibits Hopf bifurcation around at provided , where
Proof. If , then (72) has exactly one positive root , i.e., from (69) and (70), , are obtained as functions of :If is locally asymptotically stable, the stability behaviour of will remain unaltered for (using Butler’s Lemma [51]).
To check the transversality condition, , let us differentiate (67) and (68) with respect to and set and . The following equations are obtained:whereSolving (77) and (78),Now, we have , if . Hence, the transversality condition is satisfied and a Hopf bifurcation occurs around when passes through its critical value .
Now, linearize system (60) using the transformations , , and :where ,
The characteristic equation corresponding to (81) iswhere
If , of system (60) is LAS provided equation (83) has no purely imaginary roots and it is LAS for . Furthermore, it has to be noted that changes of stability occur at . Already, it has been discussed that is LAS when provided , , and . Here, equation (83) is a transcendental equation, so it contains infinitely many eigenvalues. In this situation, we cannot apply the Routh–Hurwitz criteria to determine the stability of system (60). To understand the stability behaviour, our necessity is to check the sign of the real parts of the eigenvalues of equation (83).
Now, putting in equation (83), we have
Equating respective real and complex parts from both sides, we get
To check whether (83) has purely imaginary roots or not, set ; then, (86) and (87) become
Eliminating from (88) and (89) (squaring and adding), we get
Putting , we have
This is a cubic equation of . It is noted that . So, equation (91) has exactly one positive root if , i.e., if .
Let be a positive root of (91); then, .
Let us study the existence of Hopf bifurcation around with as bifurcation parameter.
Theorem 19. Suppose exists and is locally asymptotically stable for system (2) when . If , then there exists a critical value such that of system (60) is LAS when and unstable when , whereand (minimum value). Also, a supercritical Hopf bifurcation is exhibited around at provided , where
Proof. Proof is similar to that in Theorem 18.
Theorem 20. Suppose interior (coexistence) equilibrium point exists and is locally asymptotically stable for system (2) when . Let equation (91) have exactly two positive roots when and irrespective of sign of . Moreover, letthen there is a positive integer such thatand there are switches from stability to instability to stability; that is, whenthen is locally asymptotically stable and whenthen is unstable. Further, at , system (60) experiences Hopf bifurcation providedwhere
Proof. Proof is similar to that in Theorem 18.
6. Numerical Computations
Here, we have illustrated numerical simulations to verify the analytical findings of the proposed system (2). We select a parameter set: {, , , , , , , , , , , , , }. Under this set of parametric values, the stable nature of is shown in Figure 1. If we take subsidy input rate and other parametric values are chosen from the data set of Figure 1, then the subsidy only equilibrium exists and stable nature of with time is depicted in Figure 2. Now, we choose another parameter set: {, , , , , , , , , , , , , }. Under this set of parametric values, the prey only equilibrium exists and stable behaviour of is presented in Figure 3. Let us choose the parameters as follows:
If we take subsidy input rate and other parameters are taken from set (100), then subsidy-free equilibrium point exists and is locally asymptotically stable. Stable time series and stable phase diagram are represented in Figure 4. In the same manner, if we change the value of the parameter and others are the same as in the data set of Figure 4, then it is observed that is unstable accompanied with a limit cycle (see Figure 5). From Figures 4 and 5, it can be easily noted that there must exist a threshold value of , say for which unstable behaviour of changes to stable spiral. Since the vector fields for and are qualitatively different, a Hopf bifurcation is created around taking as bifurcation parameter (see Figure 6). For the set of parameter values {, , , , , , , , , , , , , }, prey free equilibrium point exists and is stable (see Figure 7). Next, let us take a different set of parameters of system (2): {, , , , , , , , , , , , , }. Then, predator free equilibrium point is locally asymptotically stable. The stable behaviour with time is shown in Figure 8.
(a)
(b)
(a)
(b)
(a)
(b)
If we take subsidy input rate and others are fixed as in the data set of Figure 4, then exists and is locally asymptotically stable. Figure 9 depicts the stable behaviour of . Comparing Figures 4 and 9, it is observed that subsidy input rate enhances the value of component of and decreases the value of component of . Also, from Figure 10, it is noticed that the prey population is leading towards extinction and the predator population has enormous growth (due to huge supply of food source) at high subsidy input rate (when ) in the presence of fear felt by prey population. So, it can be concluded that it is not possible to control prey population from extinction in presence of nonlinear prey refuge because they cannot get enough time to protect themselves from predation risk. After extinction of prey species, predator can easily survive with the help of resource subsidy. Thus, the parameter has great importance in the proposed population dynamics.
(a)
(b)
Moreover, Figure 11 represents the unstable nature of when and other parametric values are the same as in Figure 9. So, the parameter has an interesting nature because there exists a threshold value of for which unstable nature (limit cycle) of switches to stable behaviour (stable spiral) when passes through its critical value ; i.e., the vector fields for and are topologically different. Hence, a Hopf bifurcation occurs around and Figure 12 depicts the corresponding bifurcation diagram taking as bifurcation parameter. Also, it has to be noted from Figure 13 that, in the absence of fear effect, the oscillatory behaviour of changes to stable nature when subsidy input rate crosses its critical value (approximately) and since predator population has huge growth rate at very large value of subsidy input rate , the prey population cannot persist in ecosystem in presence nonlinear prey refuge. This phenomenon is very interesting because the prey refuge cannot control the prey population from extinction due to enormous growth of predator when subsidy input rate is very high. In this manner, system (2) is not persistent.
(a)
(b)
(c)
(a)
(b)
(c)
(a)
(b)
(c)
Further, Figure 14 depicts that the nature of steady states and when . Here, the predator population cannot go extinct for large value of coefficient of prey refuge parameter. Also, Figure 15 shows the changes of predator’s growth at the steady states and for three different fear levels when varies from 0 to 1. Here, also the predator can persist for large . From here, it may be concluded that system (2) is always persistent for small subsidy input rate in the presence of nonlinear prey refuge function. Again, Figures 16 and 17, respectively, show that, in the absence of fear effect , the equilibria and are approaching towards stable state by excluding the existence of oscillatory behaviour taking as the bifurcating parameter. In this manner, predator population also survives in ecosystem for large coefficient of prey refuge parameter .
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
(c)
A “transcritical bifurcation” (BP) occurs at around . At this point, exactly one eigenvalue of the Jacobian matrix is zero and others have negative real parts. Bifurcation diagram 32 depicts that when , then is unstable, and when , then is stable. Also, Figure 18 exhibits that when , then is unstable, and when , then is stable. So, a transcritical bifurcation is exhibited at around . Similarly, Figures 19 and 20 depict the transcritical bifurcation diagrams around equilibrium and taking and as bifurcation parameter, respectively.
6.1. Effect of Time-Delay on Population Dynamics
Now, let us perform the numerical computations to validate the analytical results of the delayed model (60). For the parameter set {, , , , , , , , , , , , }, equation (72) has exactly one positive root 0.2246 (correct up to four decimal places). So, from Theorem 18, the planer equilibrium point is stable when and unstable when . Stable time series and stable phase trajectory of are shown in Figure 21 when . Also, Figure 22 depicts the corresponding unstable behaviour of when . Moreover, Figure 23 presents the supercritical Hopf bifurcation diagram around taking as bifurcation parameter.
(a)
(b)
(a)
(b)
(c)
(a)
(b)
For the parameter set {, , , , , , , , , , , , }, equation (91) has exactly two positive roots, 0.1127 and 0.0526 (correct up to four decimal places). Then from Theorem 20, we have calculated , , , , and . The interior equilibrium is locally asymptotically stable when and unstable when and . At and , Hopf bifurcation appears around . Figures 24–26 depict the stable nature of for respectively. Also, unstable behaviour of is presented in Figures 27–29 for and , respectively. The corresponding bifurcation diagrams are depicted in Figures 30 and 31 .
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(a)
(b)
(c)
(a)
(b)
(c)
7. Conclusion
We have analyzed a system for generalist predator which utilizes more than one food source: predator-prey-subsidy model of non-Kolmogorov form introducing nonlinear prey refuge function and the effect of fear felt by prey population. Our main interest is to find the situations such that dynamical stability and instability appear so as to make out more fully how subsidy may influence the predator and their prey. It has been shown that the solutions of system (2) remain positively invariant always and they are asymptotically uniformly bounded. These, in turn, imply that system (2) is biologically well-behaved. Existence criteria and stable behaviour of all the biologically meaningful equilibria have been discussed. It has to be noted that Hopf bifurcations are exhibited around (subsidy free) and (interior) of system (2) considering as a bifurcating parameter (see Figures 4–6, 9, 11, and 12). Also, observing Figures 6 and 12, it can be concluded that high levels of fear can stabilize system (2) by excluding the existence of periodic solutions. These phenomena are biologically significant because prey species are aware after a certain level of fear; i.e., after a certain level of fear, they are not affected as they are aware and show signs of habituation.
Moreover, this work derives transcritical bifurcations (local bifurcation of codimension 1) at the various equilibrium points , , , and , respectively (see Figures 18–20 and 32).
Also, we have discussed numerically the influences of coefficient of prey refuge parameter on the nature of the equilibrium points (zero subsidy input rate) and (fixed small subsidy input rate) irrespective of fear level . Noting Figures 14–17, it is observed that both the prey and predator species always persist in ecosystem due to continuous increment of coefficient of prey refuge. But Figures 10 and 13 depict that, irrespective of fear level, a highly subsidized predator should indeed drive the prey population towards extinction regardless of whether the prey and subsidy arise in the same habitat. This phenomenon is ecologically meaningful because the prey population cannot get enough time to protect themselves from predation risk for enormous growth of predator at high subsidy input rate. So, the prey population is leading towards extinction, but the predator species can easily survive in ecosystem with the help of resource subsidy. Thus, the study of system (2) is ecologically very significant.
In reality, fear effect does not instantaneously reduce the birth of a prey population, but some time lag should be needed to create an impact on the birth rate of the prey population. We have considered that there is a time-delay on the impact of fear to the birth rate of prey, from the instance it perceives the fear of predator through any means. So, the incorporation of time-delay makes system (60) more realistic. It is noted that delay parameter has a significant role because there exists a threshold value such that stable behaviour of planer equilibrium point (in the absence of subsidy input rate) switches to oscillatory nature when passes through its threshold value ; i.e., the vector fields for and are qualitatively different. So, system (60) exhibits a supercritical Hopf bifurcation around considering as bifurcation parameter (see Figures 21–23). Also, a rigorous study of the stability and bifurcation of interior equilibrium point has been performed. Our analysis describes that the delay within a certain specified range could maintain the stable behaviour of . On the other hand, the delay could drive the system into an unstable state. Hence, the study of the time-delay parameter has a regulatory impact on the whole system.
Data Availability
The data used to support the findings of the study are available within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors are grateful to Dr. Eberhard O. Voit (editor) for his careful reading, valuable comments, and helpful suggestions, which have helped them to improve the presentation of this work significantly. This research was completed during the visit of Prof. G.P. Samanta to the Institute of Mathematics of the University of Santiago de Compostela. This work has been partially supported by the Agencia Estatal de Investigacion (AEI) of Spain, cofinanced by the European Fund for Regional Development (FEDER) corresponding to the 2014–2020 multiyear financial framework project MTM2016-75140-P and by Xunta de Galicia under grant ED431C 2019/02.