A discrete-time dynamical model of prey and stage-structured predator with juvenile hunting incorporating negative effects of prey refuge

This paper examines a discrete predator-prey model that incorporates prey refuge and its detrimental impact on the growth of the prey population. Age structure is taken into account for predator species. Furthermore, juvenile hunting as well as prey counter-attack are also considered. This paper provides a comprehensive analysis of the existence and stability conditions pertaining to all possible fixed points. The analytical and numerical investigation into the occurrence of different bifurcations, such as the Neimark-Sacker bifurcation and period-doubling bifurcation, in relation to various parameters is discussed. The impact of the parameters reflecting prey growth and prey refuge is thoroughly addressed. Numerous numerical simulations are presented in order to validate the theoretical findings.


Introduction
One of the most pivotal mechanism in maintaining the ecological balance of the ecosystem is the dynamic interplay between prey and predator.In recent years, mathematical models have gained significant traction and utility in explicating population dynamics.The predator-prey model has garnered significant attention from researchers in the field of ecology, following the groundbreaking contributions of Lotka [1] and Volterra [2] to the field.Subsequent to that, numerous kinds of enhancements for the predator-prey model have been suggested [3][4][5][6][7][8].
Predominantly, within predator-prey systems, it is commonly assumed that predators within a given population possess uniform predation capacity and fecundity.However, in the real world, several authors believed that predators residing within a specific population can be classed by two fixed ages: juvenile predator and adult predator.Numerous scholarly articles have been dedicated to examining the dynamics of populations structured by stages.In recent times, several researchers have employed stage structure in prey species [9,10], while numerous scholars have utilised stage structure in predator species [11][12][13] as well.This paper will focus solely on the stage structure of predator species.Morever, during this juvenile phase, predators develop the essential predatory skills required for their survival.
Predation during the juvenile phase poses a significant challenge, given that early juvenile predators lack the requisite abilities and expertise in foraging and hunting.In numerous instances, when juvenile predators partake in attacks on their prey, the prey species respond by initiating counter-attacks, leading to the killing of young and inexperienced predators as a means of self-defense [14][15][16].Only a small number of scholarly articles have explored the intricacies of juvenile hunting through the application of mathematical models [17][18][19].So, juvenile hunting and in response counter-attacks by the prey are also addressed in this work.
Numerous studies [20,21] have confirmed that predators employ diverse tactics in order to capture their prey .Likewise, prey species employ diverse strategies [22][23][24] to mitigate the rate of predation.Prey refuge is one of them.The utilisation of refuges can afford a certain degree of protection to prey species.Despite its benefits, prey refuge can have negative consequences on prey's growth.The utilisation of refuges by prey carries significant costs, particularly in terms of potential reductions in feeding or mating success due to increased time spent in refuges [25].As a result, the presence of prey refuge may result in a decrease in growth of prey [26][27][28].In this paper, we consider that prey uses refuge, which has a detrimental effect on their growth.
In cases where there are populations with overlapping generations, the birth processes take place in a continuous manner.As a result, the interaction between predator and prey is typically represented through the use of ordinary differential equations [29].However, it is worth noting that in reality, there are other types of species, such as monocarpic plants and semelparous animals [30], that exhibit discrete non-overlapping generations and their births occurexclusively during regular breeding seasons.Their interactions are characterised by difference equations or as discrete-time mappings.The dynamics of discrete-time predator-prey models might reveal greater complexity compared to their continuous-time models [31].When the population size is relatively small, it is appropriate to use discrete models to represent populations, even if some species have a long lifespan and overlapping generations.Additionally, population change is typically examined on a yearly (or monthly, or daily) basis.Hence, it is imperative to examine discrete population dynamical models [32].In recent years, there has been a significant increase in collaboration among ecologists studying discrete dynamical ecological models [33,34].
Kaushik et al. [17] considered a mathematical model which is as follows: where x, y, z represent population sizes of prey, juvenile predator, and adult predator respectively.The objective of this study is to analyse the aforementioned model, taking into account the prey's refuge behaviour and its detrimental impact on prey population growth, within the context of a discretetime framework.This analysis aims to provide an in-depth investigation of the advantages and disadvantages associated with the prey refuge, a topic that has not been thoroughly explored in the existing scholarly literature.To the best of the authors' knowledge, there has been no prior investigation into a mathematical model that elucidates the adverse effects of prey refuge on the population dynamics of prey.
This paper is structured in the subsequent fashion: section 2 presents the mathematical framework of the system.Section 3 presents the parametric conditions that pertain to the existence and stability of the equilibrium points.The theorems pertaining to Neimark-Sacker bifurcation and period-doubling bifurcation are presented in sections 4 and 5.In section 6, Numerical simulations are provided to validate the analytical results.Finally, section 7 draws a quick conclusion.

Mathematical Modelling
Kaushik et al. [17] examined a mathematical model that describes the dynamics of a prey species and a predator species with stage structure, as presented in equation (1).It is postulated that both juvenile and adult predators exhibit Holling type I functional response when interacting with prey species.Morever, it is considered that juvenile predators cannot reproduce, only adult predators can.The model ( 1) is further modified to incorporate the anti-predator effect, specifically the utilisation of prey refuge.We assume that prey uses prey refuge to mitigate adult predator attacks.The prey refuge effect has no impact on juvenile hunting.This may be due to the fact that the size structure or desire to hunt causes juveniles to exert extra effort to thwart the prey refuge effect, or because the prey species has no fear response to juvenile hunting, and thus the anti-predator behaviour towards juvenile predators is not the prey refuge but the aggressive counter-attack.
Let n represents the prey refuge constant, such that nx represents the number of prey species that are unaccessible to the adult predator, and adult predators do not concern themselves with the pursuit of this quantity of prey.Consequently, (1 − n)x denotes the quantity of prey that are available for consumption by the adult predators.We employ a negative effect of this strategy in terms of diminished growth.We assume that the unavailable preys do not take part in the growth of prey population.Therefore, the total growth of the prey population at any instant of time is r(1 − n)x(1 − x/k).We also assume that ϕ = 0 in the model (1) , therefore, the modified system of equation becomes We discretize (2) by Euler's forward method, and obtain the discretized model Considering here, the variables x t , y t , and z t denote the population sizes of prey, juvenile predator, and adult predator at generation t, where t ∈ N. r indicates the prey's growth rate, and k represents the system's environmental carrying capacity.The predation rates of juvenile and adult predators are denoted by α 1 and α 2 respectively, µ is the conversion efficiency or reproduction rate of the adult predators, α 3 is the prey counter-attacking rate to juvenile predators, β is the juvenile predators' natural death rate, γ is the juvenile predators' maturation rate, and m is the adult predators' depletion rate in the absence of prey, and n ∈ (0, 1) is the coefficient of prey refuge.

Equilibrium points and their stability
This section discusses the existence and stability of all biologically viable equilibrium points.After performing some calculations, all of the equilibrium points that are biologically feasible have been determined.These are vanishing equilibrium point E 1 (0, 0, 0), axial equilibrium point E 2 (k, 0, 0) and the coexisting equilibrium point E 3 = (x * , y * , z * ),where, The vanishing equilibrium is E 1 (0, 0, 0).E 1 exists for all biologically possible parameter values.It is unstable in nature, as proven by the following theorem.
Theorem 1.The vanishing equilibrium E 1 is not stable.
Proof.The eigenvalues of the Jacobian matrix at E 1 (0, 0, 0) are given by and m < γ.Hence proved.

Axial equilibrium(E 2 ):
The axial equilibrium point is given by E 2 (k, 0, 0).Clearly E 2 exists for all possible parameter values of the system (4).The following theorem demonstrates that, under certain parametric conditions, the axial equilibrium point E 2 exhibits stability.
Proof.The Jacobian matrix of the model (( 4)) at the axial equilibrium point E 2 (k, 0, 0) is Now, the eigenvalues of the Jacobian matrix The stability of the fixed point E 2 (k, 0, 0) is reliant on the absolute values of the eigenvalues of the Jacobian matrix evaluated at the axial equilibrium point E 2 .The axial equilibrium
The existence of the coexisting fixed point E 3 is evident from figure (1a).The x-nullclines, ynullclines, and z-nullclines of the model (4) are depicted in figure (1a), with the x-nullclines shown in brown, the y-nullclines in blue, and the z-nullclines in red.The parameter values used for this illustration are provided in table (1).The coexisting fixed point E 3 (0.232013, 0.0154654, 0.0419775) represents the intersection point of these nullclines.The Jacobian matrix of the system (4) at any point (x, y, z) is given by Now, the Jacobian matrix of the system (4) at the interior equilibrium point E 3 = (x * , y * , z * ) is The characteristic equation of the matrix J e3 is as follows: where, , and The subsequent theorem provides proof for the stability of the fixed point, E 3 = (x * , y * , z * ).
Theorem 3. The coexisting equilibrium E 3 is locally stable if and only if 3 , and Proof.Please refer to Theorem 3.2 in [35].

Neimark-Sacker bifurcation
The Neimark-Sacker bifurcation is a well-known bifurcation phenomenon that occurs in dynamical systems when a stable limit cycle experiences a loss of stability, leading to the emergence of an invariant torus or periodic cycles.In order to analyse the Neimark-Sacker bifurcation phenomenon concerning the coexisting equilibrium E 3 , it is necessary to utilise the explicit criterion [36] outlined below.
Theorem 4. [36] Given a discrete dynamical system of l dimensions: , where v ∈ R denotes a bifurcation parameter.Assume that Z * is a fixed point of f v .Then, the characteristic polynomical for Jacobian matrix J(Z * ) = (a ij ) l×l of l-dimensional map f v is as follows: where, b i = b i (v, c),i=1,2,3,...l, c represents either a control parameter or another parameter that requires determination.Let us consider, a sequence of determinants of the type ,where Furthermore, it is assumed that the subsequent criteria are true: for i=l-3,l-5,...,2 (or 1), when l is odd or even, respectively.
The following theorem offers criteria that establish the occurrence of Neimark-Sacker bifurcation for system (4) with respect to the bifurcation parameter r.
Theorem 5.The fixed point E 3 experiences Neimark-Sacker bifurcation at the critical value r = r ns depends on the satisfaction of the specified conditions.
Proof.Let us consider, r as a bifurcation parameter and l = 3.Now, following the theorem (4) and using the equation ( 5), we compute the following values .
Other parameters can also be taken into consideration as the bifurcation parameter, leading to similar results.

Period-doubling bifurcation
The period doubling bifurcation is a notable occurrence in discrete dynamical systems, wherein the system experiences a series of bifurcations that lead to the doubling of the period of its orbits.To conduct an analysis of the Period-doubling bifurcation for the map (4) about the fixed point E 3 , a specific criteria [37] is required, as described in the following section.Theorem 6.Given a discrete dynamical system of l dimensions: Z u+1 = f r (Z u ), where r ∈ R denotes a bifurcation parameter.Assume that Z * is a fixed point of f r .Then, the characteristic polynomical for Jacobian matrix J(Z * ) = (a ij ) l×l of l-dimensional map f r is as follows: where, b i = b i (r), i=1,2,3,...l .Let us consider, a sequence of determinants of the type (∆ ± i (r)) l i=0 such that ∆ ± 0 (r) = 1, and , where K 1 and K 2 are same as given in theorem (4).Then, a period-doubling bifurcation occurs at a critical value r = r pb if and only if the following requirements are fulfilled.(i) Eigenvalue requirement: P r pb (−1) = 0, P r pb (1) > 0 , ∆ ± l−1 (r pb ) > 0, ∆ ± i (r pb ) > 0, i=l-2, l-4,...., 2(or 1), when n is even or odd, repectively.
(ii) Transversality requirement: By employing the aforementioned theorem, we determine the conditions that lead to the occurrence of period-doubling bifurcation in relation to the parameter r.Theorem 7. [35] The fixed point E 3 of the map (4) exhibits a period-doubling bifurcation at r = r pb when the subsequent specified conditions are satisfied.
where, the values of p 1 , p 2 , and p 3 are provided in equation ( 5).  6 Numerical simulation In this section, we present numerical simulations to validate the theoretical findings previously discussed in the preceding sections.The hypothetical parameter values depicted in table (1) are taken into consideration.The Mathematica software is employed for conducting numerical simulations to facilitate the analysis of the obtained results.At first we consider the parameter values h = 0.1, r = 0.5, n = 0.01, α 1 = 10, α 2 = 10, µ = 0.275944, α 3 = 0.03, β = 5.3, γ = 9.5, k = 100, and m = 3.5.In order to validate the stability requirements of the axial equilibrium point E 2 as given in theorem (2), these parameter values are employed.By utilising the given parameter values, the eigenvalues of the Jacobian matrix J e2 can be determined.These eigenvalues are |λ 4 | = 0.9505 < 1, |λ 5 | = 0.678371 < 1, and |λ 5 | = 0.845371 < 1.As a result, as illustrated in the figure (2), the fixed point E 2 is stable.We now take the parameter values listed in the following table (1)  Using these parameter values, we compute the characteristic equation of the Jacobian matrix J e 3 which is given by Comparing equation ( 6) with equation ( 5 In order to validate the outcome presented in theorem (5), we examine the parameter values r ∈ (0.2, 0.7) , µ = 3.07227, while keeping the remaining parameters consistent with those specified in table (1).Parameter r is used as the bifurcation parameter in this case.In the vicinity of the parameter value r = 0.539 = r ns , the fixed point (0.179337, 0.0118774, 0.0322387) undergoes a transition in stability, transitioning from a stable fixed population to a stable periodic population as a result of a Neimark-Sacker bifurcation.For a given value of r = 0.539, µ = 3.07227, and assuming all other parameters are as specified in table (1), we find the characteristic equation of the Jacobian matrix at the fixed point here, p 1 = −1.15989,p 2 = −0.645119,and p 3 = 0.827696.Now, we find 1 )) r=r ns = 0.0000654356 ̸ = 0, and using the equation cos( 2π j ) = 0.993794, one obtains j = ±56.3685,therefore, the non-resonance criterion is also satisfied i.e, all the necessary conditions for the occurrence of the Neimark-Sacker bifurcation have been satisfied, as stated in theorem (5).The visual representations for the same can be observed in the diagrams depicted in Figure (4).Furthermore, by selecting a value of r that is less than r ns , specifically r = 0.48, and µ = 3.07227, while keeping all other parameters as specified in table (1), it is found that all the eigenvalues of the Jacobian matrix J e3 are λ 7 = 0.994474+0.104977i,λ 8 = 0.994474−0.104977i,and λ 9 = −0.828008i.e., |λ 7 | < 1, |λ 8 | < 1, and |λ 9 | < 1.This confirms the stability of the fixed point E 3 .Although, with a value of r = 0.6 > r ns and all other parameter values remaining the same as previously stated, the eigenvalues of J e3 are found to be λ 7 = 0.998289 + 0.105108i, λ 8 = 0.998289 − 0.105108i,   To examine the conditions for the occurrence of period-doubling bifurcation, as laid out in theorem (7), we consider the parameter values r ∈ (22,25) , µ = 0.59977, and the remaining parameters are maintained in accordance with the values provided in the table (1).Here , r is taken as the bifurcation parameter.The stability of the coexisting fixed point (0.920016, 0.0509269, 0.13823) undergoes a transition from stability to instability at the value r = 23.7137= r pd , resulting from a period-doubling bifurcation.At the period-doubling bifurcation point r = r pd , the fixed point E 3 undergoes destabilisation, resulting in the emergence of two points that constitute the period-2 solution.The characteristic polynomial of J e3 with r = r pd , µ = 0.59977 and the other parameters as stated previously, λ 3 + 0.992639λ 2 − 0.951366λ − 0.944005 = 0 (8) here, p 1 = 0.992639, p 2 = −0.951366,and p 3 = −0.944005.
Next, we proceed with the computation of the expression 1 − p 2 + p 3 (p 1 − p 3 ) = 0.123164 > 0, 1 + p 2 − p 3 (p 1 + p 3 ) = 0.094544 > 0, 1 + p 2 = 0.0486336 > 0, 1 − p 2 = 1.95137 > 0, 1 + p 1 + p 2 + p 3 = 0.0972673 > 0, and −1 + p 1 − p 2 + p 3 = 0 which implies the fact that, as stated in theorem ( 7), all the requirements for a period-doubling bifurcation are met in the vicinity of the coexisting fixed point (0.920016, 0.0509269, 0.13823) at the critical value of the bifurcation parameter r = r pd .The figure (5) illustrates the period-doubling bifurcation diagram associated with the parameter r.Furthermore, when r = 22 < r pd and all other parameter values remain unchanged as previously discussed, the eigenvalues of the matrix J e3 are found to be 0.975264, −0.90591 + 0.124984i, and −0.90591 − 0.124984i.These eigenvalues have modulus less than 1, indicating that the fixed point E 3 is stable.However, considering the value of r is 25, which is greater than r pd , and assuming that the remaining parameter values are the same as those discussed earlier, the eigenvalues of the Jacobian matrix J e3 are determined to be | − 1.18461| > 1, |0.975302| < 1, and | − 0.900491| < 1, thus confirming the unstable nature of the coexisting fixed point E 3 .
Figure 5: Period-doubling bifurcation is depicted in relation to the bifurcation parameter r.The figure is constructed using the parameters r = 23.7137 and µ = 0.59977, with the remaining parameters are kept consistent with the values specified in table (1).Now, considering all the parameter values specified in table (1), with the exception of µ = 3.3125, and by varying the prey refuge parameter n, it becomes apparent that the stability of the coexisting fixed point experiences a change near the value n = 0.0555353 = n ns .More precisely, stable periodic cycles arise as a result of the manifestation of a Neimark-Sacker bifurcation at the critical value n = n ns .Numerical verification of this claim can be carried out by applying the theorem (5).After some calculations, we get p 1 = −0.969343,p 2 = −0.985856,and p 3 = 0.983486.Subsequently, the following calculations are obtained: 1 )) n=n ns = 0.114834 ̸ = 0, and by utilising the equation cos( 2π j ) = 0.988932, it is found that j = ±42.1926,consequently, it can be concluded that the non-resonance criterion is also met.Therefore, based on theorem (5), it can be concluded that all the necessary conditions for the occurrence of a Neimark-Sacker bifurcation have been satisfied.In addition, the figure (6) confirms the same.

Conclusion
The primary objective of our study is to investigate a predator-prey model incorporating stage structure with juvenile hunting, paying special attention to the adverse effects of prey refuge behaviour on their own population dynamics, utilising a discrete-time mathematical model that accurately represents this specific scenario, taking into account nonoverlapping generation for the species under consideration.The model under investigation in this paper is a discrete counterpart of the continuous model proposed by Kaushik et al. [17] with certain modifications made to the original model.In this paper, the existence conditions of all ecologically relevant fixed points are found and a stability analysis of these fixed points are done.Under certain parametric conditions, it is observed that both the axial fixed point and the coexisting fixed point exhibit stability.A comprehensive bifurcation analysis is performed.It is found that the population dynamics in this model are significantly influenced by the intrinsic growth rate of the prey (r).Different bifurcations of codimension 1, such as the Neimark-Sacker bifurcation and the period-doubling bifurcation, can be observed when the growth rate of the prey-related parameter r is varied.The figures (4) and ( 5) illustrate the manifestation of the Neimark-Sacker bifurcation and the period-doubling bifurcation, respectively, in relation to the parameter r.Furthermore, the parameter n, which is related to the prey refuge, plays a crucial role in sustaining population stability within the system under study.It is observed that varying the value of the parameter n can result in the coexistence of all species or the extermination of predator species.The occurrence of a Neimark-Sacker bifurcation is also observed, leading to the destabilisation of the system when the value of the prey refuge parameter varies, as portrayed in the figure (6).These illustrate the significance of prey refuge in the system under consideration.Various numerical simulations and graphical representations are presented in this paper to demonstrate the intricate dynamics of the system model.
(a) Existence of the coexisting fixed point E 3 , (b) Stability region of the coexisting fixed point E 3 in nrµ-space

Figure 1 :
Figure 1: The existence and stability of the coexisting fixed point E 3 are displayed using the parameter values stated in table (1).
where b ′ i represents the first derivative of b i with respect to r at r = r pb .

( a )
Time series of prey (b) Time series of juvenile predator (c) Time series of adult predator

( a )
Time series of Prey (b) Time series of juvenile predator (c) Time series of adult predator

Figure 3 :
Figure3: The stability of the coexisting fixed point E 3 is demonstrated using parameter values from table(1)

( a )
Time series of prey (b) Time series of juvenile predator (c) Time series of adult predator (d) Phase portrait

Figure 4 :
Figure 4: The occurrence of the Neimark-Sacker bifurcation is demonstrated when the parameter r is varied.These figures have been generated utilising the parameter values r = 0.539 and µ = 3.07227, while the remaining parameter values are sourced from table (1).

( a )
Time series of prey (b) Time series of juvenile predator (c) Time series of adult predator (d) Phase portrait

Figure 6 :
Figure 6: The occurrence of Neimark-Sacker bifurcation is shown when the prey refuge parameter n is altered.The figure is generated utilising the parameter values n = 0.055535 and mu = 3.3125, while the remaining parameters are maintained in accordance with the values indicated in table (1).

Table 1 :
(3)validate the stability criteria of the coexisting equilibrium E 3 as mentioned in the theorem(3).Parameter values of the system (4) for the purpose of numerical simulation