A delay non-autonomous model for the combined effects of fear, prey refuge and additional food for predator

In this paper, we investigate the combined effects of fear, prey refuge and additional food for predator in a predator–prey system with Beddington type functional response. We observe oscillatory behaviour of the system in the absence of fear, refuge and additional food whereas the system shows stable dynamics if anyone of these three factors is introduced. After analysing the behaviour of system with fear, refuge and additional food, we find that the system destabilizes due to fear factor whereas refuge and additional food stabilize the system by killing persistent oscillations. We extend our model by considering the fact that after sensing the chemical/vocal cue, prey takes some time for assessing the predation risk. The delayed system shows chaotic dynamics through multiple stability switches for increasing values of time delay. Moreover, we see the impact of seasonal change in the level of fear on the delayed as well as non-delayed system.


Introduction
Rapid change in climatic conditions, human activities, upwelling resource pulses etc., have caused extinction of many species around the globe. To maintain the ecological balance and biodiversity, there is an urgent need of conservation of these species. In these situations, only statistical data is not enough to maintain the ecological balance. Indeed, mathematical models can help to better understand the behaviour of ecosystem. The most important factor in a predator-prey model is to observe the type of interactions between prey and predator, affecting dynamical behaviour of the system. There are two types of interaction − one prey dependent which is a function of prey only and the other one is predator dependent which depends on both species. The most widely used functional response in describing predator-prey interactions is the Holling type II functional response, in which per capita predation rate is an increasing, smooth, and saturating function of prey density. However, the prey-dependent functional responses fail to model the interference among predators and have been facing challenges from the biological as well as physiological communities. Some biologists have argued that in many situations, especially when predators have to search for food, the functional response in a predator-prey model should be predator dependent. There are much significant evidences which suggest that predator dependence in the functional response occurs quite frequently in laboratory and natural systems. In order to mediate theoretical and experimental perspectives, Beddington [8] and DeAngelis et al. [21] considered a functional form of prey consumption rate that depends on both the species. They proposed the form, g = px/(ax + by + c), which is similar to the Holling type II functional response with the extra term by in the denominator, which interprets the mutual interference among predators. This functional form is known as Beddington-DeAngelis functional response.
It is well documented that Beddington-DeAngelis type functional response is more acceptable than other available response functions [14,64], and it has the ability to take care of a number of ecological mechanisms. Several mathematical models including the Beddington-DeAngelis type functional response have been investigated [22,26,43], and the results showed that this type of functional form produces very rich and biologically reasonable dynamics. Cantrell and Cosner [14] have studied a predator-prey model with Beddington-DeAngelis functional response, and analysed various dynamical properties like stability, limit cycle, etc. Dimitrov and Kojouharov [22] found that mutual interference between predators can alone stabilize predator-prey interactions. Growing biological and physiological evidences suggest that predator-dependent functional response is more appropriate when predators have a highly competitive food searching process in nature due to limited resources [6,24]. In view of these facts, in the present article, we consider Beddington-DeAngelis functional response as the interaction term between prey and predator.
Over the last few decades researchers paid attention on direct predation only in a predator-prey model as it is easy to observe in real world [1][2][3]68]. But, recent studies have shown that indirect impact of predator greatly affect the population dynamics [36,37,52,54,61,65,66]. Fear effect plays a pivotal role among the indirect effects. In any habitat, prey feel predation pressure as a result they exhibit anti-predator behaviour like foraging, habitat change, different physiological change, etc. Scared species forage less. Due to fear of predator, prey acquires different mechanisms like reduction in mating, starvation, migration from vulnerable place to invulnerable which actually lessen their growth rate. For example, scared Mule deer opt for lesser foraging to avoid the acute predation risk of mountain lions [5]. Creel et al. [17] observed that elks alter their reproductive physiology due to the fear of wolves. Such complex strategies play a central role in the evolutionary processes. Candolin [13] observed that threespine stickleback males regulate their reproductive decisions by assessing the threat of predation. The empirical evidences by Zanette et al. [74] established that offspring generation of song sparrow is demographically impressed and reduced by 40% due to the perceived threat of predation. In their experiment, they actively prevented the chances of direct killing with the help of electric fencing and seine netting, and imitated the predation risk by broadcasting predator call playbacks over an entire breeding season of song sparrows. Modelling the impact of fear effect on the dynamic interaction of prey and predator has gained commendable attentions among researchers after the seminal work of Wang et al. [69]. Panday et al. [53] investigated a three species food chain model by considering that the growth rates of prey and middle predator are reduced due to fear of middle and top predators, respectively. In [40], authors discussed the control of chaotic dynamics in a three species food chain model with fear. Zhang et al. [75] showed that the fear effect not only reduces the predator population but also stabilizes the system by excluding periodic solutions.
If predator can not reach to a portion of prey, it means prey species are able to take refuge. Refuge is an anti-predator behaviour which can help to prolong a predator-prey interaction by reducing the chance of extinction due to over predation. When prey experience attack cue from predator, they take refuge by changing their functioning like their body colour, shape, size, habitat, etc. Most of the adult preys are capable to take refuge. The first demonstration of this effect was given by Gause in his experiments with the protozoan, Didinium (predator) and Paramecium (prey) [28]. Mite predator-prey interactions often exhibit spatial refugia which afford the prey some degree of protection from predation and reduce the chance of extinction due to predation. In 1970, Connell [16] found that adults of the barnacle, Balanus glandula, were safe from predatory snails of the genus, Thais, in intertidal areas which are submerged only for a short period of time during high tide. The mature barnacles were restricted to refuge in the areas where Thais was present but not where Thais was absent, suggesting that the refuge prevents Thais from eliminating Balanus. Dolbeer and Clark [23] found that snowshoe hare (prey) prefers forested habitat with dense undercover, where they are relatively safe from their predator (e.g. Canada lynx). Additional food is any kind of food or species other than focal/basal prey. Additional food to the predator acquire less attention than focal prey. Due to fear and refuge, prey availability reduces to a great extent to the predator and as a consequence predator population goes to extinction. To maintain ecological balance or to sustain predator population in the system, additional foods are provided to predator. Joint effects of reduction in prey refuge and the presence of appropriate amount of additional food can control the prey population [7,15,20,30,59,62]. Das and Samanta [19] studied the effects of fear and additional food for predator in a fluctuating environment. The findings of Wang et al. [67] showed that for a certain level of fear, prey refuge is beneficial for the coexistence of prey and predator populations. Recently, Mondal and Samanta [51] studied the dynamics of a delayed predator-prey system with nonlinear prey refuge under the influence of fear effect and additional food for predator. Note that in these studies, authors have considered that a constant proportion of prey population is taking refuge, but, in a realistic scenario, refuge depends on predator density [49,63].
It is well understood that time delay occurs often in almost every biological situation. Several studies have been done with retarded systems which explored different dynamics of the systems including periodic oscillation, persistence, bifurcation and chaos of populations [9,10,18,29,46]. The process of perceiving predator signals through chemical and/or vocal cues and the changes in life history and behavioural responses in the prey population conceals some time delay. Panday et al. [55] investigated the impact of such time delay in a two dimensional predator-prey system. They observed that for the gradual increase in the magnitude of delay, the system dynamics switches multiple times between stable focus and limit cycle oscillations, and ultimately enters into the chaotic regime. Mondal et al. [50] investigated the impact of gestation delay in a predator-prey system with fear effect and additional food for predator. Kumar and Dubey [39] studied the fear effect and prey refuge in a predator-prey system with gestation delay; they observed that the system enters into chaotic regime for larger values of gestation delay. Seasonality can be defined as the regular and periodic changes of variable on an annual time scale [70]. Seasonal variables relevant in ecological systems not only include temperature, photoperiod but also include rainfall, wind, human activity, etc. [45,60]. The recognition of these varied drivers of ecological systems shows the ubiquity of seasonality. When the environmental fluctuation is taken into account, a model must be non-autonomous. The level of fear may vary due to change in predation pressure, intraspecific competition or resource distribution among the individuals [25,32]. A predator-prey model with seasonal variation of fear effect is studied by Roy and Alam [58]. Sk et al. [62] investigated the joint effects of fear, refuge and additional food on the dynamic interactions of prey and predator by allowing the respective parameter to vary seasonally; complex dynamics such as higher periodic solutions and bursting patterns were observed due to variations in the seasonally forced parameters.
In the present investigation, our first goal is to explore the combined effects of fear of predator, prey refuge and additional food for predator with modified Beddington type functional response. The fraction of prey population taking refuge is assumed to depend on the predator density. Our another goal is to investigate the effect of time delay involved in the process of perceiving predator signals and behavioural responses by the prey population. Finally, we see the impact of seasonal variation of fear of predator in the systems with and without time delay.

The mathematical model
To construct the model, we consider water bugs as prey and fish as predator [44]. In the presence of fish, the bugs become scarce, and measurable densities (a proportion of water bugs depending on predator density) are restricted to areas with thick aquatic vegetation. As a result, scarce prey which are spending time in refuges also typically costs prey in terms of reduced feeding, mating opportunities, competition, etc. So, for fish, available water bugs reduced in a great amount. Due to less prey (water bugs), predators (fish) interfere with one another. Therefore, fish need additional food to persist in the system. We can also consider snowshoe hare (prey) and Canada lynx (predator) as another example of species for our proposed model. Let at any time t > 0, N (number per unit area) be the density of prey population and P (number per unit area) be the density of predator population in the region under consideration. In the modelling process, the following assumptions are made.
(1) Prey population grows with r 0 as its intrinsic growth rate and r 1 as mortality rate due to intraspecific competition when there is no predation and fear effect. The reason behind considering the intraspecific competition among the species of prey population is that for high density of prey population and less availability of resources, the individuals of prey population compete with each other for the available resources. (2) The effect of fear not only influence the reproduction of prey but also has the potential to affect the death rate of the prey population [42,72,73]. Thus, the growth rate of prey population should be multiplied by a factor that decreases with the increase in level of fear and the density of predator population. In view of this, the growth rate (r 0 ) is multiplied by a factor 1 1+kP , which accounts for the cost of anti-predator behaviour of prey due to fear of predation [62,69]. Here, the constant k measures the level of fear.
(3) The prey population has capability to take refuge, which provides them a degree of protection against predation. However, the refuge behaviour of prey population is affected by the number of predator population present in the area. So, the proportion of prey population taking refuge is considered to be a function of predator population [49,63]. (4) The interactions between prey and predator populations are assumed to be governed by modified Beddington type functional response [15]. (5) The predator population is supposed to be supplied with an additional food of biomass A, which is assumed to be distributed uniformly in their habitat. The additional food is assumed to be non-reproducing and is supplied at a constant rate. The number of encounters per predator with additional food is proportional to the density of additional food. Moreover, providing additional food to predator population reduces the rate at which predators consume the prey items. (6) The density of predator population decreases in the region due to natural death at a constant rate d.
Keeping these assumptions in mind, our proposed model is represented by the following system of differential equations: Throughout this paper, we consider the acceptable range 0 ≤ (1 − mP) ≤ 1, i.e. P ≤ 1 m . Parameters involved in system (1) are assumed to be constant and have positive values. The initial conditions are given by In Table 1, we have mentioned the biological meanings of the parameters involved in system (1).

Positivity and boundedness of system's solutions
Positivity: The positivity of state variables (system's solutions) of system (1) implies that the population never go to extinction and always survive in the system. To prove the positivity of the solutions of system (1), we integrate both equations of system (1) over the limit 0 to t, and get Since N(0), P(0) > 0, we have N(t), P(t) > 0 for all t ≥ 0. Therefore, any solution initiating in the positive quadrant of R 2 always remains therein. Hence, the interior of the first quadrant of R 2 is an invariant set. Boundedness: In theoretical ecology, boundedness of a system implies that the system is well behaved. Boundedness of the solutions entails that none of the interacting populations grow exponentially for a long-time interval. Proof: We note from the first equation of system (1) that Therefore, lim sup t→∞ N(t) ≤ r 0 r 1 .
Thus, for an arbitrary small 1 > 0, there exists T 1 > 0 such that for some M 1 Again, from the second equation of the system (1), for all t > T 1 , we have Applying the theory of differential inequality [41], we obtain Hence, as t → ∞, for m > 0, we have 0 ≤ P(t) ≤ β(αM 1 +ηA) Therefore, solutions of the system (1) are bounded above. Thus, the region of attraction for all solutions of system (1) initiating in the positive quadrant is given by the following set: which is closed and bounded in the positive cone of the two dimensional plane.

System's equilibria
Due to nonlinearity of model system (1), it is not possible to find exact solutions to the system. Instead, we determine the long-term behaviour of the system. In general, a nonlinear system either gravitates towards an equilibrium point or it blows up. An equilibrium point represents the rest state of a dynamical system. These points can be obtained by putting the growth rate of different variables of model system (1) equal to zero. The equilibria of system (1) are obtained as follows.
(2) The predator-free equilibrium (4) The coexistence equilibrium E * = (N * , P * ), where and P * is the positive solution(s) of the following biquadratic equation: where Clearly, A 4 is positive. Thus, Equation (4) has at least one positive root if A 0 is negative.

Remark 3.1:
In the equilibrium E 0 , where there is no species in the system, it should not be possible to observe it in the real ecosystem. That is, this equilibrium point should be unstable. In the equilibrium E 1 , where there is no predator population in the system, it should be possible to observe it in some special circumstances. That is, this equilibrium can be stable under certain conditions. The equilibrium E 2 represents the case when ecosystem is free from prey population and only predators survive. This is possible because the predators depend upon additional food apart from their focal prey. For the existence of this equilibrium, the amount of additional food must be in accordance with the ability to detect and handle additional food, conversion efficiency, and death rate of predator. Finally, the equilibrium E * , where both prey and predator are present, is very common in natural ecosystem. Thus, it should be stable in some ecosystems. It is worthy to note that the stability of coexistence equilibrium will be helpful for the biological conservation programs.

Stability analysis
In this section, we perform the local stability analysis of the equilibria of system (1). The local stability analysis characterizes whether or not the system settles to the equilibrium point if its state is initiated close to, but not precisely at a given equilibrium point. The equilibrium point is said to be locally asymptotically stable if there is a neighbourhood of the equilibrium point such that for all initial starts in that neighbourhood, the system approaches to the equilibrium point as t → ∞ [57]. The local stability of all equilibria of system (1) is stated in the following theorem.
(2) The equilibrium E 1 is stable if the following condition holds: (3) The equilibrium E 2 , if exists, is stable if the following conditions hold: Proof: Stability of the equilibria E 0 , E 1 and E 2 can be shown by evaluating the Jacobian of system (1) at these equilibria and determining signs of the real parts of the eigenvalues. At the equilibrium E * , the Jacobian of system (1) reduces to The corresponding characteristic equation is obtained as where B 1 = −(a 11 + a 22 ) and B 2 = a 11 a 22 − a 12 a 21 . By Routh-Hurwitz criterion, we claim that roots of Equation (7) are either negative or have negative real parts if and only if B 1 , B 2 > 0.

Remark 3.2:
The first point of above theorem indicates that it is not possible to observe the equilibrium without prey and predator. The second point of this theorem shows that the predator-free ecosystem is achieved only if the growth rate of predator due to consumption of prey and additional food is lower than its natural death. It is apparent from the third point of the above theorem that for stability of prey-free equilibrium, the fear of predator and the growth rate of predator due to additional food should be high, predator consume more prey, growth rate of prey is less due to fear and also prey take less refuge. The last point of Theorem 3.2 tells that if the initial state of system (1) is near the equilibrium point E * , then the solution trajectories not only stay near the equilibrium E * for all t > 0, but also approach to the components of equilibrium E * as t → ∞. Thus, if the initial values of state variables N and P are close to N * and P * , respectively, then the system (1) will eventually get stabilized provided B 1 > 0 and B 2 > 0. That is, small perturbations in system's variables do not affect stability of the system at the coexistence equilibrium point.

Hopf-bifurcation analysis
Now, we consider m as a bifurcation parameter and investigate the occurrence of Hopfbifurcation through the coexistence equilibrium point E * , keeping the remaining parameters fixed. Hopf bifurcation with respect to m means that there exists a critical value of m at which systems's stability switches and a periodic solution arises from stable state. Existence of Hopf bifurcation ensures survivability of all species in oscillatory mode, which actually fit with our ecological subsystem. In this context, we have the following theorem.

Theorem 3.3: If there exists a critical value of m (say m
] m=m * = 0. Then, the system (1) undergoes Hopf-bifurcation through the equilibrium E * at m = m * .

Proof:
To determine the nature of equilibrium E * , we require the sign of the real parts of the roots of the characteristic Equation (7). Let λ(m) = u(m) + iv(m) be one of the eigenvalues of the characteristic Equation (7). Substituting the value of λ in Equation (7), and separating real and imaginary parts, we get 2uv A necessary condition for the change of stability through the equilibrium E * is that the characteristic Equation (7) should have purely imaginary roots. We set m = m * such that u(m * ) = 0, and put u = 0 in Equations (8) and (9). Then, we have From Equations (10) and (11), The eigenvalues of the characteristic Equation (7) are Here, B 1 and B 2 are the functions of the parameter m, when other parameter values are fixed. Moreover, we assume there exists some m = m * such that B 1 (m * ) = 0 and B 2 (m * ) > 0. Therefore, the positive real parts of these eigenvalues change the sign when m passes through m * . Consequently, the system switches its stability provided that the transversality condition is satisfied.
Differentiating Equations (8) and (9) with respect to m, and put u = 0, we have

Effect of time lag involved in the fear of predation
In this section, we modify our previous model system (1) by incorporating a discrete time delay. After sensing the chemical or vocal cue, prey takes some times for assessing the predation risk. Therefore, the prey does not respond instantaneously to the population size of the predator species, rather there must be some time lag [11,55]. In order to incorporate this time lag, we consider that at time t, the fear of predator is in accordance with the predator density at time t − τ (for some τ > 0). Incorporating this time lag in system (1), we get the following system of delay differential equations: The initial conditions for the system (12) take the form +0 and denotes the norm of an element κ ∈ C + by κ = sup −τ ≤ζ ≤0 {|κ 1 (ζ )|, |κ 2 (ζ )|}. For biological feasibility, we further assume that κ i (0) ≥ 0 for i = 1, 2.

Stability analysis in the presence of time delay
Now, we study the stability behaviour of the system (12) in the presence of time delay around the coexistence equilibrium E * . Linearizing the system (12) about the coexistence equilibrium point E * = (N * , P * ), we get where R = r 11 r 12 Here, n and p are small perturbations around the equilibrium point (N * , P * ).
The characteristic equation for the linearized system (14) corresponding to equilibrium E * is given by where For τ = 0, all the roots of Equation (15) are either negative or have negative real parts provided b 1 > 0 and b 0 + c 0 > 0 while for τ > 0, it has infinitely many roots. By Rouche's Theorem and continuity in τ , stability of Equation (15) will change if the real part of the leading roots crosses imaginary axis, i.e. if Equation (15) has purely imaginary roots. Hence, putting λ = iω (ω > 0) in Equation (15), and separating real and imaginary parts, we get Squaring and adding Equations (16) and (17), we get Substituting ω 2 = ψ in Equation (18), we obtain the following equation: where a 1 = b 2 1 − 2b 0 and a 0 = b 2 0 − c 2 0 . Note that Equation (19) has exactly one positive root if a 0 is negative, and two positive roots if a 0 is positive and a 1 is negative. Theorem 4.1: Suppose that the equilibrium E * exists and is locally asymptotically stable for τ = 0. Also, let ψ 0 = ω 2 0 be a positive root of Equation (19). Then, there exists τ * = min n≥0 τ * n such that the equilibrium E * of the system (12) is asymptotically stable for 0 ≤ τ ≤ τ * and unstable for τ > τ * , where Proof: Since ψ = ω 2 0 is a solution of Equation (19), the characteristic Equation (15) has a pair of purely imaginary roots ±iω 0 . It follows from Equations (16) and (17) that τ * n is a function of ω 0 for n = 0, 1, 2, 3, . . .. Thus, if the system (12) is locally asymptotically stable around the coexistence equilibrium E * for τ = 0, then by Butler's lemma, the equilibrium E * will remain stable for τ < τ * and unstable for τ > τ * such that τ * = min n≥0 τ * n provided the following transversality condition holds: Differentiating Equation (15) with respect to τ , we get Clearly, (ω 2 0 ) = 0, since ω 2 0 is a simple positive root of Equation (19). Therefore, the transversality condition is verified and hence Hopf-bifurcation occurs at τ = τ * , i.e. a family of periodic solutions bifurcate from the equilibrium E * as τ passes through τ * [31]. (19) has two positive real roots ψ + (corresponds to ω 2 + ) and ψ − (corresponds to ω 2 − ), where ψ + > ψ − , i.e. the characteristic Equation (15) has two pair of purely imaginary roots ±iω ± , then we have

Remark 4.1: If Equation
for n = 0, 1, 2, 3, . . .. If the following transversality conditions are satisfied, then there exists a positive integer q such that there are q switches from stability to instability and back to stability. In other words, the equilibrium E * is unstable. Therefore, Hopf-bifurcations occur around the equilibrium E * at τ * n ± for n = 0, 1, 2, 3, . . ..
The system (21) is of the formẊ and a bilinear inner product where γ (σ ) = γ (σ , 0). Clearly, D(0) and D * are adjoint operators. We know that ±iω 0 τ * are eigenvalues of D(0). So, they are also eigenvalues of D * . Now, we search for the eigenvector of D(0) and D * corresponding to iω 0 τ * and −iω 0 τ * , respectively. Now, we assume that q(σ ) = (1, u) T e iω 0 τ * σ and q * (s) are the eigenvectors of D(0) and D * for the eigenvalues iω 0 τ * and −iω 0 τ * , respectively. Then, we have D(0)q(σ ) = iω 0 τ * q(σ ). Using the definition of D(0) and from (24), we have We need to determine the value of H in such a way that q * (s), q(σ ) = 1, q * (s), q(σ ) = 0. Therefore, which implies H = 1 1+u * u+τ * u * s 12 ue −iω 0 τ * . To describe the centre manifold C 0 at ν = 0, we compute the coordinates by using the same notations and procedures as proposed by Hassard et al. [34]. Let X t be the solution of Equation (21) at ν = 0. Define On the centre manifold C, we have U and U are local coordinates for centre manifold C 0 in the direction of q * and q * . Here, V is real when X t is real. For solution X t ∈ C 0 of Equation (21), since ν = 0, we havė Then, from Equation (28), we obtain Thus, from Equation (29), we have Comparing the coefficients of above equations with Equation (29) To calculate the value of g 21 , we need to compute the values of V 20 (σ ) and V 11 (σ ). From Equations (25) and (28), we havė where Expanding the above series and comparing the corresponding coefficients, we get From Equation (32), for −1 ≤ σ < 0, we have Comparing the coefficients with Equation (33), we obtain Using Equation (34) and first equation of (36), we havė where E 1 = (E (1) 1 , E (2) 1 ) ∈ R 2 is a constant vector. Similarly, from Equation (34) and second equation of (36), we obtain where E 2 = (E (1) 2 , E (2) 2 ) ∈ R 2 is a constant vector. The constant vectors E 1 and E 2 should be chosen appropriately.
Thus, from Equations (37) and (38), we can easily find out V 20 and V 11 . Hence, we can find each g ij defined in Equation (31) in terms of delay and other biological parameters. Finally, we can compute the following quantities: (12) Regarding global stability of the equilibrium E * in the presence of time delay, we have the following theorem.

Theorem 5.2:
The sufficient condition for global asymptotic stability of the equilibrium E * of system (12) is min{l 1 m 1 , l 2 m 2 } > 0, where Proof: To show the global stability of the equilibrium E * , we consider the following transformations: , P(t) = P * e y(t) . Thus, the system (12) becomes Therefore, the equilibrium point E * converts into trivial equilibrium (x, y) = (0, 0) for the transformed system (50). Let V 1 (t) = |x(t)|. Now, calculating the upper right Dini derivative of V 1 (t) along the solution of system (12), we get We consider the functional as P * |e y(s) − 1| ds dp.
Calculating the Dini derivative of V 11 (t), we get Thus, Let V 2 (t) = |y(t)|. Calculating right hand Dini derivative of V 2 (t) along the solution of system (12), we have Therefore, where l 1 and l 2 are given by (46) and (47), respectively. As system (12) is permanent, we have for all t > T Using the mean value theorem, we have where N * e ξ 1 (t) lies between N * and N(t), and P * e ξ 2 (t) lies between P * and P(t). Therefore, where l = min{l 1 m 1 , l 2 m 2 }. Hence, the trivial equilibrium point of the transformed system (50) is globally asymptotically stable if l > 0. That is, the coexistence equilibrium point (N * , P * ) of system (12) is globally asymptotically stable if l > 0.

Effect of seasonality in delayed system
In system (12), we have assumed that the parameter representing the effect of fear is constant. But, in realistic scenario, this parameter is not constant. Indeed, the fear of predator depends upon several ecological and environmental factors, and hence vary with time [58,62,63]. Therefore, we extend our delayed system (12) by letting seasonal effect in the fear parameter k. Thus, we have the following delay non-autonomous system: Here, we assume that k(t) is positive continuous and bounded function with positive lower bound. We neglect phase shift and assume for simplicity the seasonal variation in the level of fear by letting the parameter k as a periodic function of time with period ω.

Existence of periodic solution
In this section, we show that system (51) possesses at least one positive periodic solution.
To this, we use the following lemma from [27].
then system (51) has at least one positive ω-periodic solution, where L 1 , H 1 , L 2 and H 2 are defined in the proof.
Proof: Let us assume first that k(t) be a continuous ω-periodic function on R, and we denote We consider the following transformations: Then, system (51) becomes It is obvious that, if (y 1 (t), y 2 (t)) T is an ω-periodic solution of the system (54), then (N(t), P(t)) T = (e y 1 (t) , e y 2 (t) ) T is a positive ω-periodic solution of the system (51). Now, we define the set and the norm (y 1 (t), y 2 (t)) T = max Note that Y and Z are both Banach spaces with respect to the above norm . . Then, and dimKer(L) = 2 =codimIm(L).
Since Im(L) is closed in Z, L is a Fredholm mapping of index zero. It is easy to show that P and Q are continuous projections such that, Im(P) = Ker(L), Im(L) = Ker(Q) = Im(I − Q).

Now, from Equations
Therefore, for any t ≥ 0, from Equation (59), we get for some H 2 Again, from Equations (56) and (60), we have Thus, for any t ≥ 0, from Equation (58), we obtain for some L 1 Finally, from Equations (57) and (60), we obtain Thus, for any t ≥ 0, from Equation (59), we have for some L 2 Now, for some C 1 and C 2 , we have max t∈[0,ω] Clearly, C 1 and C 2 are independent of λ. Define C = C 1 + C 2 + C 3 , where C 3 is sufficiently large such that for each solution (y * 1 , y * 2 ) T of the following algebraic equations: satisfies (y * 1 , y * 2 ) T < C provided that system (61) has one or a number of solutions. Let us define = {(y 1 , y 2 ) T ∈ X : (y 1 , y 2 ) T < C}. Clearly, satisfies the first condition of Lemma 6.1.
Let (y 1 , y 2 ) T ∈ ∂ ∩Ker(L) = ∂ ∩ R 2 , (y 1 , y 2 ) T is a constant vector in R 2 with (y 1 , y 2 ) T = |y 1 | + |y 2 | = C. Then, Hence, the second condition of Lemma 6.1 is also satisfied. Now, define the homomorphism J : Im(Q) → Ker(L). Since Im(Q) =Ker(L), J is the identity mapping, and hence J is an isomorphism. Therefore, by direct calculation, we have Therefore, the third condition of Lemma 6.1 is verified. Thus, by Lemma 6.1, we can conclude that system (54)

Numerical simulations
For the numerical simulations of systems (1), (12) and (51), we choose the parameter values as given in Table 1 and show qualitative results. Unless it is mentioned, the values of parameters used for the numerical simulations are the same as in Table 1. For solving the ordinary differential equations (1) and the system (51) in the absence of time delay, we have used MATLAB solver ode45 whereas for solving the delay differential equations (12) and the non-autonomous delay differential Equations (51), we have used MATLAB solver dde23.

Impacts of fear, refuge and additional food on the equilibrium densities of prey and predator
At first, we see the impacts of fear, refuge and additional food on the equilibrium values of prey and predator populations in system (1). For this, we vary k, m and A in fixed intervals, and plot the equilibrium densities of prey and predator populations in system (1), Figure 1. It is apparent from Figure 1(a) that both prey and predator populations decrease with increase in fear factor. If fear effect is very high, the abundances of prey and predator become very low in the ecosystem. The refuge property of prey although increases prey density, but high refuge in prey population leads to decline in the predator density (see  Table 1. Figure 1(b)). We observe that moderate values of refuge can balance the prey as well as predator populations. Moreover, providing additional food to the predator population help to boost its density whereas diminish the prey with little bit fluctuations in the latter (see Figure 1(c)). Next, we vary two parameters of system (1) at a time viz. (m, k), (A, k) and (A, m), and see the evolutions of prey and predator populations, Figure 2. The pictures show surfaces representing the equilibrium level or the maximum and minimum values of prey and predator populations when they oscillate. These plots indicate that on varying m and k simultaneously, the prey population suddenly increases just after a fixed value of m. The predator population decreases on enhancement in prey refuge, but after a certain value of refuge coefficient, the predator population suddenly increase. We observe that although fear factor causes rapid decline in prey and predator populations, providing additional food accelerates their growths. Coupling A with m, we note that providing additional food fails to balance the growths in prey and predator whereas high refuge prevents the extinction of prey population.  Table 1.

Sensitivity analysis
In order to characterize the connection between the parameters of system (1) and its solutions, the differential sensitivity analysis of system (1) is performed [12,47,48]. The semi-relative and logarithmic sensitivity solutions with respect to three most sensitive relevant parameters k, m and A corresponding to the state variables (prey and predator) are computed. The former provide the amount the state will change when a parameter is doubled whereas the latter indicate percentage change in the state variable on doubling a model parameter. The sensitivities of model variables are plotted in Figure 3. In Figure 3 Table 1 except m = 0.9 and A = 0.7.
for predator. This may be due to the fact that additional food to the predator contains low nutrients and also acquire less attention than focal prey. The fear of predator and the tendency of taking refuge highly influence the prey density. From Figure 3(b), which depicts logarithmic sensitivity solutions, it is evident that doubling the parameter k leads to 39.37% and 1.81% decrements in the prey and predator populations, respectively, in 80 days. On doubling the refuge parameter m, we observe 46.43% increase and 97.43% decrease in the prey and predator populations, respectively, in 80 days. The densities of prey and predator populations are boosted by 0.33% and 0.14%, respectively, in 80 days on doubling the value of A. As in the case of semi-relative sensitivity solutions, here also both the populations are least sensitive to the additional food, however, it fosters their growth. Both of these sensitivity solutions indicate that when fear of predator affects prey and hence the predator, prey refuge supports the survival of prey whereas providing additional food to the predator enhances the survivorships of prey and predator populations. (1) At first, we keep system (1) in the absence of fear, refuge and additional food, i.e. we set k = m = A = 0. We see the impact of theses three ecologically important parameters by introducing them in the system one-by-one. We find that for k = m = A = 0, the system (1) is in oscillatory state (see Figure 4(a)) whereas the oscillations are terminated and the system returned to stable coexistence by incorporating only the fear factor (see Figure  4(b)) or prey refuge (see Figure 4(c)) or additional food (see Figure 4(d)). These results explored the stabilizing roles of fear [69], prey refuge [59] and additional food [59] in a real predator-prey system. Now, we see the dynamics of system (1), where the fear, prey refuge and additional food are simultaneously present. The phase portraits of system (1) are depicted in Figure 5. We observe that the system exhibits limit cycle behaviour for higher values of fear factor whereas stability is achieved on decreasing the value of k. That is to say, fear effect destabilizes the system. In contrast, prey refuge and additional food stabilize the system, i.e. on increasing the value of any of them, the system returns to stable state from persistent oscillations which resemble the findings of Samanta et al. [59]. In Figure 6, we draw bifurcation diagrams of system (1) with respect to the three interesting parameters, k, m and A. Rising of periodic solutions from stationary solution in N−P plane is apparent from Figure 6(a). This figure clearly shows the occurrence of Hopf-bifurcation on increasing the value of k. On the other hand, the existing periodic solution disappear and the system reaches to a  stable state on increasing the value of m (see Figure 6(b)). Similar effect is observed for the additional food; increasing values of A evacuated the persistent oscillations and settle the system to a stable coexistence (see Figure 6(c)). Thus, fear of predator destabilizes the system through Hopf-bifurcation whereas prey refuge and additional food push back the system to stable state by terminating the limit cycle oscillations.

Impact of delay in fear of predator
Here, we report the simulations performed to investigate the behaviour of delay differential equations (12). At first, we fix the values of parameters at which the system (12) is at stable state in the absence of time delay and gradually increase the magnitude of delay parameter τ . We plot the solution trajectories of system (12) in Figure 7. We see that at τ = 4, the system shows limit cycle oscillations (see Figure 7(a)) whereas at τ = 15, the system exhibits stable focus (see Figure 7(b)). On increasing the value of τ to τ = 25, we observe that the system again showcases limit cycle oscillations (see Figure 7(c)). Finally, we pick τ = 65 and find that the dynamics of system (12) is chaotic (see Figure 7(d)); in the figure incommensurate limit cycles indicate the chaotic behaviour of the system [33,35]. Next, we set system (12) at oscillatory state in the absence of time delay and see the behaviour of system at τ = 4, 15, 25 and 65 (see Figure 8). It is apparent from the figure that the system exhibits limit cycle oscillations at τ = 4 and τ = 15 while at τ = 65, the nature of system is chaotic. Interestingly, at τ = 25, the system demonstrates periodic oscillations of period two. Thus, for a certain range of delay parameter, the peak of oscillations jumps between two values each year. To get a clear idea of stability change for increasing values of τ , we draw bifurcation diagram of system (12) for τ ∈ (0, 70], Figure 9. It is clear from the figure that for initial values of τ , the system shows limit cycle oscillations and becomes stable as τ increased. We observe several stability switches on gradual increase in the values of τ , and the system ultimately enters into chaotic regime through period doubling oscillations. Figure 10 demonstrates the stability regions of system (12) in (k,τ ) and (A,τ ) parametric planes. In the figures, blue regions represent the domains of stability for the coexistence equilibrium E * and green regions represent the domains of instability. From Figure 10(a), it is clear that for lower values of k, system (12) is always stable at the coexistence equilibrium for all values of time delay. Also, for higher values of k, the coexistence equilibrium is unstable irrespective of the values of τ . However, for moderate values of k, there exist stability switches for increasing values of τ . It is inferred from Figure 10(b) that for higher values of A, stable coexistence is achieved, i.e. delay does not affect stability behaviour of the system. On the other hand, for lower values of A, system switches its stability properties finite number of times with gradual increase in the values of time delay.  Table 1, and the initial conditions are chosen as (2, 1).
To confirm chaotic nature of system (12), we plot the maximum Lyapunov exponent of the system at τ = 65 by using the algorithm of [56,71], Figure 11. The maximum Lyapunov exponent has been proven to be most useful dynamical diagnostic for chaotic systems. It represents the average exponential rate of convergence/divergence of nearby orbits in the phase space. The idea of calculating its value is to follow two nearby orbits and calculate their average logarithmic rate of separation. Whenever they get too far apart, one of the orbits has to be moved back to the vicinity of the other along the line of separation. If the attractor is chaotic, the maximum Lyapunov exponent is positive; for a bifurcation point the maximum Lyapunov exponent is zero; the negative values of maximum Lyapunov   Table 1. exponent correspond to a fixed point/periodic attractor. To plot the maximum Lyapunov exponent, the delayed system (12) is first simulated, then by considering the time series solutions of each component of the system, the maximum Lyapunov exponent is computed. In Figure 11, the positive values of maximum Lyapunov exponent demonstrate the chaotic regime of system (12).   Table 1, and the initial conditions are chosen as (2, 1). In the figure, positive values of the maximum Lyapunov exponent confirm the occurrence of chaotic oscillation.

Impact of seasonality
Now, we incorporate seasonal variation in the fear parameter by considering k(t) as a periodic function of time. The functional form of k(t) is taken as k(t) = k + k 11 sin(ωt) with period of one year. In this way, the periodic function in time encompasses both high and low seasons, that correspond to the periods in which sin(ωt) is positive and negative, respectively. The parameter k 11 (0 < k 11 < k) controls the strength of the seasonal forcing. We investigate the impact of seasonality in the system with and without time delay. We plot solution trajectories for the predator population only (the prey population is found to exhibit similar behaviour), Figure 12. First, we fix τ = 0 in system (51) and see how seasonal forcing is regulating the dynamics of system (see Figure 12( a-c)). For a particular parameter set up, the autonomous system (1) shows stable focus whereas the non-autonomous system (51) without delay shows periodic solution (see Figure 12(a)). Figure 12(b) shows that solutions initiating from three distinct points ultimately go to a unique periodic solution, i.e. the positive periodic solution is globally attractive. It is apparent from Figure 12(c ,d) that when systems (1) and (12) show limit cycle oscillations, the non-autonomous system (51) shows higher periodicity. Figure 12(e) depicts that chaotic nature of system (12) turns into higher periodic behaviour when seasonal effect is introduced in the fear factor. For a particular parameter set up, the chaotic behaviour of system (12) remain unchanged although the fear parameter vary seasonally (see Figure  12(f)).

Conclusion
In every ecological system, prey populations are always threatened by predators and perceived predation risk. The fear of predation alters the behaviour of prey population. In predator-prey systems, the prey population shows a variety of different defense mechanisms to escape from predation; one of them is prey refuge, which decreases the predation rate. In contrast, the additional food always provides benefits to the predator population and helps the predators to persist in the system. In this paper, we have studied the combined effects of fear, refuge and additional food on the interacting dynamics of a predator-prey system. Our sensitivity results showed that both prey and predator populations decrease with increase in the fear factor, and the abundances are too low for higher level of fear of predator. However, the joint effects of reduction in prey refuge and the presence of appropriate amount of additional food can control the prey population. Numerical results showed that fear induces oscillation in the system, which can be terminated by increasing the refuge coefficient or amount of additional food to the predator. However, if the system exhibits persistent oscillations in the absence of fear, prey refuge and additional food for the predator, then stability can be regained by including the role of any of these three ecological factors. The existence of limit cycles in predator-prey systems can be used to explain many real-world oscillatory phenomena, such as the Canadian lynx snowshoe hare 10-year cycles.
Further, we have incorporated time delay involved in the process of perceiving predator's signals through chemical and/or vocal cues and the changes in life history and behavioural responses in the prey population. We have analysed the delayed system for the existence of Hopf bifurcation by taking time delay as a bifurcation parameter. We also obtained analytical conditions determining nature of the bifurcating periodic solutions. Moreover, sufficient conditions are derived under which the coexistence equilibrium is globally asymptotically stable in the presence of time delay. We gradually increase the delay parameter and observed that our delayed system experiences multiple stability switches through chaotic dynamics; for larger values of time delay, the system remains chaotic. Ecologically, presence of chaos affects ecosystem features such as predictability, species persistence [4], and biodiversity [38]. Moreover, we observed that presence of delay in fear of predator has no effect on the system's dynamics if the level of fear is too low or too high. However, for moderate level of fear or additional food, the system experienced finitely many stability switches on increasing time lag in the fear of predator.
Furthermore, we consider seasonal variation in the level of fear, and assume that it is periodic function of time with a period of one year. We saw the impact of seasonal changes in fear effect on the dynamics of systems with and without time delay. Our nonautonomous system exhibits positive periodic solution for the set of parameter values at which the corresponding autonomous system is stable. We showed that the positive periodic solution is globally attractive. However, if the autonomous system with or without time delay shows limit cycle oscillations, the non-autonomous counterpart generates higher periodic solutions. The outcomes of the present investigation are very important for understanding complex behaviours of dynamical systems, and have wide spread applications in conservation biology.

Disclosure statement
No potential conflict of interest was reported by the author(s).