A systematic study of autonomous and nonautonomous predator-prey models for the combined effects of fear, refuge, cooperation and harvesting.

In the present study, we investigate the roles of fear, refuge and hunting cooperation on the dynamics of a predator-prey system, where the predator population is subject to harvesting at a nonlinear rate. We also focus on the effects of seasonal forcing by letting some of the model parameters to vary with time. We rigorously analyze the autonomous and nonautonomous models mathematically as well as numerically. Our simulation results show that the birth rate of prey and the fear of predators causing decline in it, and harvesting of predators first destabilize and then stabilize the system around the coexistence of prey and predator; if the birth rate of prey is very low, both prey and predator populations extinct from the ecosystem, and for a range of this parameter, only the prey population survive. The fear of predators responsible for increase in the intraspecific competition among the prey species and the refuge behavior of prey have tendency to stabilize the system, whereas the cooperative behavior of predators during the hunting time destroys stability in the ecosystem. Numerical investigations of the seasonally forced model showcase the appearances of periodic solution, higher periodic solutions, bursting patterns and chaotic dynamics.

Introduction

One of the major driving forces behind evolution of species in the ecosystem is the interaction among its individuals. The interaction between an organism and its natural enemies is termed as the predator–prey interaction that determines the mortality of prey and the birth of a new predator, and plays a major role in the movement of energy through food chains. These interactions are one of the central themes in mathematical ecology [1]. Differential equations have been extensively used to model real-life problems [2-4]. After the seminal works of [2, 3], a plethora of mathematical models on predator–prey interactions have been investigated with the inclusion of various types of functional responses which depict realistic situations [5-10]. It is worthy to note that these functional responses depict only the effect of the direct killing of prey by predator. However, the predators not only affect the ecology of prey directly by consuming them, but they also inflict their behavior and physiology indirectly [11, 12]. In order to avoid the fear of predation, the prey species opt to change their habitat to a safer place. This in turn increases the survivability of prey for a short period of time, but has negative influence in the long-run. For instance, the fear of wolves alters the rate of reproduction of Elk to a greater extent [13]. Indeed, manipulation of the fear factor plays a vital role in shaping the structure of ecological community [11, 14, 15]. In [11], authors have found that neither the high population density of Snowshoe hare nor food level is the main reason behind the decreasing reproduction rate; this happens only due to the high level of fear of predators. Fear of predation can reduce the reproduction of Song sparrows by 40% [14]. At first, Wang et al. [16] constructed a mathematical model by letting the fear factor to suppress the reproduction of prey population. Many studies have demonstrated that the predator–prey systems with fear effect exhibits comparatively a more complex dynamics [17-19]. Role of space in predator–prey models with fear effect is considered in [20, 21]. In [22], authors have investigated the impact of media induced fear on the recent COVID-19 pandemic situation.

Fear-induced refuge is a fascinating key factor acquired by prey population in an ecological community. It greatly affects the dynamics of a food web system and plays a crucial role in maintaining the ecological balance of prey and predator [23]. Firstly, the effect of refuge was demonstrated by Gause [24] in his experiment with the protozoan, Didinium (predator) and Paramecium (prey). The Snowshoe hare (prey) prefers the forested habitat with dense undercover, where they are comparatively safe from their predator (e.g., Canada lynx) [25]. Mathematical models considering the effect of prey refuge have been investigated by many researchers in the recent past [26-29]. Mondal et al. [30] have investigated the complex dynamics of a predator–prey system in the presence of resource subsidy by considering the influence of nonlinear prey refuge and fear effect. Recently, authors have assumed that the proportion of prey taking refuge is a function of predator density [31, 32], a more realistic consideration. In these studies, stabilizing role of prey refuge has been observed on the dynamics of predator–prey systems.

For many species, social communication is an integral part of the life history traits of a population and regulates the dynamics of system, and also affects the population densities [33]. For instance, the cooperation within a population is a well-known important phenomenon in ecological systems [34]. Nudds [35] observed a clear connection between wolf pack size and food acquired per wolf, by using sparse data from the literature. Smaller packs acquired substantially less food per wolf than those of optimal size. There are many other living organisms that cooperate among themselves for the hunting process. Males pursue less, though they refrain more than females; refraining while a group hunt is more common during hunts of prey that seem to be easier to capture [36]. During hunts of wart hog, lions are more likely to refrain Phacochoerus aethiopicus and less likely to refrain at the hunting time of zebra, Equus burchelli, and buffalo, Syncerus caffer. Aplomado falcons (Falco femoralis) sometimes hunt in pairs while chasing the birds; in eastern Mexico, 29% of 349 hunts were observed involving mated pairs of falcons simultaneously chasing the same prey animal; and 66% hunts of birds were tandem pursuits [37]. Some mathematical models have been investigated to explore the role of hunting cooperation in predator–prey systems [31, 38, 39]. In [38], it was observed that the hunting cooperation destabilizes a predator–prey system. Mondal et al. [40] explored complex dynamics such as bistability along with local and global bifurcations in a predator–prey system where predators are generalist and also cooperate during hunting. Recently, Roy et al. [41] have investigated the combined effects of fear, refuge and hunting cooperation in a seasonally forced eco-epidemic predator–prey system in which the predators have the ability to detect between susceptible and infected prey. They found that hunting cooperation can induce chaos in the ecological system.

Harvesting of the species is another major issue in a predator–prey system. Harvesting is very important for the conservation of the resources and also for the social management, on the ecological as well as economical grounds. The over-exploitation of an ecological system can be controlled by harvesting of its species. Several studies have explored the role of harvesting in ecological systems [42-44]. Clark [43] proposed the optimal management on harvesting in shery system for the first time. An optimal harvesting policy is explored by Zhang et al. [44], on a stage-structured two species model. Sahoo et al. [45] studied the role of alternative resources on the dynamics of a harvested predator–prey system. They showed the capability of suitable alternative resources to reduce the risk of extinction of predator species if the rate of harvesting is too high. Sk and Pal [46] have explored the dynamics of an eco-epidemic system with infection in prey species by considering that the predators are generalist. In this study, they have deeply investigated the effects of fear, refuge and harvesting. Further, they have seen the impact of stochasticity on the survival as well as extinction of species in the ecological system. In [47], authors have explored the dynamics of a predator–prey system with fear effect based state-dependent harvesting. Mondal et al. [48] explored different dynamical features of an imprecise predator–prey system with fear effect and nonlinear harvesting of predator in an uncertain environment.

As many environmental factors affecting the survival of species in the ecological community are seasonally forced, considerations of parameters in ecological models as periodic functions of time rather than constants would mimic comparatively more realistic situations [49, 50]. For instance, seasonal variation appears in the birth rate of many species that is caused due to humidity, temperature, rainfall, abundance of food, change in daylight period, etc., [51]. The seasonality in birth rate is also due to the fact that some seasons being evolutionary beneficial and favorable for giving birth and raising the young. It has been well documented that the level of fear varies due to change in predation pressure, intraspecific competition or distribution of resources among the individuals of prey species [52, 53]. The study of Sk et al. [54] has documented the appearances of higher periodic solutions and bursting patterns in a seasonally forced predator–prey system with fear, refuge and additional foods for predators. Periodic changes of refuge and hunting cooperation in a predator–prey system were considered by Sk et al. [31]. In [32], authors have investigated seasonal effects of fear, refuge and hunting cooperation in a predator–prey system. Tiwari et al. [55] have deeply investigated the effects of fear, migration and switching in an ecological system with two prey and one predator. They have studied the dynamics of the system in the absence as well as presence of the seasonal forcing in the costs of fear.

In most of the above-mentioned studies on predator–prey systems, the fear effect has been considered only on the birth rate of prey population; a little bit attention has been paid on the effect of fear of predator on the death of prey species [41, 46, 56]. Here, we investigate a predator–prey system where the birth rate of prey reduces due to fear of predator, while the fear of predators enhances the intraspecific competition among the individuals of prey species. The prey population are assumed to take refuge; the proportion of prey population taking refuge is a function of predator’s density. The predators are considered to cooperate each other at the time of hunting. We also take into account the harvesting of predators to reduce the predation pressure and maintain the bio-diversity. Later, we consider seasonal patterns of some of the ecologically important model parameters and investigate their impacts on the dynamics of the system. Thus, the aim of this study is two folds: firstly, we see how the fear effects, prey refuge, hunting cooperation and harvesting of predators regulate the dynamical interactions of prey and predator in an ecological community, and secondly, we investigate how the seasonal influences on model parameters affect the dynamics of ecological systems.

The remaining parts of this paper are organized in the following way. We introduce our autonomous system for the combined actions of fear, refuge, cooperation and harvesting in the next section and analyze the proposed model in the following section. In Sect. 3, we extend our autonomous model by letting seasonal variations of some of the model parameters. The nonautonomous system is analyzed in the following section. Numerical simulations are performed in Sect. 4, which explore rich dynamics of autonomous as well as nonautonomous systems. We wrap up the paper with conclusion and discussion in Sect. 5.

The autonomous model

Here, we consider an ecological community which consists of a single prey and a single predator populations. At any instant of time , we denote the densities of prey and predator populations by N and P, respectively, that are measured in terms of number per unit area. We construct our mathematical model for predator–prey system of an ecological community based upon the following ecological facts. In view of the above said assumptions, we propose the following mathematical model for the interacting dynamics of prey and predator populations in ecological system:Recently, Mondal and Samanta [59] explored the impacts of fear and its carry-over effect in a predator–prey system in which a constant proportion of prey population are able to take refuge. Our model system (1) is different from the model studied by Mondal and Samanta [59] in the following sense. Further, our model system (1) is different from the model proposed and studied by Guin et al. [60]. In [60], authors have investigated the dynamics of a predator–prey system in which the prey population are harvested at a linear rate, by considering ratio-dependent functional response. They have also assumed that the proportion of prey population taking refuge depends on the predator’s density. The effects of fear factor and the hunting cooperation of predators were not considered in [60]. In [60], the main objective was to investigate the effects of diffusion on the predator–prey system, and the generated spatial patterns.

The growth of prey is logistic with as its birth rate and as mortality due to intraspecific competition in the absence of predation and fear effect; the prey population dies naturally at the rate .

Anti-predator behavior of prey population due to the fear effect causes reduction in the reproduction rate of the prey species. Moreover, the fear of predators increases the intraspecies competition among the prey species [41, 46, 56]. Let and be the levels of fear of predators causing a decline in the reproduction rate of prey population and incline in the intraspecies competition among the individuals of prey species, respectively.

Refuge taken by prey to avoid the predation risk depends upon the frequency of encounter of prey by a predator and depends on the densities of both prey and predator population in the ecological system. So, it is reasonable to consider the number of prey taking refuge to be proportional to the direct interaction of prey and predator, i.e., NP. Thus, the available prey for predation is given by , where , and must satisfy [32].

At the hunting time, predators cooperate, which in turn boosts up the predation rate. Thus, the attack rate is the linear increasing function of density of predators in the ecological system [32]. Let c denotes the degree of hunting cooperation among the predator species.

By considering Michaelis-Menten type harvesting, several researchers have showed that harvesting of species has a strong impact on the population dynamics [57, 58]. Here, we introduce nonlinear harvesting of predators by the term where q is the catching capability coefficient, E is the effort applied to harvest the predators, and and are positive constants.

The predator population is assumed to consume the prey items at a constant rate following the Holling type-II response function However, the predation term gets modified due to prey refuge and hunting cooperation of predators, and becomes .

The growth of the predator population depends on the consumption of prey items. Let be the conversion efficiency of prey biomass into the biomass of predator population. Natural death takes place for the predator population at the constant rate .

In [59], authors have considered the predator’s fear and its carry-over effects in the predator–prey model. But, here we focus only on the effect of fear induced by predator and do not consider its carry-over effect in the modeling process.

In [59], authors have considered that the fear of predators suppresses the birth rate of prey species only. But, here we assume that besides reducing the birth rate of prey, the fear of predators enhances the intraspecies competition among the individuals of prey species.

In [59], authors have considered that a constant proportion of the prey species is taking refuge whereas in the present paper, it is assumed that the portion of prey population that take refuge is in accordance with the predator’s density in the ecological system, a more realistic assumption.

Here, we assume that the predator population cooperate with each other during hunting time while the cooperative behavior of predators was not considered in [59].

In the present investigation, it is considered that the predator species are being harvested at a nonlinear rate to avoid the exploitation of resources in the ecological system. But, in [59], harvesting of species was not considered.

In [59], authors have investigated the effects of environmental stochasticity on their proposed predator–prey system. Here, we do not consider the impact of environmental fluctuations rather we focus on the seasonal variations of some the ecologically important model parameters.

All the parameters utilized in the model system (1) are assumed to be positive constants, and the model system (1) is to be analyzed with initial conditions having non-negative values. In Table 1, we have mentioned the ecological meanings of all the parameters appearing in the system (1).

Table 1

Biological meanings of parameters present in system (1) and their numerical values used for simulations

Parameters	Descriptions	Values	Units
r_0	Birth rate of prey	1.1	1/time
r_1	Natural death rate of prey	0.09	1/time
r_2	Death rate of prey due to intraspecific competition	0.06	unit area/number/time
k_1	Level of fear responsible for the reduction in the birth rate of prey	5	unit area/number
k_2	Level of fear responsible for enhancement in the intraspecific competition among prey population	1.8	unit area/number
α	Rate of consumption of prey by predators	0.9	1/time
m	Coefficient of refuge	0.8	unit area/number
c	Degree of hunting cooperation	0.5	unit area/number
a	Saturation constant	1.2	number/unit area
α_1	Conversion efficiency of prey biomass into that of predator	0.8	–
d_1	Natural death rate of predator	0.01	1/time
q	Catching capability	0.25	1/time
E	Harvesting effort	1.5	number/unit area
c_1	Saturation constant	1.1	–
c_2	Proportionality constant	1.5	–

On integrating both equations of the system (1) with the non-negative initial conditions, we getSince , we have for all . Therefore, any solution of system (1) initiating in the positive quadrant of always remain in the same quadrant showing that the interior of the first quadrant of is an invariant set.

Feasibility and stability of system’s equilibria

All the ecologically meaningful equilibria of the model system (1) are listed as follows. We state the following theorem in respect of local stability of the three equilibria of the system (1).

The population-free (trivial) equilibrium point , which always exists. In this equilibrium condition, neither prey nor predator population exists in the ecosystem.

The predator-free (axial) equilibrium point which exists if . Thus, for the feasibility of this equilibrium point, the birth rate of prey population should always be greater than its natural death rate.

The coexistence (interior) equilibrium point , where and are the positive solution(s) of the following isoclines: This equilibrium is ecologically important as both prey and predator populations survive here.

Theorem 1

The population-free equilibrium always exists and is stable (unstable) whenever the predator-free equilibrium does not exist (exists).

The predator-free equilibrium exists if and is stable if the following conditions hold:

The coexistence equilibrium , if exists, is locally asymptotically stable if and only if the following conditions hold: where and are defined in the proof.

Proof

Jacobian matrix of system (1) is obtained as , where

Evaluating the Jacobian matrix J at the trivial equilibrium point gives the two eigenvalues as and Note that the second eigenvalue is always negative while the first one is negative if and positive for . Thus, the equilibrium is stable if , which is opposite to the condition for the feasibility of the equilibrium . Therefore, the equilibrium is stable (unstable) if the equilibrium does not exist (exists). This indicates that these two equilibria are linked via a transcritical bifurcation. Biologically, it interprets that it is possible to observe population-free equilibrium in a predator–prey system of an ecological community until the birth rate of prey population exceeds its natural mortality rate. However, if the birth rate of prey population exceeds its natural mortality rate, then the population-free equilibrium becomes unstable and the equilibrium with prey only appears in the ecosystem.

The two eigenvalues obtained after evaluating the Jacobian matrix J at the predator-free equilibrium are and In view of feasibility of the equilibrium , one eigenvalue is always negative whereas the other eigenvalue will be negative provided the inequality (2) holds. Thus, we can visualize the equilibrium with prey only in some ecosystems.

Jacobian matrix J after evaluating at the equilibrium becomes , where are the entries evaluated at the components of the equilibrium . Characteristic equation of matrix is obtained as follows: where and . In view of conditions stated in (3) and employing Routh-Hurwitz criterion, both the roots of Eq. (4) are either negative or have negative real parts. Thus, the equilibrium , if exists, is locally asymptotically stable if and only if .

Existence of Hopf bifurcation

Here, we investigate the possibility of Hopf-bifurcation through the equilibrium point by taking the birth rate of prey population () as the bifurcation parameter and keeping other parameters fixed. In this regard, we present the following theorem.

Theorem 2

The necessary and sufficient conditions for the occurrence of Hopf-bifurcation through the equilibrium are that there exists a critical value of , say such that

, .

Proof

To have an idea about the nature of equilibrium , we determine the signs of the real parts of the roots of the characteristic Eq. (4). Let be a root of the characteristic Eq. (4). Putting this value in Eq. (4), and separating real and imaginary parts, we getNote that if the characteristic Eq. (4) has a pair of purely imaginary roots, the stability of the system changes through the equilibrium . Setting so , and putting in Eqs. (5) and (6), we getFrom Eqs. (7) and (8), we have and , which implies .

The eigenvalues of the characteristic Eq. (4) areNote that both and are functions of the parameter , when the remaining parameters are assigned fixed values. Let there exists some such that and . Therefore, the positive real parts of these eigenvalues change their sign whenever passes through . As a consequence, system (1) switches its stability provided that the transversality condition is satisfied.

Differentiating Eqs. (5) and (6) with respect to , we get the following after putting From the above two equations, we get thatwhich is non-zero if

Now, we discuss the direction and stability of the bifurcating periodic solutions originating from the equilibrium point via Hopf-bifurcation. For this, we calculate the first Lyapunov coefficient using the approach of [66].

At first, we shift the equilibrium point at the origin by introducing the new variables and . Then, system (1) gets the following form:Taylor’s series expansion of the above system at the point up to terms of order 3, on neglecting higher-order terms, giveswherewithAfter neglecting the higher-order terms of degree four and above, system (9) can be written aswhereThe eigenvector of the matrix corresponding to the eigenvalue at is , where .

Now, we defineLet and , where . Under this transformation, system (10) becomesThis can be written aswhere and are nonlinear in and , and are given bywithNow, based on the normal form (11), we determine the first Lyapunov coefficient aswhere all the partial derivatives are evaluated at the bifurcation point, i.e., . According to [66], the Hopf-bifurcation is supercritical if and subcritical if . For , generalized Hopf-bifurcation occurs at which the Jacobian matrix corresponding to the equilibrium has purely imaginary eigenvalues, and the bifurcation point separates branches of subcritical and supercritical Hopf-bifurcation in the parametric plane.

The nonautonomous model

In system (1), we have assumed that all the parameters describing the model equations are constant. However, in realistic scenarios, parameters involved in ecological models are not constants rather they depend upon several ecological and environmental factors that induce seasonality in the parameters. Here, we extend our autonomous system (1) by allowing the parameters representing the growth rate of prey, fear of predators, refuge behavior of prey, cooperation for hunting and harvesting of predators to vary with time. Thus, after considering the parameters , , , m, c and E as functions of time, the autonomous system (1) is changed into the following nonautonomous form:In system (12), the seasonally forced parameters , , , m(t), c(t) and E(t) are considered to be positive, continuous and bounded functions having positive lower bounds. Also, these parameters are assumed to be periodic. For simplicity, we ignore the phase shifts and consider that these time-dependent parameters have a period of 365 days.

For a continuous -periodic function f(t) on , we introduce the following notations:

Permanence

Theorem 3

System (12) is permanent if the following conditions hold:

From the first equation of the system (12), we getBy the comparison theorem, we have for some From the second equation of the system (12), we haveUsing differential inequality, we getThus, for some , we haveAgain, from the first equation of system (12), we obtainBy the comparison theorem, we have for some Again, from the second equation of system (12), we getAgain, from comparison theorem, for some , we haveObviously, and . Now, we defineThus, we haveTherefore, the system (12) is permanent in view of inequalities in (13).

Global attractivity of positive periodic solution

Lemma 1

For a real number , let f be a nonnegative function which is integrable and uniformly continuous on the interval , then [61].

Theorem 4

A positive periodic solution of system (12) is globally attractive provided,

Let system (12) possesses at least one positive periodic solution , and also we haveSuppose (N(t), P(t)) be any positive periodic solution of system (12). Consider the following functional:Calculating upper-right Dini’s derivative, we getWe have,Thus, we haveWe can write,whereClearly, the function V(t) is monotonic decreasing on the interval . Integrating the above inequality over the interval [0, t], we getHence, from Lemma 1, we haveThus, the positive periodic solution of system (12) is globally attractive.

Next, we show that the globally attractive periodic solution is unique. To this, we assume that is another globally attractive periodic solution of the system (12) with the same period. If this solution is different from the previous one, then there exists at least one such that . This means that . Thus, we havewhich contradicts global attractivity of the periodic solution . Therefore, . Similar arguments can be used for the component also. Hence, the system (12) has a unique positive periodic solution that is globally attractive.

Numerical observations

Here, we report the simulations performed to investigate the behaviors of the autonomous system (1) and the associated nonautonomous system (12) using the MATLAB variable step Runge–Kutta solver ode45. The numerical observations of systems (1) and (12) will reinforce the analytical findings and provide some more insights into the dynamical properties of these systems. Unless otherwise stated, the set of parameter values used for the numerical simulations will be the same as provided in Table 1.

Simulation results of system (1)

Sensitivity of parameters

Following [62], we employ two statistical approaches viz. Latin Hypercube Sampling (LHS) and the Partial Rank Correlation Coefficients (PRCCs) to overcome the uncertainties involved in choosing the values of parameters in system (1). By considering a uniform distribution of our parameters of interest, , , , c, m and E, we run 50 simulations of the model (1) per LHS. For this, we use the baseline values of parameters as given in Table 1, and let them to deviate from their nominal values. The PRCC values using the population of prey as response function are displayed in Fig. 1. It can be seen from the figure that the parameters that have negative correlations with prey species in the ecosystem are and c whereas the parameters with positive correlations are identified as , , m and E. The most influential parameters are coming out be , , m and E that have significant correlations with the prey population of the ecological system.

Fig. 1

Effect of uncertainty of the model (1) on prey population (N). The baseline values of parameters are same as given in Table 1. In the figure, parameters with significant PRCCs with the prey population are marked as (p-value)

In Fig. 2, we have plotted the abundances of prey and predator populations in the model system (1) by varying two parameters at a time viz. , (m, c) and . The figures show surfaces representing the values of the populations at a (dynamic) equilibrium, i.e., the steady state or persistent oscillations. If two surfaces appear in the picture, they indicate the minimum and maximum values that these state variables attain in the limit cycles. If the two surfaces collide, it depicts that a stable steady state is attained, while when they differ, the solution shows oscillatory dynamics. On varying and simultaneously, we find that both prey and predator populations increase in the ecological system as the value of increased whereas increasing values of cause decline in the prey population (see Fig. 2a). The predator population is found to decrease with an increase in prey refuge while its density climb up with increment in the degree of hunting cooperation (see Fig. 2b). The predator’s density also decreases if they are harvested at an increased rate that can be observed in Fig. 2c.

Fig. 2

Variations in the equilibrium abundances of prey and predator populations in system (1) as functions of (a) and , b m and c, and (c) E and . Rest of the parameters are the same values as in Table 1, and the initial conditions are chosen as (1, 1)

Bifurcation results

Now, we see how the dynamics of system (1) changes on variations in the values of some of the ecologically important parameters. We pick here , , , m, c and E as parameters of our interest. Figures 3 and 4 demonstrate phase plots of system (1) for different parametric set ups. At first, we fix and choose different values of . We observe that at , system (1) exhibits stable focus (see Fig. 3a) while on increasing the value of from 0.2 to 1, stability in the system is lost and it becomes unstable by producing limit cycle oscillations (see Fig. 3b). But, the system (1) becomes stable again as the value of increased to 1.5 (see Fig. 3c). Now, we set , and pick different values of the parameters , , m and c for the phase plots. Again, we see changes in stability behavior of the system (1) as the parameter is varied. At , the system (1) exhibits stable focus (see Fig. 3d) whereas the system shows unstable nature at (see Fig. 3e). The system (1) again showcases stable nature for (see Fig. 3f). The effect of parameter can be seen in Fig. 4a and b. We observe from the figures that the system (1) is in unstable mode at , but stability is achieved by the system at . On increasing the value of refuge coefficient (m) from 0.6 to 0.9, the system (1) changes its behavior from unstable mode to stable state (see Fig. 4c and d). We observe unstable behavior of system (1) at while the system achieved stability on decreasing the value of c to 0.1 (see Fig. 4e and f). Finally, we see the behavior of system (1) for different values of the parameter E (see Fig. 4g and h). Here, we have chosen and kept rest of the parameters at the same values as in Table 1. We observe stable nature of system (1) at while instability appears in the system at .

Fig. 3

Phase portraits of system (1) for different parametric set ups: a , , b , , c , , d , e and f . The remaining parameters have the same values as in Table 1

Fig. 4

Phase portraits of system (1) for different parametric set ups: a , b , c , d , e , f , g , and h , . The values of other parameters are the same as mentioned in Table 1

As the birth rate of prey species has a major role in shaping the structure of an ecological community, we further explore the impact of parameter on the dynamics of system (1). In Fig. 5, we draw bifurcation diagram of system (1) by varying the parameter in the interval [0, 1.5]. We plot only the equilibrium value of the predator population and explore stability nature of the system’s equilibria. In the figure, regions divided by vertical dashed lines represent different stability behavior of the system. The figure shows that for lower ranges of , the system settles at the stable population-free equilibrium point . As the value of increases, we observe that the population-free equilibrium becomes unstable and a stable predator-free axial equilibrium appears in the ecological system. This shows the existence of a transcritical bifurcation between equilibria and at the critical point . We further increase the value of parameter , and noted that after a certain value of , the equilibrium loses its stability and the stable coexistence equilibrium comes into the picture. Here, again a transcritical bifurcation occurs between equilibria and at the threshold value . We see that after this critical value of , the trivial and axial equilibria remain unstable for all values of . We keep on increasing the value of and find that there is again a critical value of at which the dynamics of system (1) changes drastically. We find that the system experiences supercritical Hopf bifurcation at ; the coexistence equilibrium becomes unstable and limit cycle oscillations appear around it. But, these limit cycle oscillations do not appear if the parameter exceed another critical value at which the system undergoes another Hopf bifurcation that is subcritical in nature. At this point, the persistent oscillations are killed out and the system manifests a stable coexistence of prey and predator. Thus, we find that lower or higher growth rate of prey population might be beneficial for stability in the ecological system whereas very low (intermediate) growth rate of prey species causes extinct of species (unstable coexistence of prey and predator).

Fig. 5

Bifurcation diagram of system (1) with respect to the growth rate of prey population (). Rest of the parameters are at the same values as in Table 1

In Figs. 6 and 7, we have drawn the bifurcation diagrams of system (1) with respect to the parameters , , m, c and E by varying them one-by-one in a fixed interval. The effect of fear parameter can be seen in Fig. 6a. This figure depicts stable coexistence of prey and predator for lower ranges of . It is apparent from the figure that the system destabilizes on increasing the value of , but for larger values of , the system gets back to a stable state. The critical values of at which the behavior of system changes from stable to unstable to stable are obtained as and , respectively. The fear parameter potentially suppresses the persistent oscillations and drives the system to a stable state as apparent from Fig. 6b. The change in stability is observed at . Figure 7a exhibits that on gradually increasing the values of parameter m, the dynamics of system (1) changes from unstable state to stable state at . In contrast, the cooperative behavior of predator population during the hunting time has tendency to destabilize the system by inducing persistent oscillations as apparent from Fig. 7b. We obtain the critical value of c at which the dynamics of system (1) changes from stable state to unstable mode as . Our simulation results show interesting effect of the harvesting of predator species on the dynamics of system (1), Fig. 7c. The figure demonstrates the stable coexistence of prey and predator populations if the latter are harvested at a lower rate, while the system exhibits limit cycle oscillations for intermediate ranges of the harvesting effort. However, if the predators are harvested at a higher rate, then the existing oscillations are killed out and the system returns to a stable state. In this case also, we obtain the threshold values of the parameter E at which the dynamics of system (1) changes from stable to unstable and unstable to stable modes as and , respectively.

Fig. 6

Bifurcation diagrams of system (1) with respect to (a) and (b) . Rest of the parameters are at the same values as in Table 1. Here, green and black dots represent the upper and lower limits of the oscillation cycles, respectively

Fig. 7

Bifurcation diagrams of system (1) with respect to a m, b c and c E. Rest of the parameters are at the same values as in Table 1 except in c . Here, green and black dots represent the upper and lower limits of the oscillation cycles, respectively

Instead of varying one parameter at a time, simultaneous variations of two parameters will provide more information about the behavior of a dynamical system. Thus, we vary two parameters of system (1) at a time, viz. and (m, c), and draw two-parameter bifurcation diagrams in Fig. 8. In the figures, cyan and magenta regions correspond to the domains for stable and unstable coexistence equilibrium of the system (1), respectively. It is evident from Fig. 8a that for lower values of , the system is at stable state for all values of . We also observe that fixing the value of and increasing the value of , the system’s dynamics changes from stable to unstable to stable modes. Further, it is apparent from the figure that simultaneous increments in the values of parameters and lead to increase in the region of instability of the system around the coexistence of prey and predator. From Fig. 8b, it is clear that there is a range of the parameter m for which the system is in unstable mode irrespective of the values of c. Also, we find a range of the parameter c for which the system (1) changes its dynamics from unstable mode to a stable state on increment in the coefficient of refuge m.

Fig. 8

Two-parameter bifurcation diagrams of system (1) in (a) and (b) (m, c) planes. Here, cyan and magenta regions correspond to the domains for stable and unstable coexistence equilibrium of the system (1), respectively. Rest of the parameters are at the same values as in Table 1

Simulation results of system (12)

Now, we study a wide spectrum of the dynamical behaviors of the nonautonomous system (12), in which some of the model parameters are seasonally forced. By considering periodicity of the rate parameters in an ecological model, we actually incorporate the seasonality of the environment. For our case, we consider that the following model parameters have sinusoidal forms while the others have constant values:The reason behind considering these forms of time-dependent parameters is that by doing so, the functions in time encompass both high and low seasons, that correspond to the periods in which is positive and negative, respectively. Here, the parameters , , , , and , control the strength of the seasonal forcing in the respective parameter. For simplicity, we neglect the phase shift and assume that these seasonal parameters have period of 365 days. Unless otherwise mentioned, the seasonal forcing terms, , , , , and , are assumed to have zero values throughout the numerical simulations.

We set the autonomous system (1) at stable and unstable states, and plot the solution trajectory of the corresponding nonautonomous system (12) for the predator population only (see Fig. 9). We see that if the autonomous system (1) is in unstable mode, then the seasonal variation of the parameter induces higher periodic solutions for , as depicted in Fig. 9a. Now, we keep system (1) at stable state and let the parameter to vary seasonally. We observe that at and , the nonautonomous system (12) exhibits 1−periodic and 2−periodic solutions, respectively (see Fig. 9b and c). We find that if the autonomous system (1) is in oscillatory state, then seasonality in the fear effect responsible for reduction in the birth rate of prey population causes the appearance of bursting patterns in the nonautonomous system (see Fig. 9d). On the other hand, if the autonomous system is at stable state, consideration of the seasonal pattern in the parameter leads to the existence of positive periodic solution (see Fig. 9e).

Fig. 9

Time series solutions of system (12) for different parametric setups: a , , b , , c , , d , , e , , f , , g , , h , , i , , j , , k , , l , , m , , and n , . Other parameters are at the same values as in Table 1

Next, we introduce seasonality in the fear parameter , keeping the autonomous system (1) at stable/unstable mode. We observe the appearances of higher periodic solutions and simple periodic solution in the nonautonomous system if the corresponding autonomous system is at unstable and stable states, respectively (see Fig. 9f and g). Next, we see the effect of seasonal forcing in the pattern of refuge taken by the prey population. We observe that in the unstable scenario, time varying refuge by prey species induces bursting patterns in the nonautonomous system (see Fig. 9h); in the stable case, seasonal forcing of refuge causes appearance of 2−periodic solutions (see Fig. 9i); the higher strength of seasonality leads to bursting patterns (see Fig. 9j). The nonautonomous system (12) exhibits higher periodic solutions and a simple periodic solution due to periodicity in the hunting cooperation if the autonomous system is at unstable and stable states, respectively. In Fig. 9m, we see that if the autonomous system (1) is already in oscillatory state, then harvesting of predators in a periodic manner induces chaos in the ecological system. The occurrence of chaotic oscillations may be explained through incommensurate limit cycles [65]. A simple periodic solution is observed in the nonautonomous system (12) due to seasonal harvesting of predator species if the system without seasonality is at stable state (see Fig. 9n).

The maximum Lyapunov exponent has been identified as an instrumental for the estimation of stability of a nonlinear dynamical system. The main idea is to calculate the average logarithmic rate of separation of two nearby orbits. If the two orbits get too far apart, one of them has to be shifted back to the vicinity of the another one along the line of separation. It has been documented that a chaotic attractor has a positive maximum Lyapunov exponent; a zero maximum Lyapunov exponent corresponds to a bifurcation point; a negative maximum Lyapunov exponent represents a fixed point or a periodic attractor [63, 64]. Here, we compute the Lyapunov exponents by simulating the seasonally forced model (12), and considering the time series solutions of each of the component. Corresponding to Fig. 9m, the maximum Lyapunov exponent of the nonautonomous system (12) is computed and constructed in Fig. 10. The figure confirms chaotic nature of the nonautonomous system (12) as the maximum Lyapunov exponent is positive here.

Fig. 10

Maximum Lyapunov exponent of the nonautonomous system (12) for the parameter values in Fig. 9m. In the figure, positive values of the maximum Lyapunov exponents indicate that the system is in chaotic state

Figure 11 illustrates the global asymptotic stability of the coexistence equilibrium point of the autonomous system (1) and that the positive periodic solution of the nonautonomous system (12) is globally attractive. We can observe in Fig. 11a that the solution trajectories of system (1) starting from four different initial values ultimately converge to the components of the coexistence equilibrium point . This shows that the equilibrium is globally asymptotically stable. Further, the global attractivity of positive periodic solution is apparent from Fig. 11b as the solution trajectories of system (12) originating from four different initial starts ultimately merge with a single positive periodic solution. Thus, Theorem 4 is numerically illustrated in Fig. 11b.

Fig. 11

a Global stability of the coexistence equilibrium of the autonomous system (1), and b global attractivity of the positive periodic solution of the nonautonomous system (12). Parameters are at the same values as in Table 1 except in a and b ,

Conclusion

In the present work, we have investigated the combined effects of fear of predators, prey refuge, hunting cooperation and harvesting of predators on the dynamics of a predator–prey system of an ecological community. We have analyzed the proposed model mathematically, and performed some numerical simulations to explore different dynamical features of the system. Our sensitivity results suggested that a strategy that reduces (enhances) the parameters with negative PRCC values (i.e., and c) will adequately enhance (reduce) the density of prey population in the ecosystem. Simultaneously, a strategy that increases (decreases) the parameters with positive PRCC values (i.e., , , m and E) will be effective in boosting up (lowering down) the density of prey population in the ecological community. Thus, by increasing the refuge tendency of prey and harvesting of predators, the density of prey population can be maintained. Instead, reduction in the fear of predators affecting the growth rate of prey and hunting cooperation of predators can increase the prey’s density in the ecological system.

Our simulation results showed that the birth rate of prey population has potential to regulate the behavior of system around the extinction of species, survival of prey only, and coexistence of prey and predator populations. The fear effects of predators causing reduction in the birth rate and enhancement in the intraspecific competition of prey population have capability to stabilize the system. The property of prey population to take refuge can alter the behavior of system by bringing it back to a stable state from an unstable one. The cooperative behavior of predators has tendency to destroy stability and induce instability in the system. Our results also indicate that harvesting of predator species may play an important role in maintaining stability of the ecological system by terminating persistent oscillations. Since cooperation by predators and refuge by prey are self-adjusted behaviors, manipulating fear by artificial vocalization (or reinforcing predators) and harvesting of predators could be useful for conserving the biodiversity in the ecological system [15].

As seasonal changes of parameters make the system more realistic, we have also studied the behavior of system by letting the parameters standing for the birth rate of prey, fear effects causing decrease (increase) in the birth rate (intraspecific competition) of prey population, prey refuge, hunting cooperation and harvesting of predators to vary with time. The nonautonomous system is analyzed mathematically, and numerical investigations explored the effects of seasonal forcing on the system’s behavior. Note that the presence of positive periodic solution in an ecological model describes an equilibrium situation consistent with the environmental conditions such that the interacting populations survive. We observed that if the autonomous system is in unstable mode, then the seasonal variation of the growth rate of prey induces higher periodic solutions whereas in the stable case, 1−periodic and 2−periodic solutions were seen. We found that if the autonomous system is in oscillatory state, then seasonality in the fear effect responsible for the reduction in the birth rate of prey population causes the appearance of bursting patterns; in the stable scenario, only a simple periodic solution is observed. Seasonal variation in the level of fear responsible for the increase in intraspecific competition of prey species induces higher periodic solutions and simple periodic solution if the autonomous system is at unstable and stable states, respectively. We observed that in the unstable scenario, time varying refuge by prey species induces bursting patterns in the system while in the stable case, seasonal forcing of refuge causes appearance of 2−periodic solutions whereas higher strength of seasonality leads to bursting patterns. Our nonautonomous system exhibited higher periodic solution and a simple periodic solution due to strength of seasonal forcing in the hunting cooperation if the autonomous system is at unstable and stable states, respectively. We have seen that if the autonomous system is already in oscillatory state, then harvesting of predators in a periodic manner induces chaotic dynamics in the system, which is explained through the appearance of incommensurate limit cycles. A simple periodic solution is observed in the nonautonomous system due to seasonal pattern in the harvesting of predators if the autonomous system is at stable state.

The chaotic nature of nonautonomous system is confirmed by the positivity of the maximum Lyapunov exponent. We also illustrate the global stability of the coexistence equilibrium (in autonomous case) and positive periodic solution (in nonautonomous case). Overall, we observed that our autonomous as well as nonautonomous systems exhibited complex dynamics due to changes in some model parameters, which can be managed by regulating the strengths of seasonal forcing in the respective parameters. Overall, our findings demonstrated that the combined effects of fear of predators, prey refuge, hunting cooperation and harvesting of the predator species together with seasonal patterns of the parameters can play an important role in preserving the biodivesrity in the ecological system.