Non-smooth dynamics emerging from predator-driven discontinuous prey dispersal

In ecology, the refuge protection of the prey plays a significant role in the dynamics of the interactions between prey and predator. In this paper, we investigate the dynamics of a non-smooth prey–predator mathematical model characterized by density-dependent intermittent refuge protection of the prey. The model assumes the population density of the predator as an index for the prey to decide on when to avail or discontinue refuge protection, representing the level of apprehension of the prey by the predators. We apply Filippov’s regularization approach to study the model and obtain the sliding segment of the system. We obtain the criterion for the existence of the regular or virtual equilibria, boundary equilibrium, tangent points, and pseudo-equilibria of the Filippov system. The conditions for the visibility (or invisibility) of the tangent points are derived. We investigate the regular or virtual equilibrium bifurcation, boundary-node bifurcation and pseudo-saddle-node bifurcation. Further, we examine the effects of dispersal delay on the Filippov system associated with prey vigilance in identifying the predator population density. We observe that the hysteresis in the Filippov system produces stable limit cycles around the predator population density threshold in some bounded region in the phase plane. Moreover, we find that the level of apprehension and vigilance of the prey play a significant role in their refuge-dispersion dynamics.


Introduction
In an ecosystem, the prey population exhibit a variety of mechanisms to avoid predation [29]. In general, prey prefers to stay away from their predators by utilizing a secondary habitat as refuge protection which is otherwise inaccessible to their predators [36,37]. The prey refuge protection can reduce the predation pressure and subsequently protects the prey population from depletion [15]. Apart from this, a sufficiently large prey refuge can even drive the predator population to extinction [5]. Since refuges are safe for the prey but restrict their foraging time or mating opportunities than in open habitat, there is a reduction in the growth rate in the prey population during refuge [1,18,30,35]. To balance the trade-off between food and safety, prey mostly prefers to stay in the open habitat when the predator population is relatively low and switches to the refuge at a sufficiently high predator population [30]. When the predator population density becomes low, the prey emerges from refuge to take advantage of feeding and mating opportunities in the open habitat [18,30]. In this paper, we investigate the anti-predator behaviour of prey in response to predation risk versus available resources.
The environmental attributes that portray the behaviour of the prey such as habitat usage and foraging strategy are the consequences of natural selection [17]. The behavioural responses of the prey to its predator examine how a prey incorporates the risk of predation in its foraging decisions. The predators create stimuli that alert the target prey of their presence, creating certain apprehension in the prey population [10,11,33]. While different behavioural traits of prey are favoured in different environments, certain behaviour type affects the survival fitness of the prey under the threat of predation [8]. Some prey appears to ignore predator cues when there is a substantial presence of the predators, referred to as bold prey [4,16]. An innately bold prey often stays in the open habitat for a longer period by prioritizing resource acquisition at the cost of enhanced predationdriven mortality [16,38], whereas an extremely apprehensive individual being overtly cautious, often overstays in a refuge by missing out the forage and mating opportunities [38]. Regardless of the level of apprehension, some prey may require more information on the predators to make decisions for switching the habitat [14]. In the context of safeguarding against the possible predation, the prey needs some time to ascertain the threat employing vigilance [32]. The vigilance of prey is a behaviour related to monitoring the potential threats caused by their predators. Vigilance can be pre-emptive or reactive. Prey with pre-emptive vigilance is overtly vigilant, whereas for reactive vigilance, prey chooses a course of action only after the predator density crosses some threshold value. Overt vigilant prey is expected to make faster habitat-choice decisions compared to their less vigilant counterpart.
Prey vigilance is mostly attributed to the apprehension of the prey [20,21,41]. Researchers [2,20] observed that uncertainty in the risk of predation can lead to an increase in the level of vigilance among prey. In an unpredictable environment, apprehensive prey not being able to access the risk of predators, often overestimates the degree of predation risk, thereby increasing their vigilance level. However, the researcher [2] stated that prey vigilance may not always be an indicator of the state of apprehension. As observed by the researcher [2], in a predictable environment, apprehensive prey when faced by predators with reliably known cues, becomes less vigilant, whereas an unknown cue perceived in a predictable environment will increase the level of vigilance of the prey without any change in prey apprehension. It has also been observed that hunger can keep the vigilance of prey at its low ebb even if the prey is highly apprehensive to its predator [3]. Apart from this, highly apprehensive prey when moving in a group can become less vigilant with an increase in their group size [3]. In this case, we focus on studying the cases when bold prey becomes hypervigilant or highly apprehensive prey becomes less vigilant.
Many researchers have studied prey-predator interactions by assuming that a constant proportion of prey is protected from predation by the refuge [12,13,27,28,39,42]. In reality, the prey dispersal can be densitydependent [24]. To the best of our knowledge, the preydensity-dependent refuge has been studied in relatively few papers [25,31], whereas the study of predatormediated density-dependent refuge is relatively new. In this paper, we assumed that the prey remains in some open habitat whenever the predator population density is less than some threshold value, representing the loss of prey apprehension when the predators are rare; above the threshold value, a density-dependent refuge of the prey takes place. It seems reasonable to assume that the switching of the habitat of the prey is not instantaneous, rather, there can be some lag due to the prey's behaviour and associated dispersal time. It is important to study the impact of prey behaviour on the dynamics of the system. Here we further modify the model to include the hysteresis effect and an associated switching delay and address the implications of the habitat switching delay in the dynamics of the system. The incorporation of the hysteresis effect in the system allows identifying the level of vigilance of the prey in identifying the threshold population density of the predator.
The main emphasis of our study will be in exploring the dynamic behaviour of the non-smooth system with prey and predator density-dependent refuge and identify the refuge time spent by the prey based on predator-mediated different behavioural responses of the prey.

Model equations
In this study, we have considered a two-species predatorprey model in an open habitat under the assumption that a portion of the prey population takes refuge from the predators when the density of the predators goes above a certain threshold level L; below this threshold value, the prey does not go to the refuge and prefers to stay in the open habitat which is more favourable compared to the refuge. The smaller values of L represent prey to be more apprehensive of predation risk, whereas the larger values of L imply that the prey is less scared to predation. Let x(t) and y(t) represent the densities of prey and predator population, respectively, at time t in the open habitat. For y > L, we assume the amount of the prey x r (t) that takes refuge depends on the prey population density in the open habitat at time t given by x r (t) = αx(t) β+x(t) , where α (0 ≤ α ≤ 1) is the maximum capacity of the refuge and β (β ≥ 1) is the half-saturation coefficient of the refuge. Therefore, the per-capita fraction of prey in the refuge is The prey in the open habitat is assumed to have a logistic growth in the absence of the predators with intrinsic growth rate r and carrying capacity K . In the presence of the predators, the rate of consumption of the prey by the predators is assumed to follow linear functional response with the per-capita consumption rate c. The growth of the predator can be split up into three portions; conversion rate (with conversion efficiency e), density-dependent death rate γ due to intraspecific competition (overcrowding rate), and the natural mortality rate δ.
For 0 < y < L, the prey prefers to stay in the open habitat. In this case, the equations governing the dynamics are: where x(0) ≥ 0 and y(0) ≥ 0. For y > L, a density-dependent refuge protection of the prey is applied. In this case, the density of the prey available for the consumption of the predator in the open habitat at time t is (x(t) − x r (t)), and is governed by the following equations: where x(0) ≥ 0 and y(0) ≥ 0.
x ≥ 0, y ≥ 0 . Therefore, if the predator density falls below the threshold L, that is, H (Z ) < 0, then the prey stays in the open habitat. If H (Z ) > 0, then prey starts taking refuge. The mapping H : R 2 + → R is smooth and the normal vector H Z (Z ) perpendicular to the switching line [22] with discontinuous righthand sides given by (1) and (2) can be written as: where The non-overlapping regions S 0 and S 1 are separated by the switching line such that S 0 ∪ ∪ S 1 = R 2 + . In the region S 0 , the prey does not move to the refuge, while in the region S 1 , a portion of the prey population goes to the refuge. The switching line divides system (3) into two subsystems defined in the regions S 0 and S 1 . The dynamics of the vector fields F (Z ) ( = 0, 1) is analysed using the classical methods to study systems of ordinary differential equations.
3 Dynamics in the region S ( = 0, 1) Clearly, system (3) has two different structures, a refuge-free system (for = 0) and a system with refuge protection for the prey (for = 1). First, we show that the solutions of (3) are nonnegative and contained in a compact subset of the phase space if the initial values are positive.
From (3) we have dx dt = x f 1 S (x, y; t) and dy dt = y f 2 S (x, y; t), where Therefore, This implies, all solutions of (3) remain within (x, y) ∈ R 2 : x ≥ 0, y ≥ 0 starting from an interior point of it. Therefore, R 2 + = (x, y) ∈ R 2 : x > 0, y > 0 is an invariant region, and as long as x(t) > 0 and y(t) > 0 for all t, the local existence and uniqueness properties hold in R 2 + . We now prove that the solutions of (3) with initial values in R 2 + are bounded. From the prey equation we have the differential By comparison, it follows that lim t→∞ sup x(t) ≤ K , and so, corresponding to 1 > 0, there exists t 1 > 0 such that x(t) < K + 1 , for all t ≥ t 1 . Thus, we get . Therefore, lim t→∞ y(t) ≤ ecK −δ γ . The following result states the boundedness of the solutions of system (3).

Lemma 3.2 All the nonnegative solutions of system (3) lie in the region
The results from Lemma 3.1 and 3.2 imply that system (3) is dissipative. Next, we find the equilibria of system (3) and establish the parameter conditions for which stability of the subsystems can be ensured. System (3) defined in the regions S has the non-trivial nullclines f i S = 0 (i = 1, 2 and = 0, 1) (cf. Fig. 1).
An equilibrium E of system (3) is said to be a regular H (E) < 0, then the equilibrium E of system (3) is referred to as a virtual equilibrium. In other words, if any steady-state solution E S of system (3) defined in the region S lies in the region S , it becomes a regular equilibrium point, otherwise a virtual equilibrium point. The intersection of the nullclines of system (3) yields the following equilibria: The interior equilibrium of system (2) Since the boundary equilibria E 0 S 1 and E 1 S 1 of the subsystem (2) are always virtual, these have no biological significance and have nothing to do with the dynamics of system (3). Therefore, by using eigenvalue analysis, we investigate the local behaviour of the subsystems (1) and (2) The eigenvalues of the Jacobian matrix of the subsystem (1) at E 0 S 0 are r and −δ. Therefore, E 0 S 0 is always a saddle point. Since E 0 S 1 is always virtual, it follows that the organisms of system (3) are unlikely to go to extinction.
At E 1 S 1 , the eigenvalues of the Jacobian matrix of the subsystem (2) are −r and δ * − δ, where δ * = ecK . Therefore, the predator-free equilibrium E 1 S 1 is locally asymptotically stable if δ > δ * . This implies, a high whereas E * S0 is virtual, and so no sliding mode exists. The prey and predator-nullclines are represented by blue and red curves, respectively. The solid blue and red curves represent the nullclines within the domain of definition of the system, whereas the dashed curves represent the position of the nullclines outside their domain of definition mortality rate of the predators can lead to its extinction. The stability at E 1 S 1 implies the non-existence of E * S 1 . Lemma 3.3 System (3) defined in the region S ( = 0, 1) has no limit cycles.
Proof We use the Bendixson-Dulac criterion [34] to verify the non-existence of limit cycles in (3). For this, let us consider a Dulac function ψ(x, y) = 1 x y , where x = 0 and y = 0.
Then, we have If possible, let there exists a closed trajectory C of the plane autonomous system in S ( = 0, 1). Let A ⊂ S be the interior of C. Then, by Green's Theorem we have This, however, leads to a contradiction as ∇. (ψ F ) < 0, and so, there can be no such closed trajectory C. Therefore, system (3) has no limit cycles.
We note that E 0 S 0 is always a saddle point and if δ < δ * is satisfied, then the equilibrium E 1 S 0 is a repeller. Also, system (3) is dissipative. Therefore, for the persistence of system (3) at E * S ( = 0, 1), we must have δ < δ * . Since system (3) has no periodic solution and is persistent for δ < δ * , a sufficient condition for the stability at E * S is given by the following Lemma:

Lemma 3.4 System (3) defined in the region S is locally asymptotically stable at E
The boundary equilibrium E 1 S 0 collides with the interior equilibrium E * S 0 when δ = δ * and exchanges stability.
Therefore, we investigate the qualitative changes in the dynamic behaviour of system (3) defined in the region S 0 by varying the parameter δ. At δ = δ * , the Jacobian matrix J 1 S 0 at E 1 S 0 has a zero eigenvalue. Thus, at δ = δ * , the equilibrium point E 1 S 0 becomes non-hyperbolic and the usual linearization technique fails to infer the criteria of stability at E 1 S 0 . In the following Lemma we show that under some certain conditions a transcritical bifurcation occurs when δ crosses δ * . Lemma 3.5 If r = cK , the subsystem (1) undergoes a transcritical bifurcation at E 1 S 0 when δ crosses δ * . Proof The eigenvectors of subsystem (1) corresponding to the zero eigenvalue for J 1 S 0 and J 1 T S 0 are given by Therefore, if r = cK , system (3) defined in S 0 satisfies the Sotomayor theorem [34] for transcritical bifurcation at E 1 S 0 when δ crosses δ * . From Lemma 3.5 it follows that if r = cK is satisfied, the subsystem (1) switches its stability from E * S 0 to E 1 S 0 as δ is increased through δ * . A one-parameter bifurcation plot with δ as a bifurcation parameter is shown in Fig. 2. We observe that, for small values of δ, the equilibrium E * S 0 becomes virtual and E * S 0 is regular and stable. With an increase in δ, the predator population density decreases and once the predator population y * 1 drops below L the regular stable equilibrium E * S 1 becomes virtual. In this case, system (3) stabilizes at the pseudo-equilibrium E P 0 on the sliding domain S ⊂ (as defined in Sect. 4) until δ reaches Fig. 2). In the following section, we discuss the existence of a sliding domain on the switching boundary and analyse the dynamics of the Filippov system (3) on the sliding domain.  (3) with δ as a bifurcation parameter Filippov system (3) on the switching line , we need to employ Filippov's convex method [22] and Utkin's equivalent control method [40] and determine sufficient conditions for the occurrence of a sliding mode on .
By the Lie derivative method, the directional derivative of H along the vector field F (Z ), can be written ., . denotes the inner product in R 2 .
A sliding mode exists if the vectors F (Z ) ( = 0, 1) on of the two subsystems are pointed toward one another. The sliding domain S can be expressed as S = {Z ∈ : Since follows that e = φ. Therefore, the system (3) has no escaping region, the uniqueness of trajectories of (3) is ensured.
The sliding region is given by s = {Z ∈ : ec and Thus, we have S = e s = s and as long as x L 0 < x < x L 1 holds, any trajectory that hits the sliding segment S , stays there. The motion on S repre-sents that the preys are hesitant to take refuge protection.
The crossing region is given by c = c− c+ , where c− = Z ∈ : 0 ≤ x < x L 0 is the crossing region pointing left to the sliding segment and c+ = Z ∈ : x > x L 1 is the crossing region pointing right to the sliding segment.
On the crossing region c , one of the vector fields points toward the switching boundary and the other pointing away from . We use Filippov's convex method to extend the vector field of the system (3). For this, we define Then the sliding vector field F S (Z ) on S becomes a convex combination of the two vector fields F 0 (Z ) and F 1 (Z ) of the subsystems. The convexification of system (3) leads to a system of convex differential inclusion [23,26] The existence of solutions of the differential inclusion (4) can be guaranteed in the sense of Filippov. We note that on the sliding segment S , λ(Z ) depends on x only. We re-write the Filippov system (4) and use the equivalent control method [19,40] to determine the sliding mode dynamics on S as follows: where the control function is given by On the sliding plane S , we have H (Z ) = 0 and therefore we obtaiṅ Solving for λ(x), from (6) we get λ(x)= (ecx−γ L−δ)(β+x) ecαx . The governing equations in the sliding region S , obtained by using the equivalent control method [19,40] In the event that the vector fields F 0 (Z ) and F 1 (Z ) of the system (3) point toward , a trajectory hitting is compelled to evolve within the sliding region S governed by the differential equation (7) along the vector field F S (Z ) until reaching some point on S where one of the two vector fields F 0 (Z ) or F 1 (Z ) changes its direction.
There are several different types of equilibria of the Filippov system (4) that includes regular (or virtual) equilibrium, pseudo-equilibrium, boundary equilibrium, and tangent point. We will focus on studying the dynamics of system (3) on the sliding segment S .
The dynamics defined on the sliding segment S is determined by the differential equation (7), which admits pseudo-equilibrium of the Filippov system (4) (cf. Fig. 4). The pseudo-equilibria E P i = (x P i , L) (i = 0, 1) of system (4) satisfies the differential equations (7), where  A boundary equilibria of system (4) occurs when an equilibrium point of any one of the two subsystems intersects the switching line . Therefore, the boundary equilibria satisfy F (Z ) = 0, H (Z ) = 0, and are given by E B S = x * , L ( = 0, 1). The boundary equilibria E B S exist on S when L approaches the critical value L = y * ( = 0, 1). The existence of boundary equilibrium of the system (4) is shown in Fig. 4b, d.
At the tangent points on S , one of the vectors F 0 (Z ) or F 1 (Z ) is tangent to the discontinuity boundary. The tangent points satisfy the equations L F H (Z ) = 0 ( = 0, 1) and act as a demarcation between the sliding and crossing regions (cf. Fig. 4). Solving these two equations yields the tangent points 1). In either cases the trajectories ofŻ (t) = F (Z ), Z ∈ S , passing through the point E T S stay in the region S ( = 0, 1).
In the following section we examine the discontinuityinduced bifurcation followed by the changes in the sliding mode domain.

Sliding bifurcation
Sliding bifurcations of the Filippov system (4) are related mainly to the value of L and on the equilibria of its subsystems. We investigate the sets of regular, virtual and pseudo-equilibrium bifurcations of system (4) by employing numerical methods.

Bifurcation sets of equilibria
We see that the existence and the type of interior equilibria E * S ( = 0, 1) are dependent on δ and L. We The three curves i (i = 1, 2, 3) divide the δ − L parameter plane into six regions, and the existence or coexistence of stable regular or virtual interior equilibria E * S ( = 0, 1) is pointed out in each of the regions (cf. Fig. 3). The curves 1 and 2 are used to determine the relationship between pseudo-equilibrium E P 0 and the sliding line S given by x L 0 ≤ x P 0 ≤ x L 1 . The curve 3 is used to study the stability of the equilibria E * S 0 and E 1 S 0 of subsystem (1). We see that for 0 ≤ δ < min δ * 0 , δ * , the equilibrium E * S 0 of the subsystem (1) becomes virtual. For 0 ≤ δ < δ * and 1 (δ, L) > 0, the equilibrium E * S 1 of the subsystem (2) becomes virtual. Therefore, for 0 ≤ δ < min δ * 0 , δ * and 1 (δ, L) > 0, the equilibria E * S ( = 0, 1) become virtual followed by the presence of an attract-ing sliding region bounded by the curves 1 and 2 . The boundary equilibria E B S ( = 0, 1) appear on the curves 1 and 2 , respectively.
When an ordinary equilibrium of (4) collides with the sliding segment , it gives either an equilibrium transition bifurcation or a non-smooth fold bifurcation. For the equilibrium transition bifurcation, an ordinary equilibrium collides with a virtual pseudo-equilibrium, giving rise to a virtual ordinary equilibrium and a pseudo-equilibrium. For the non-smooth fold bifurcation, an ordinary equilibrium collides with a pseudoequilibrium, giving rise to a virtual ordinary equilibrium and a virtual pseudo-equilibrium. Boundary equilibrium bifurcation Boundary equilibrium bifurcation (BEB) in the Filippov system (4) involves the possible topological changes when a real or a virtual equilibrium point collides with the discontinuity boundary due to the changes in some parameter of the system. When E * S is a regular equilibrium, a boundary node bifurcation occurs for system (4) once the tangent point E T S and the regular equilibrium E * S collides at L = y * . When E * S is a virtual equilibrium, a boundary equilibrium bifurcation occurs when the pseudo-equilibrium E P 1 and the tangent point E T S  [7] merge with the equilibrium E * S at L = y * . Therefore, we take L as a bifurcation parameter.
The choice of the values of L is dependent on behavioural pattern of the prey in the presence of the predator. The lower values of L will indicate that the prey is more apprehensive to the predators, whereas the higher value of L indicates that the prey is relatively bold. Under the set of parameter values as given in Table 1, we obtain E * S 0 = (2.1266, 0.5357) and E * S 1 = (2.4669, 0.3724). From Fig. 4a it is observed that for L < y * ( = 0, 1), the virtual equilibrium point E * S 0 , the regular equilibrium E * S 1 , the visible tangent point E T S 1 and the invisible tangent point E T S 0 coexist. This indicates that if the prey is very apprehensive to their predators, they will take refuge protection even at a low population density of the predators. In this case, the prey will have to take refuge continuously. Due to this, the predation pressure on the prey will become low, resulting in a high prey population density.
For relatively bold prey, the threshold population density of the predators to start taking refuge by the prey becomes relatively high. As L is increased through L = y * 1 , the stable regular equilibrium E * S 1 collides with the visible tangent point E T S 1 , generating the stable boundary equilibrium E B S 1 (cf. Fig. 4b). For y * 1 < L < y * 0 , E * S becomes virtual ( = 0, 1) followed by the emergence of pseudo-equilibrium E P 1 on the sliding segment (cf. Fig. 4c). In this case, both the tangent points E T S ( = 0, 1) are virtual. In this case, the prey takes refuge protection intermittently and so, the preydensity does not get compromised even with a relatively high population density of the predators. This suggests that the bold prey can prevent their population density from possible depletion even with a high population density of the predators by taking refuge protection intermittently.
Another boundary equilibrium bifurcation of Filippov system (4) occurs when L is increased through L = y * 0 . In this case, the virtual equilibrium E * S 0 collides with the pseudo-equilibrium E P 1 at the invisible tangent point E T S 0 , generating the stable boundary equilibrium E B S 0 (cf. Fig. 4d). Thereafter, for any further increase in L, the regular equilibrium E * S 0 , the virtual equilibrium E * S 1 , and the visible tangent point E T S 0 coexist, whereas the pseudo-equilibrium E P 1 disappears from the sliding segment. This indicates that the bold prey (the prey that does not take refuge protection even at a high population density of the predators) is mostly prone to sustained loss in their population driven by high predation pressure.
In order to show the changes in the sliding mode domain with the changes in L and δ, we plot a codimension two bifurcation diagram with L and δ as bifurcation parameters (cf. Fig. 5a). The bifurcation plot confirms the possible existence of a shorter sliding mode domain when the prey is bold and mortality rate of the predator is low. For a more apprehensive prey, a relative larger sliding domain can exist for high mortality rate of the predators. To investigate the sensitivity of the parameters L, c, δ and K for the existence of a pseudoequilibrium on the sliding domain, we plot bifurcation diagrams as given in Fig. 5b-d. From Fig. 5b and c it follows that the interval of existence of the pseudoequilibrium increases with the decrease in the predation rate, while from Fig. 5d we see that an increase in the prey carrying capacity can lead to longer sliding domain even with low apprehension of the prey. Fig. 4 A boundary node bifurcation with L as a bifurcation parameter for δ < δ * . a For 0 < L < y * , E * S0 is virtual and E * S1 is regular. b For L = y * 1 , a boundary node bifurcation occurs when E * S1 merges with the tangent point E T S1 , generating the boundary point E B S1 . c For y * 1 < L < y * 0 , E * becomes virtual followed by the existence of stable pseudo-equilibrium point E P 1 on the sliding line. d For L = y * 0 , a boundary node bifurcation occurs when E * S0 merges with E P 1 at the boundary point E B S0 on the switching line Pseudo-saddle-node bifurcation The existence of a pair of pseudo-equilibria on the sliding segment and their subsequent annihilation due to the changes in some threshold parameter values occurs via a saddle-node bifurcation. Taking K = 8 and L as a bifurcation parameter, we observe that a generic pseudo-saddle-node E P SN appears on the sliding segment when L crosses the critical value L K sn = 1.2142 (cf. Fig. 6b, d. For L < L K sn , a pseudo-saddle E P 0 and a pseudo-node E P 1 appear on the sliding segment S (cf. Fig. 6a). In this case, the orbits initiating from I 0 S 0 = (4.05, 1.165) and I 0 S 1 = (3.9, 1.25) converge to the pseudo-node E P 1 from the left, whereas the orbits initiating from I 1 S0 = (4.25, 1.117) and I 1 S 1 = (4.2, 1.25) converge to the pseudo-node E P 1 from the right. All the orbits plotted in Fig. 6a, initiated from the left of I 0 S 0 and I 0 S 1 , converge to the stable regular equilibrium E * S 0 . At L = L K sn , the generic pseudo-saddle-node E P SN is halfstable. The orbits starting from I 0 S 0 and I 0 S 1 converge to the regular equilibrium E * S 0 , whereas orbits starting from I 1 S 0 and I 1 S 1 converge to the pseudo-saddle-node E P SN from the right (cf. Fig. 6b). For L > L K sn , pseudo-  Table 1. a A pseudo-saddle (E P 0 ) and a pseudo-node (E P 1 ) exist on the sliding segment S for L < L K sn . b For L = L K sn , the generic pseudo-saddle-node E P SN exists on S . The evolu-tion of the sliding modes (grey regions) with the changes in L for c all parameters values as in Table 1, where a pseudo-saddlenode bifurcation occurs on the crossing segment c− and for d for K = 8 and other parameters values as in Table 1, where a pseudo-saddle-node bifurcation occurs on the sliding segment S equilibrium does not exist (cf. Fig. 6c, d) and the regular equilibrium E * S 0 becomes globally asymptotically stable. The effect of L on the variation of the sliding mode is shown in Fig. 6c, d. Next, we choose γ = 0.015 and fix all other parameters as in Table 1. By varying L, for L < L γ sn = 1.4148, it is observed that a pseudo-saddle E P 0 close to the tangent point E T S 0 and a pseudo-node E P 1 coexist on S Fig. 7 Pseudo-saddle-node bifurcation of system (4) for γ = 0.015 and L as a bifurcation parameter, other parameters values as in Table 1. a The existence of a pseudo-saddle (E P 0 ) and a pseudo-node (E P 1 ) on the sliding segment S for L < L γ sn . b For L = L γ sn , the generic pseudo-saddle-node E P SN exists on S . The evolution of the sliding modes (grey regions) with respect to the changes in L for (c) γ = 0.015 and for d c = 0.5, other parameters values as in Table 1 (cf. Fig. 7a, c). As L is increased through the critical value L = L γ sn , E P 0 and E P 1 collide and generate the pseudo-saddle-node E P SN (cf. Fig. 7b). As L is increased further, the pseudo-saddle-node E P SN disappears (cf. Fig. 7c) and the regular equilibrium E * S 0 becomes globally asymptotically stable. The changes in the length of the sliding segment for different values of L is shown in Fig. 7c. Similarly, for c = 0.5 and other parameter values as in Table 1, the changes in the sliding mode due to the changes in the values of L is given in Fig. 7d.

Prey vigilance and hysteresis effect
In the Filippov system (4), the dispersal of the prey from the unprotected zone to refuge protection is assumed to be instantaneous whenever the predator population crosses some critical threshold value. In reality, the movement of prey species between the two habitats involves some travel time. Apart from this, the reaction time of the prey plays a pivotal role in the switching delay. The slow or less vigilant prey will take more dispersion time than the fast or more vigilant prey. Due to this delay in dispersion, the population threshold of the predators to enter the prey-refuge differs from the population threshold of predators to return to the unprotected zone. This delay in the dispersion of the prey leads to increased predation pressure. To monitor for such possibilities, we introduce the hysteresis effect on system (4) which allows a period of over-exploitation of prey species without compromising its density. The modified Filippov system with hysteresis effect is given bẏ with the control U h H (Z ) defined as where σ > 0 is the hysteresis parameter. The hysteresis parameter σ , representing the level of vigilance of the prey, satisfies 0 < σ < |L − y * | so that the interior equilibria E * S are always virtual ( = 0, 1). The smaller values of σ represent hypervigilant prey, whereas the larger values of σ implies that the prey is less vigilant. For σ > |L − y * |, the collapse of prey stock becomes irreversible due to sustained predation pressure. The system (8) has the two switching lines 0 = Z ∈ R 2 + : H (Z ) = −σ and The switching lines 0 and 1 generate the following regions in the state space, 0 = Z ∈ R 2 + : H (Z ) < −σ , 1 = Z ∈ R 2 + : H (Z ) > σ and 2 = Z ∈ R 2 + : −σ ≤ H (Z ) ≤ σ . The trajectories that start from 0 and evolve toward 1 has the threshold 1 , whereas the trajectories that start from 1 and evolve toward 0 has the threshold 0 (cf. Fig. 8). Let the vector H Z (Z ), normal to 1 and 2 , is oriented in the direction from 0 to 1 . The curves V S (Z ) = L F H (Z ) = 0 ( = 0, 1), subdivide the state phase into the following four regions: We use the geometrical approach proposed by Boukal and Krivan [9] to prove the existence of a trapping region in the state space where the trajectories of system (8) enter with a transversal motion.
On the switching line 0 , L F 0 H (Z ) > 0 implies y > ecx−δ γ and so, the trajectories starting from 0 and moves below the line segment AH 2 , enter the region (cf. Fig. 8a). The points that satisfy L F 0 H (Z ) > 0 lie on the right of the curve V S 0 (Z ) = 0 in the region R 1 (R 3 R 4 ) (cf. Fig. 8c).
Let B and C be the points of intersection of the curve V S 1 (Z ) = 0 with the switching lines 0 and 1 , respectively, where B = γ (L+σ )+δ ec , L + σ , C = (x L+σ , L + σ ) and On the switching line 1 , and so, the trajectories starting from 1 and evolve above the line segment C H 4 , enter the region (cf. Fig. 8b). The points that satisfy Fig. 8c).
The trajectories that enter into the region 2 with a transversal motion and satisfies Fig. 8a it is seen that the trajectories of (8) with initial points below 0 intersect the switching line 1 at a point on the line segment C D and remains within = 2 R 1 R 4 . The trajectories with an initial point above 1 intersect the switching line 0 at some point on the line segment AB and stays in thenceforth. This implies that the trajectories with a transversal motion enter of nonzero measure and stays there forever. Therefore, the region is a trapping region and so, a hysteresis loop around the predator population threshold L stabilizes system (8) through low-amplitude oscillations in the bounded region .
From Fig. 8c we observe that during the period of prey-refuge, the trajectories of system (8) stay in the region 2 and are directed toward the switching line 1 . Therefore, the time taken by the trajectories of the system (8) to reach the switching line 1 after hitting the switching line 0 can be regarded as the refuge time period of the prey. At this time, the predation pressure on the prey is quite low owing to the prey refuge protection. Similarly, the time taken by the trajectories of system (8) to reach the switching line 0 after hitting 1 refers to the non-refuge period. The non-refuge period refers to the time of over-exploitation of the prey by the predators. To determine the changes in the refuge time period owing to the changes in the threshold population density of the predators, we vary L subject to the condition y * 1 + σ < L < y * 0 − σ while all other parameter values remain unaltered. From Fig. 9a, b it follows that with an increase in L, there is a drop in the time period in the refuge. Therefore, with an increase in the predator population density, the more apprehensive prey is most likely to avail refuge protection, while the population of the bold prey will decrease due to the increased predation pressure (cf. Fig. 9c).
By increasing the parameter σ in the domain 0 < σ < |L − y * | and keeping all other parameters of the system fixed, from Fig. 10a, b we observe that the prey prefers to take refuge protection for a longer time period. This implies that the relatively less vigilant prey will stay in refuge protection for a longer period even when the predator population density drops significantly and will stay in the open habitat for a longer period even when the predator population becomes high. Frequent switching of habitat of the prey from the refuge to the open habitat happens when the hysteresis parameter value is low. In this case, the predator population density fluctuates from L − σ to L + σ . Therefore, hypervigilant prey species are most likely to move to and fro from the open habitat to the refuge depending on the changes in the population density of the predators.
To study the combined effect of L and σ on the refuge time period, we divide the Lσ -parameter plane into the following four regions: The regions I, II, III and IV in the σ − L parameter plane (cf. Fig. 11) represent the region of intermittent refuge protection, continuous refuge protection, no refuge protection, and compromised refuge, respectively. We know that the switching mode is possible only in the region I and the predator population fluctuates in between L − σ and L + σ . At a threshold population density L of the predators, the maximum value of the hysteresis parameter is σ L = Fig. 9 The time evolution of the population when the prey species is a apprehensive (L = 0.41) or b bold (L = 0.5), but with equally vigilant (σ = 0.03). c For σ = 0.03, the changes in the refuge stay time period with L as an active parameter. The shaded regions indicate the time spent by the prey in refuge min{L − y * 1 , y * 0 − L}. When the threshold predator population value is 1 2 (y * 0 + y * 1 ), the hysteresis parameter reaches its global maximum σ max = 1 2 (y * 0 − y * 1 ); and consequently, the refuge time period becomes maximum. In the region II, the predator population is greater than L + σ , whereas in the region III, the predator pop- c For L = 0.45, the changes in the refuge stay time period with σ as an active parameter. The shaded regions indicate the time spent by the prey in refuge ulation is less than L − σ . In the region IV, we have σ > σ max and so the least vigilant prey will continue to stay in the same habitat irrespective of the changes in the predator population density. In this case, the prey will be deprived of foraging or mating opportunities due to their continuous stay at the refuge or the prey Fig. 11 The time spent in refuge as a function of prey apprehension (L) and vigilance (σ ). The colour bar represents the refuge time period population gets severely compromised due to their continuous stay at the open habitat. It is observed that prey can be hypervigilant as well as bold (or apprehensive). In either case, the prey frequently switches their habitat due to the changes in the predator population. A bold and less vigilant prey continues to stay in the open habitat by compromising their stock (region III), whereas a more apprehensive and less vigilant prey prefer staying at the refuge by compromising their foraging and mating opportunities in the open habitat (region II). We also observed that the more apprehensive prey (in the region I) prefers to stay for a longer time at the refuge even if they are hypervigilant.

Discussion
In this paper, we have studied a non-smooth system with density-dependent intermittent refuge protection of the prey by employing Filippov's convex method and Utkin's equivalent control method. The qualitative analyses for the subsystems of the Filippov system were carried out, confirming the non-existence of periodic solutions of the subsystems and giving the threshold conditions for the existence of transcritical bifur-cations. We obtained the conditions for the existence of the regular and virtual equilibrium, boundary equilibrium, tangent points and pseudo-equilibrium of the Filippov system. We verified the existence of a pair of pseudo-equilibria when the level of apprehension of the prey is less than some threshold value. Under this level of prey apprehension, we have shown that one of the pseudo-equilibria is a saddle point and the other one is a stable node. We investigated regular or virtual equilibrium bifurcation, boundary-node bifurcation and pseudo-saddle-node bifurcation. We studied the effect of hysteresis in the Filippov system by introducing switching delay arising due to different behavioural traits of the prey species. We identified the evolution of stable limit cycles around the threshold population density of the predators in some bounded region in the phase space due to the delay in the switching of habitat. From numerical simulations, we observed that a highly apprehensive prey has a longer refuge stay time compared to a bold prey. We also observed that a hypervigilant prey frequently switches their habitat compared to their less vigilant counterpart, whereas a less vigilant prey overstays either in the refuge or in the open habitat depending upon the population density of the predators. Irrespective of the level of apprehension of the prey, the hypervigilant prey frequently switches their habitat due to the changes in the predator population, decoupling prey vigilance and apprehension. We also observed that a bold and less vigilant prey continues to stay in the open habitat by compromising their stock, whereas a more apprehensive and less vigilant prey prefers staying at the refuge by compromising their foraging and mating opportunities in the open habitat.