Stability and bifurcation in a predator-prey system with effect of fear and additional food

1.   Introduction

Predation is one of the most common phenomena for survival in the wilderness that connects all the elements of a multispecies food web. It plays an essential role in maintaining stability and balance in the nature. Mathematical models are the most effective tools to illustrate the relationship between two or more species. Initially, the mathematical models for prey-predator interactions were explored by two mathematicians, Alfred J. Lotla ^[1] and Vito Volterra ^[2], independently in the mid 1920s. They found periodic fluctuations of species around their steady states. Since then, ecologists and mathematical biologists have been in practice with various models considering different essential parameters ^[3,4,5,6]. The dynamics of the prey-predator system are significantly influenced and become more realistic by the incorporation of some notions such as functional response ^{[7,8,9,10,11]}, stage structure ^[12,13], effect of fear ^[14,15] and time delay ^[5].

Understanding the complex dynamics of a population model depends on the functional response of predators, which quantifies the prey captured by each predator during certain periods of time. It significantly influences the dynamical features of the system. From an ecological perspective, incorporating the functional response enhances the realism and accuracy of the model. In recent years, numerous studies have been conducted to explore the effects of different functional responses on the prey-predator system. These studies include the Holling type-II ^[7,16,17], Holling type-III ^[8,18], Holling type-IV ^[6,12,19], Beddington-DeAngelis ^[9,20,21,22] and Crowley-Martin functional responses ^{[10,11,13,23,24,25,26,27]}.

In recent decades, several researchers have deeply investigated the qualitative behavior of prey-predator systems when supplementary nourishment sources are available for the predators ^{[17,20,28,29,30,31,32]}. In particular, Sahoo ^[32] put forward and explored a model to investigate the significance of supplementary nourishment within an eco-epidemiological system. The research focused on managing chaotic population dynamics of predator-prey interactions by introducing supplementary nourishment resources for top predators. To accomplish this, they developed a three-species predator-prey model that specifically incorporates the provision of supplementary food to the top predators. In another study, Thirthar et al. ^[20] focused on prey, predator and super predator within a predator-prey system. In this study, predation interactions between prey and predator are regulated by the Beddington-DeAngelis functional response, and the impact of predation fear in prey species was also carried out. The authors assumed that the predator species offers protection to the super predator species. It was anticipated that the super predators, driven by fear of predators and their refuge behavior, would receive a constant supply of additional food, which in turn would benefit both the super predators and the predators. Furthermore, the latest research by Kumar and Dubey ^[17] examined the intricacies of a prey-predator model with the effect of supplementary nourishment, multiple delays, and the interaction between the prey and predator was modeled using a Holling type-II response. The provision of extra nourishment was specifically targeted toward the predator population, aiming to reduce its dependence on the prey species. Srinivasu et al. ^[30] examine a modified predator-prey model with the type II functional response to study how providing supplementary nourishment to the predator impacts the behavior of the system. The research emphasizes that the handling times of the available foods significantly influence the eventual state of the ecosystem. Ghosh et al. ^[31] investigate the influence of providing supplementary nourishment to predators within an ecological model that includes a prey refuge. The authors illustrate theoretically and numerically that species can co-exist in oscillatory patterns when there is an availability of both abundant and high-quality supplementary nourishment and prevent predator extinction in systems with high prey refuge, offering insights for conservation biologists in real-world ecological systems.

Since predation events are obvious, our knowledge of the prey-predator system in forest provinces has mostly concentrated on the direct killing of predators on prey populations. The effects of predators on their prey can be physical as well as mental. Some recent studies have revealed a new perspective, showing that predators not only reduce prey populations by directly killing prey but also by inducing fear in prey individuals, which results in a fall in the prey birth rate. Zanette et al. ^[33] carried out research involving song sparrows, manipulating predation risk by playing predator sounds to some species and non-predator sounds to others. Despite protecting nests from direct killing, the fear of predators alone led to a significant $40\%$ reduction in offspring production due to decreased egg laying, lower hatching success and increased nestling mortality. Various anti-predator responses contributed to these demographic effects. In another study by Cresswell ^[34], the impact of predators on population size extends beyond direct lethal effects, as fear of predation leads to costly behavioral changes with potentially significant population and community effects. The dominance of lethal or non-lethal effects depends on factors such as resource constraints, prey vulnerability, trophic interactions, and evolutionary history, making cases without any population effects due to predation or cases dominated solely by lethal effects rare and specialized.

As a result of the fear induced by predators, prey individuals experience a restricted environment for their activities, including mating, preventing them from freely accessing open habitats. The fear inhibits them from venturing into open spaces, which reduces their ability to reproduce ^[33]. We construct a function that represents the fear effect on the prey population. The reduction of birth rate in the prey population depends on the number of predator individuals $y$ and their nature. We quantify the nature (aggressiveness, agility, etc.) of the predator in terms of a parameter $k$, which is known as the cost of fear. Thus, the fear function $f(y, k)$ is a function of $y$ and $k$, which follows the following characteristics:

$(1)$ In the absence of predators or when the cost of fear is negligible, the impact on the reproduction rate of the prey species remains same

$(2)$ The reproduction rate of prey approches zero at the high predator density or at the large cost of fear

$(3)$ As the density of predators and the cost of the induced fear increase, the reproductive rate of prey decreases

There may be several functions which follow all the conditions of a fear function, but commonly we use following functions due to their simplicity

The proposed prey-predator model incorporates the influence of fear on prey reproduction, suggesting that high levels of fear can contribute to stabilizing the system by eliminating the periodic solutions. In the study conducted by Wang et al. ^[14] on a terrestrial vertebrate, it was observed that fear of predators alone had a significant impact on anti-predator defense, leading to a notable reduction in prey reproduction. To capture the influence of fear on the population dynamics of terrestrial animals, the authors developed a predator-prey model that incorporates the reproductive cost associated with fear. Dubey et al. ^[35] analyze the predator-prey model that accounts for fear-induced responses in prey towards predators and anti-predator behavior exhibited by the prey. Additionally, the model incorporates delays in fear response and gestation. They state that the system experiences Hopf-bifurcation in the context of both delays.

Furthermore, Chen et al. ^[36] examined how the dynamics of predator-prey model behavior is affected by the fear effect and the Leslie-Gower function. When fear intensifies, the system will experience multiple dynamic behaviors, switching until the prey population is wiped out. However, because there are availability of other prey sources, the predator population manages to endure. The concept of effect of fear on the species dynamics has been utilized by various researchers to show the efficacy of their respective models ^{[15,37,38,39]}. Zhang et al. ^[38] examined anti-predator behavior and fear on a Holling type-II predator-prey model with prey refuge. The study shows that fear reduces the predator's population density and stabilizes the system by preventing the existence of periodic solutions, while the presence of prey refuge significantly affects the persistence of the predator population.

The Beddington-DeAngelis functional response, formulated by Beddington ^[21] and DeAngelis ^[22], introduces a predator-dependent functional response that incorporates predator interference. It shares similarities with the Holling type II response, as it considers the interactions between predators in addition to prey searching and processing activities. In this response, predators form groups consisting of two or more individuals that not only engage in prey-related tasks but also interact with each other. The Beddington-DeAngelis functional response suggests that as prey abundance increases, the impact of predator interference on the feeding rate diminishes, implying that the collective behavior of predators becomes less influential. In the study conducted by Crowley and Martin ^[23], it was anticipated that predation would decrease as a result of high predator density, which can lead to interference among individual predators. This interference occurs when predator individuals interact with each other and it can reduce the overall predation rate even when there is an abundance of prey available. The Crowley-Martin functional response describes the relationship between prey biomass and the population of birds. In the Crowley-Martin functional response, predator interference is taken into account, meaning that the presence of one predator can impede the activity of another predator. This interference affects both the prey and predator population. Many authors have used the Crowley-Martin functional response to demonstrate the effectiveness of their models ^{[11,13,21,22,24,25,26,27]}. The difference between the Beddington-DeAngelis and Crowley-Martin functional responses is how they consider the effects of predator interference. The Beddington-DeAngelis type response suggests that in the case of abundant prey, the impact of predator interference on the feeding rate becomes negligible. In contrast, the Crowley-Martin functional response consistently considers the impact of interference on the feeding rate, regardless of prey abundance. Experimental evidence supports the notion that mutual predator interference decreases the feeding rate of predators ^[40,41].

Maiti et al. ^[23] investigated a prey-predator model that includes prey refuge and a Crowley-Martin type response to describe the interaction between unreserved prey and predators. The study focuses on analyzing the dynamic of the system while considering delay as a parameter that can cause bifurcation. Tripathi et al. ^[26] examined a two-dimensional predator-prey model with discrete delay and the Crowley-Martin functional response, which incorporates mutual interference by predators at high prey density. The authors analyze the impact of predator density dependence on system permanence and numerically describe Hopf-bifurcation in the non-delayed system. Moreover, Papanikolaou et al. ^[27] found that the Crowley-Martin model outperforms other functional responses, including the Beddington-DeAngelis, Hassell-Varley and Holling type-II functional responses. This finding further highlights the effectiveness and relevance of the Crowley-Martin functional response in capturing the dynamics of predator-prey interactions. Considering these insights, investigating the Crowley-Martin type interaction between prey and predator can potentially reveal rich and complex dynamics within the system.

As we mentioned earlier, fear is a psychological activity of prey in the presence of predators, which was ignored in classical models. The prey shows anti-predator nature under fear, such as reducing the foraging time, changing the living locality, adjusting reproduction strategies, etc. Due to all these, the reproduction rate decreases. Therefore, it is reasonable to consider fear in the model. The additional food is provided to predators to avoid the overexploitation of prey resources and maintain the ecological balance. In the existing models, the additional food is provided at a constant rate, but it must vary with the number of predators. In the present model, we introduce the additional food in a separate equation and keep its quantity depending on the number of predators. Since the effect of fear varies with the dependency of the predator on the provided food, we incorporated a new parameter, the degree of dependency. To the best of our understanding a mathematical model with all the above points has not been investigated. Keeping all these in view, in Section 2, we discuss the model construction. In Section 3, we show our model is biologically well-behaved by proving the positivity and boundedness of the solutions. We discuss various dynamical features of our system in Section 4. In Section 5, an extensive computer simulation is carried out to validate our results. The study concludes with a few final remarks.

2.   Model formulation

We consider an ecological habitat where two species survive, one is prey and another is predator. The prey population follows the logistic growth behavior when predators are absent. This is governed by the equation below:

where $x$ is the number of individuals in the prey population at time $t$, $r$ and $r_{0}$ are the birth rate and the natural mortality rate of prey and $r_{1}$ is the intra-specific interference among the prey individuals. In the latest research, it was discoverd that the presence of predators induces fear within the prey population which leads to a reduction in their reproduction rate. This fear-induced reduction in reproduction rate is influenced by both the predator population $y$ and a parameter $k$ representing the cost of fear. Consequently, the modified form of the prey population growth equation is as stated:

where $f(y, k)$ is referred as fear function and $rf(y, k)$ is the effective growth rate of prey under the fear effect. In this study, we consider the influence of predation fear on the prey population when predators are present. The cost of anti-predator defense is related to the predator population density. To account for this, we replace the prey reproduction term $rx$ by multiplying it with the fear function $f(y, k) = \dfrac{1}{1+ky}$, emphasizing the costs carried on by predation fear. This modification allows us to explore how predation fear affects the dynamics of the prey-predator system in the presence of predators. The following system governs the prey-predator interaction in presence of fear effect.

where $Q(x, y)$ is functional response of predator and $c_{1}$ $(0 \le c_{1} \le 1)$ is the conversion efficiency of $y$ on $x$. The parameters $\delta$ and $e$ are the natural death rate of predators and the intra-specific interference among predators.

The Holling type-II response is prey-dependent, indicating that the predator feeding rate declines with increasing prey density due to handling and searching time. The Crowley-Martin response characterizes the prey-predator interaction, predicting a decrease in predator density at high predator concentrations due to resource competition. The Crowley-Martin response is formulated as

where $\alpha$ is the capture rate, $a$ is the handling time and $b$ is the degree of interference among predator individuals. Additionally, it changes into mass action law (when $a$ and $b$ are zero) and Holling type-II response (when $a > 0, b = 0$ or $a = 0, b > 0$). With the Crowley-Martin functional response, our model implies the forthcoming form.

Additional food in a prey-predator system refers to providing extra nourishment resources to the predator population, potentially beyond what is naturally available from their primary prey. This additional food can have various effects on the dynamics of predator-prey interactions and may influence population sizes, stability and coexistence patterns in the ecosystem. As a result, the equation describing the relationship between the additional food resource and the predator population has been modified as follows:

where $A$ is the quantity of additional food provided to the predator at time $t$. The parameters $\lambda$ and $\beta$ are the maximum supply rate and the natural depletion rate of additional food, respectively. The parameter $\phi$ is the consumption rate of additional food by predators and $c_{2}$ ($0\leq c_{2} \leq 1$) is the conversion efficiency of $y$ on $A$. The three-dimensional predator-prey model with additional food can be expressed as follows:

Further, we introduce the degree of dependency in the prey-predator system, which refers to the reliance or dependence of the predator population on its primary prey as a food source. The parameter $A_{0}$ ($0\leq A_{0} \leq 1$) is the degree of dependency of predators on prey or the dependency factor. If $A_{0} = 1$, predators are fully dependent on prey and no additional food is necessary. If $A_{0} = 0$, predators are completely dependent on supplementary nourishment, and the prey species increases logistically. If $A_{0}$ is in between $0$ and $1$, then the predator partially depends on prey and supplementary nourishment. The predator depends partially on the prey and the degree of dependency $A_{0}$, which depends on the nature of the prey and provided additional food. Therefore, in the presence of additional food, the fear function $f(y, k)$ also depends on the degree of dependency $A_{0}$. To modify the fear function $f(y, k)$, we consider the following points:

$(1)$ When the predator completely depends on the additional food then no fear will be induced by the predator, i.e., when $A_{0} = 0$ then $f(y, k, A_{0}) = 1$.

$(2)$ When the predators are fully dependent on prey then induced fear is maximum, i.e., $A_{0} = 1$ then $f(y, k, A_{0})$ is maximum and same as the fear function discussed in the introduction.

$(3)$ The amount of fear monotonically increases with respect to $A_{0}$, i.e., $\dfrac{\partial f}{\partial A_{0}} > 0$.

These points lead us to the following modified fear function, $f(y, k, A_{0}) = \dfrac{1}{1+kA_{0}y}$.

Considering the points mentioned above, our mathematical model can be described as follows:

The biological description of the parameters in the system (2.3), including their units and values, is provided in Table 1.

Remark 1. From the first equation of the above system, $\dfrac{dx}{dt}\leq (r-r_{0})x$, which implies $\lim_{t\to\infty} x(t) = 0$ if $r < r_{0}$ and consequently, $\lim_{t\to\infty} y(t) = 0$ at low supply rate of additional food. This is a trivial case from a biological perspective, to avoid this case we consider $r > r_{0}$ throughtout this paper.

3.   Boundedness and persistence of the system

A biological model must carry some essential characteristics such as positivity and boundedness of the solutions, persistence of the species. The existence of the species is demonstrated by the positivity and the boundedness of the solutions represents a natural restriction on species to grow exponentially. Both these characteristics are followed by the system (2.3) and have been shown by the theorem below.

Theorem 1. The set

where $\sigma = \min\{r_{0}, \beta, \delta\}$, is a positive invariant set for all solutions of the system (2.3) that start in the interior of the positive quadrant.

Proof. We write the system (2.3) in matrix form

where $X = (x, A, y)^T \in R^3$ and $G(X)$ is given by

The fundamental theorem of ordinary differential equation ensures the local existence and uniqueness of the solution because $X(0) = X_{0} \in R_{+}^3$ and $G: R^3 \rightarrow R_{+}^3$ are locally Lipschitz-continuous in $\Omega$. It follows that $X(t) \geq 0$ for all $t \geq 0$ because $[G_{i}(X)]_{x_{i}(t) = 0, X \in R_{+}^3} \geq 0$. In fact, from the first and third equations of system (2.3), it can easily be seen that $\dot x|_{x = 0} \geq 0$, $\dot y|_{y = 0} \geq 0$ and hence $x(t) \geq 0$, $y(t) \geq 0$ for all $t \geq 0$, second, $\dot A|_{A = 0} = \lambda (1-A_{0})y\geq 0$ for all $t \geq 0$ (as $y(t) \geq 0$ for all $t \geq 0$) and hence $A(t) \geq 0$ for all $t \geq 0$. Therefore, for system (2.3) $x(t), A(t), y(t)$ all are positive.

Now, we express the first equation of system (2.3) as follows:

which implies

Now, we construct a function

which on differentiation

where $\sigma = \min\{r_{0}, \beta, \delta \}$. Thus, it is evident that

and we observe that, if $x(t) \geq \dfrac{r}{r_{1}}$ and $W(t) \geq \dfrac{1}{\sigma}\Bigg [\dfrac{c_{1}r^2}{4r_{1}}+\dfrac{c_{2}^2 \lambda^2 (1-A_{0})^2}{4e}\Bigg]$, then $\dfrac{dx(t)}{dt} \leq 0$, $\dfrac{dW(t)}{dt} \leq\; 0$. This demonstrates that the system (2.3) solutions are bounded and remain in $\Omega$ for all $t > 0$ if $(x(0), A(0), y(0)) \in \Omega$. □

The persistence of a system indicates that all species will survive in the future and none will wipe out if they are initially present in the system. In the following theorem, we will establish some requirement for the uniform persistence of the system (2.3).

Theorem 2. The system (2.3) is uniformly persistent under the following conditions and is satisfied.

where $x_{i}$ and $y_{i}$ are defined in the proof.

Proof. If $M_{1}$ and $M_{2}$ are positive constants that satisfy each positive solution $X(t) = (x(t), A(t), y(t))$ of the system with positive initial conditions, then the system (2.3) is said to have uniform persistence.

Now, we define

As a result, Theorem 1 implies that

This further demonstrates that for any sufficiently small $\epsilon > 0$, there exists a $T > 0$ such that for all $t \geq T$, the following hold

For any $t \geq T$, we can write the first equation of the system (2.3).

consequently, it follows

In the case of arbitrarily small $\epsilon > 0$, we get

The prey species of the model (2.3) is persistent, if the inequality below holds true.

Now, the third equation of the model (2.3) implies

it follows that

In the case of arbitrarily small $\epsilon > 0$, we get

The predator species of the system (2.3) persists if the forthcoming conditions holds,

The second equation of the system (2.3) implies,

consequently, it follows

Again if $\epsilon$ is arbitrarily small then we have

Taking $M_{1} = \min\{x_{i}, A_{i}, y_{i}\}$, the theorem follows. □

4.   Dynamical analysis of the model (2.3)

4.1.   Equilibrium points

The equilibrium points are the steady-state solutions for a system. From the ecological aspects, their presence may determine whether the populations coexist, one species becomes extinct, or both species collapse. We investigated that, the system (2.3) may have four equilibrium points: trivial equilibrium $E_{0}(0, 0, 0)$, axial equilibrium $E_{1}(x_{*}, 0, 0)$, planar equilibrium $E_{2}(0, \tilde{A}, \tilde{y})$ and positive equilibrium $E^*(x^*, A^*, y^*)$. The equilibrium points $E_{0}(0, 0, 0)$ and $E_{1}(x_{*}, 0, 0)$ are always present in the system without any restiction, where $x_{*} = \dfrac{1}{r_{1}}(r-r_{0})$.

● Existence of $E_{2}(0, \tilde{A}, \tilde{y})$: $\tilde{A}$ and $\tilde{y}$ are positive solutions of the forthcoming system of non-linear equations

According to the second equation of the above system

When we plug the value of $A$ into the first equation of the system (4.1), we obtain

We observed that, the Eq (4.2) has two positive roots, if

hold.

System (2.3) has two prey-free equilibriums: $\tilde{E_{2}}(0, \tilde{A}, \tilde{y})$ and $\hat{E_{2}}(0, \hat{A}, \hat{y})$ under the conditions given in the Eq (4.3). Again, if $\phi c_{2} \lambda (1-A_{0}) < \phi \delta+\beta e$, the Eq (4.2) has no positive roots. As a result, in this case $E_{2}$ does not exist.

● Existence of positive equilibrium $E^*(x^*, A^*, y^*)$: It is evident that $x^*, A^*$ and $y^*$ are the positive solutions of the forthcoming system of non-linear equations:

We can obtain from the second equation of the system (4.4)

On substituting $A$ in the first and third equations of the system (4.4) the forthcoming equations are obtained:

From the Eq (4.5), we observe the following points:

$(1)$When $y = 0$, then

$(2)$When $x = 0$, then the quadratic equation follows

where

Equation (4.7) always has exactly one positive root (say $y_{0}$). $(3)$The slope of the curve is

It is evident that, $\dfrac{dy}{dx} < 0$ exists under the following inequlity

From the Eq (4.6), we observe the following points:

$(1)$When $y = 0$, then

$(2)$The slope of the Eq (4.6) at any point $(x, y)$ is

It is evident that, $\dfrac{dy}{dx} > 0$ under the following inequlity,

From the above analysis, it is concluded that the curve of the Eq (4.5) comes into the first quadrant through the $y$-axis intersecting at $(0, y_{0})$ and it is monotonically decreasing under the inequality (4.8) and then it leaves the first quadrant by intersecting the $x$-axis at $(x_{0}, 0)$. The curve of the Eq (4.6) enters in the first qudrant through the $x$-axis intersecting at $(x_{1}, 0)$. Under the inequality (4.9), the curve increases monotonically and approaches to its horizontal asymptote which is depited in the Figure 1 by the dashed line. Both the curves (4.5) and (4.6) intersect in the interior of first quadrant if

hold. From the above discussion, the system (2.3) has a unique positive solution $E^*$ if the inequalities (4.8) to (4.10) hold.

Consequently, we can formulate the following theorem.

Theorem 3. If (4.8)–(4.10) hold, the system (2.3) has a unique interior equilibrium, $E^*(x^*, A^*, y^*)$.

Remark 2. The number of positive equilibrium for the system (2.3) depends on values of parameters, which we have chosen. Several possibilities are depicted in Figure 2.

4.2.   Stability behavior

In ecology, an equilibrium point is said to be stable if a population settle and if the initial population is near to the equilibrium point. It is known that the local behavior of a system near an equilibrium point is the same as the corresponding linearized system. Hence, we investigate the stability behavior of the equilibrium points of the non-linear system (2.3).

Theorem 4. A necessary and sufficient condition for the negativity of the real parts of all the roots of the polynomial $L(\omega) = \omega_{1}^n+a_{1}\omega_{1}^{n-1}+.....+a_{n}$ with real coefficient is the positivity of all the principal diagonals of the minors of the Hurwitz matrix

For a proof of this theorem, see Bellman ^[42].

We calculate the variational matrix of the system and conclued the following point.

● The equilibrium point $E_{0}(0, 0, 0)$ is always saddle point having stable manifold along the $A$ and $y$ axes and unstable manifold along the $x$ axis.

● The equilibrium point $E_1(x_{*}, 0, 0)$ is

- locally asymptotically stable if $\delta > \dfrac{c_{1}\alpha A_{0}x_{*}}{(1+ax_{*})}$ hold.

- saddle point having stable manifold along the $A$ and $x$-axes and unstable manifold along the $y$-axis if $\delta < \dfrac{c_{1}\alpha A_{0}x_{*}}{(1+ax_{*})}$ hold.

- stable if $\delta = \dfrac{c_{1}\alpha A_{0}x_{*}}{(1+ax_{*})}$ hold.

● The Jacobian matrix examined at the prey-free equilibrium point $E_{2}(0, \tilde{A}, \tilde{y})$ can be obtained by

The characteristic equation of the above matrix can be obtained by

The roots of the Eq (4.11) have negative real parts, under the following conditions

Hence, using the Routh-Hurwitz criterion $\tilde{E_{2}}(0, \tilde{A}, \tilde{y})$ is locally asymptotically stable under the conditions in (4.12).

The roots of the Eq (4.11) have both positive and negative real parts if

Therefore, $\tilde{E_{2}}(0, \tilde{A}, \tilde{y})$ is a saddle point under the inequlity in (4.13).

Remark 3. Similar analysis supports the stability behavior of $\hat{E_{2}}(0, \hat{A}, \hat{y})$ on replacing $\tilde{A}$ by $\hat{A}$ and $\tilde{y}$ by $\hat{y}$.

● The Jacobian matrix at the unique positive equilibrium point $E^*(x^*, A^*, y^*)$ is given by

The characteristic equation of the above matrix is,

where

By the Routh-Hurwitz criterion ^[42,43], all the characteristic roots of $J|_{E^*}$ have negative real part iff

Thus, the following theorem can be stated.

Theorem 5. The positive equilibrium point $E^*(x^*, A^*, y^*)$ is locally asymptotically stable, if the conditions in (4.15) hold.

It is shown that conditions in (4.15) hold true if

Remark 4. The model (2.3) is stable around its interior equilibrium $E^*$ if the conditions in (4.16) hold.

In the next theorem, we provide a criterion for global asymptotic stability of the unique positive equilibrium $E^*(x^*, A^*, y^*)$ of the system (2.3).

Theorem 6. The positive equilibrium $E^*(x^*, A^*, y^*)$ of the system (2.3) is globally asymptotically stable below the forthcoming inequlities:

Proof. We select an appropriate positive definite function about $E^*$ as

where $L_{1}$ and $L_{2}$ are positive constants and will be defined later. Now, we differentiate $V$ in context of $t$ along the solutions of system (2.3), we obtain

Choosing $L_{2} = \dfrac{1+ax^*}{c_{1}(1+by^*)}$ and $L_{1} = \dfrac{L_{2}c_{2}\phi}{\phi A^*-\lambda(1-A_{0})}$, we get

Using the Sylvester criterion ^[43], if the inequalities in (4.17) hold, $\dfrac{dV}{dt}$ is negative definite. Hence, $E^*$ is globally asymptotically stable above the inequalities in (4.17) ^[44]. □

4.3.   Hopf-bifurcation and its properties

In dynamical system, Hopf-bifurcation is a occurrence in which stability of an equilibrium point bifurcates into instability and periodic solutions on varing a parameter. The parameter $A_{0}$ is crucial in the system (2.3), we perform a Hopf-bifurcation analysis with $A_{0}$ as the bifurcation parameter and obtain a threshold value $A_{0} = A_{0}^*$ where bifurcation occurs. The following are the necessary and sufficient criteria for Hopf-bifurcation to occur at $A_{0} = A_{0}^*$.

(a) $A_{1}|_{A_{0}^*} > 0, \; \; \; \; A_{3}|_{A_{0}^*} > 0$,

(b) $f(A_{0}^*)\equiv (A_{1}A_{2}-A_{3})|_{A_{0}^*} = 0$,

(c) $Re\Big[\frac{d\eta_{j}}{dA_{0}}\Big]_{A_{0} = A_{0}^*}\neq 0$, where $\eta_{j}, j = 1, 2, 3$ are the roots of Eq (4.14).

We obtain an equation in $A_{0}$ from the expression $A_{1}A_{2}-A_{3} = 0$ and assume it has at least one positive root $A_{0}^*$. Second, for some $\epsilon > 0$, there is an interval containing $A_{0}^*$, $(A_{0}^*- \epsilon, A_{0}^*+ \epsilon)$ such that $A_{0}^*-\epsilon > 0$ and $A_{2} > 0$ for $A_{0} \in (A_{0}^*- \epsilon, A_{0}^*+ \epsilon)$. Therefore, Eq (4.14) is unable to have any real positive root for $A_{0} \in (A_{0}^*- \epsilon, A_{0}^*+ \epsilon)$.

Therefore, when $A_{0} = A_{0}^*$, we can write the Eq (4.14) as follows:

the above equation has the following three roots

where, $\mu = \sqrt{A_{2}}$ and $\theta = -A_{1}$.

For $A_{0}\in(A_{0}^*- \epsilon, A_{0}^*+ \epsilon)$, above roots can be written as

The transversality criterion need to be verified. On differentiation of Eq (4.14) in context of the bifurcation parameter $A_{0}$, we get

where $R = A_{1}A_{2}-A_{3}$ and $\dot A_{s}, s = 1, 2, 3$ denote the derivative of $A_{s}$ in context of time. Thus

Hence $Re\Bigg[\dfrac{d\eta}{dA_{0}}\Bigg]_{A_{0} = A_{0}^*} \neq 0$ if $\Bigg[\dfrac{dR}{dA_{0}}\Bigg]_{A_{0} = A_{0}^*}\neq 0$. This condition implies that the rate of change in the real part of characteristic roots at the bifurcation point is non-zero. The forthcoming theorem can be expressed as a result.

Theorem 7. The system experiences a Hopf-bifurcation near the positive equilibrium $E^*(x^*, A^*, y^*)$, if the necessary and sufficient conditions $(a)$–$(c)$ are satisfied. The threshold value of the bifurcation parameter $A_{0}$ is provided by the equation $f(A_{0}^*) = 0$.

From Theorem 7, the system (2.3) experiences a Hopf-bifurcation at $A_{0} = A_{0}^*$. Now we will discuss the direction and stability of bifurcating periodic orbit ^[45]. The Eigenvectors $v_{1}$ and $v_{3}$ of the variational matrix $J|_{E^*}$ corresponding to the characteristic roots $\eta_{1, 2} = \pm \iota \mu$, and $\eta_{3} = -A_{1}$

where

Now, we consider the transformation in system (2.3) in terms of $X$, $Y$ and $Z$ with $x = x^*+r_{11}X+r_{12}Y+r_{13}Z$, $A = A^*+r_{21}X+r_{22}Y+r_{23}Z$ and $y = y^*+r_{31}X+r_{32}Y+r_{33}Z$, the system (2.3) simplifies

where $\Lambda = r_{22}r_{33}-r_{32}r_{23}+r_{21}r_{32}-r_{31}r_{22}$,

The system (4.18) has an equilibrium point $(0, 0, 0)$, and the variational matrix calculated at this point is as follows

where $\theta = \dfrac{\partial \Theta_{3}}{\partial y}$. At $A_{0} = A_{0}^*$ and $(x, A, y) = (0, 0, 0)$, The values of $g_{11}$, $g_{02}$, $g_{20}$, $G_{21}$, $e_{11}$, $H_{101}$, $H_{110}$, $e_{20}$, $\mu$, $\kappa_{20}$, $\kappa_{11}$ and $g_{21}$ can be calculated by following relations

We can calculate the value for $\kappa_{11}$ and $\kappa_{20}$ by solving $(\theta-2\iota \mu)\kappa_{20} = -e_{20}$ and $\theta \kappa_{11} = -e_{11}$. Based on the above analysis, we can compute the following quantities

From the above analysis of Hopf-bifurcation, we now state the following theorem.

Theorem 8. (1) The Hopf-bifurcation curve exists for $A_{0} > A_{0}^*$ ($A_{0} < A_{0}^*$) if $\mu_{2} > 0$ ($\mu_{2} < 0$).

(2) The bifurcating perodic solutions are stable (unstable) if $\beta_{2} < 0$ ($\beta_{2} > 0$).

(3) The bifurcating perodic solution increases (decreases) if $T_{2} > 0$ ($T_{2} < 0$).

5.   Computer simulation

In this section, we will conduct computer simulations using MATLAB R2022a to confirm the validity of our theoretical finidings discussed in the preceding sections. For the simulation, we have selected the forthcoming parameter values.

For the above set of parameters, the inequalities in Theorem 3 hold true, therefore the positive equilibrium point $E^* (0.2863, 1.5289, 3.771)$ exist uniquely. Since the conditions in Theorem 5 are also followed by $E^*$ $(A_{1} = 3.0494 > 0, \; A_{3} = 0.0361 > 0, \; A_{1}A_{2}-A_{3} = 0.3917 > 0)$, hence it is locally asymptotically stable. The stability nature of $E^*$ is depicted in Figure 3. The left panel illustrates the time series evolution of population and the right panel illustrates phase portrait in the $xAy$-space. The figure illustrates that both species fluctuate around their equilibrium level for a finite time before settling. It illustrates the local asymptotic stability of the system (2.3) around the interior equilibrium.

Now, to investigate the effect of the degree of dependency $A_{0}$, we increase it keeping all other parameters unchanged. We noticed that the transient period is increasing with $A_{0}$, which indicates a route to Hopf-bifurcation. At $A_{0} = 0.2872(A_{0}^*)$ the system undergoes a Hopf-bifurcation and populations face perpetual oscillations. The necessary and sufficient conditions given in Theorem 5 are also satisfied. The bifurcation value of $A_{0}$ can be calculated by solving the equation $f(A_{0}) = 0$. For further values of $A_{0}$, the system has a limit cycle. In Figure 4, the unstable nature of $E^*$ and limit cycle are plotted. A solution trajectory, which is initiated from near to the positive equilibrium is attracted by the limit cycle around it (see Figure 4c). This shows that the Hopf-bifurcation is supercritical. If we further increase the parameter $A_{0}$ then at the higher dependency of predators on prey the system again get back to its stability near $E^*$ via another Hopf-bifurcation. Beyond another critical value of $A_{0}$, $A_{0}^{**} = 0.9291$, the limit cycles leave the system and $E^*$ becomes stable. The stable behavior of $E^*$ at $A_{0} = 0.955 > A_{0}^{**}$ is illustrated in the Figure 5.

As we vary parameter $A_{0}$, we can see that the system changes its stability twice. Thus, the Hopf-bifurcation occurs multiple times in the full range of the dependency factor. In order to provide a comprehensive representation of the Hopf-bifurcation, we have constructed a three-dimensional bifurcation diagram as depicted in Figure 6. This diagram illustrates the attractors of the system in the prey and predator plane for many values of parameter $A_{0}$. It is remarkable that the positive equilibrium point remains stable below the first critical value of $A_{0}$. However, once the critical value is exceeded, the trajectories are attracted toward a limit cycles. These stable limit cycles are represented by the colored closed orbits displayed in the diagram. When $A_{0} = 0.2872 = A_{0}^*$, we can see that all the conditions of Theorem 7 are satisfied. Also we comutate that

By Theorem 8, this result implies that the Hopf-bifurcation taking place, is of supercritical nature, and bifurcating periodic solutions are stable. When $A_{0} = 0.9291 = A_{0}^{**}$, Hopf-bifucation occurs and the limit cycle disappears. We can see that the all conditions of Theorem 7 are satisfied and

By Theorem 8, this result implies that the Hopf-bifurcation taking place is of supercritical nature, and bifurcating periodic solutions are stable.

From Figure 6, we observe that at the low degree of dependency, the predator has a higher dependency on additional food, resulting in the system (2.3) being stable. However, as the dependency factor gradually increases, the system (2.3) is unstable and forms a limit cycle. This indicates that the predator's dependence extends to both the prey and additional food resources. Finally, as we further increase the value of $A_{0}$, the predator's reliance shifts more towards the prey, leading to the system (2.3) returning to a stable state. Overall, the bifurcation diagram visually demonstrates how the system (2.3) dynamics change as the dependency of predator on prey varies, showing the changes from stable to unstable states as well as the effect of predator dependency factor on both prey and supplementary nourishment.

The fear parameter $k$ is also a vital parameter that significantly influences the dynamics of the system. For the best illustration of varing $k$, we change the parameter $A_{0} = 0.4$, whereas other parameters are the same as in (5.1). Initially, for lower values of $k$, the system is stable near its positive equilibrium. Now, we start increasing the fear parameter, and this change in $k$ opposes the stability and at $k = 2.5441(k^*)$, the system falls into the never ending oscillations via a Hopf-bifurcation. If we choose $k$ beyond $k^*$, then $E^*$ becomes unstable and the trajectory started from near by points approaches a stable limit cycle. The unstable equilibrium point and stable limit cycle at $k = 6 > k^*$ are shown in the Figure 7.

Figure 8 provides a comprehensive bifurcation diagram that illustrates the behavior changes as the fear parameter $k$ rises. It is evident from the figure that the system remains stable for lower values of $k$. However, as $k$ approaches its threshold value ($k^*$), the system destabilizes and experiences a Hopf-bifurcation. Beyond the critical value, the equilibrium point becomes unstable and a periodic solution appears in the system. When $k = 2.552 = k^*$, we can see that all the conditions of Theorem 7 are satisfied. Also, we computate that

By Theorem 8, the bifurcating periodic solution exists for $k < k^*$, the bifurcating periodic solution increases and bifurcating periodic solutions are stable.

The behavior of system (2.3) experiences Hopf-bifurcation in the context of the parameters $A_{0}$ and $k$. Hence, $A_{0}k$-plane can be divided into two types of regions, region of stability where the equilibrium point is stable and the region of instability where the equilibrium point is unstable and limit cycle presents in the system. Both the regions are separated by a curve, where the characteristic roots of variational matrix are purely imaginary, known as Hopf-bifurcation curve. Previously, we noticed that the system switches its stability twice on varying parameter $A_{0}$, while only once on varying the fear parameter $k$. Therefore, there are two stability regions and one instability region in the $A_{0}k$-plane separated by two Hopf-bifurcation curves. The stability and instability regions are represented in Figure 9. When we consider a particular value of the cost of fear $k = 3$, there are two bifurcation points for the parameter $A_{0}$: $A_{0}^{'} = 0.3735$ and $A_{0}^{''} = 0.9418$. System (2.3) remains stable until $A_{0}^{'}$ as depicted in the Figure 9. However, when $A_{0}$ exceeds $A_{0}^{'}$ ($A_{0} > A_{0}^{'}$), the system becomes unstable and enters a limit cycle. If we further increase $A_{0}$ beyond $A_{0}^{''}$ ($A_{0} > A_{0}^{''}$), the system (2.3) regains its stability.

Now, to maintain the equilibrium levels of prey ($x^* = 2$) and predator ($y^* = 2$), we vary the parameters $k$ and $A_{0}$ whereas the other parameters are fixed (same as in (5.1)). We notice that we have to decrease the value of parameter $A_{0}$ as the parameter $k$ increases to maintain the population density of prey or predator (refer the Figure 10).

6.   Conclusions

The prey-predator system in ecology describes the relationship between species, where the prey population supports the predator population as a food source, and the predator population consumes the prey population through predation. When additional food is introduced into a prey-predator system, it can lead to growth in the prey population, which subsequently supports a larger predator population due to more available food resources. This can result in a more stable or balanced system compared to a system without additional food. Numerous studies ^{[20,30,31,32]} have investigated the effects of additional food in prey-predator systems, providing valuable insights into how this factor can influence population dynamics and ecological interactions between the species involved. Recent studies ^[37,38,39] have examined the effect of predator-induced fear on prey individuals, leading to a reduction in their reproducion rate. The separation of the fear effect from direct killing by predators was made feasible due to some extensive experimental efforts by various ecologists ^[33]. The Crowley-Martin functional response model investigates how predator predation rates change with variations in prey density and shows more realistic dynamics ^[13,26,27]. In this research, we examine a prey-predator model that incorporates three important factors: Fear effect in the prey species induced by predators, predator-prey interaction by the Crowley-Martin functional response, and the influence of providing additional food to the predator. The presence of additional food may influence the predator's dependency on prey, which may lead to a significant difference in the dynamics of the system. It is essential to understand how much the predator depends on prey and additional food. We introduced a new parameter, the degree of dependency, which indicates the predator's dependency on prey. After conducting a thorough analysis, we have demonstrated that the proposed system is characterized by positivity, boundedness and persistence, highlighting its ecological significance and well-posedness. We conducted an analysis to examine how the number of interior steady states and their stability characteristics change with variations in parameters. Our system exhibits four types of equilibria: The trivial equilibrium $E_{0}(0, 0, 0)$, the axial equilibrium $E_{1}(x^*, 0, 0)$, two prey-free equilibria $\tilde {E_{2}}$ and $\hat {E_{2}}$ under condition (4.3) and a unique positive equilibrium $E^*$ under conditions (4.8)–(4.10). We have carried out the local and global stability of all existing equilibria in the system. With the stability results, we conclude that the coexistence of both species depends on the equilibrium level of provided additional food and the degree of dependency (refer to the set of inequalities in (4.16) and (4.17)). This analysis carries various conditions with the stability of the equilibrium being influenced by the degree of dependency.

We have shown that the degree of dependency plays a vital role in stabilize the system through Hopf-bifurcation. The Hopf-bifurcation phenomenon has also been validated by a computer simulation. We noticed, for a particular range of cost of fear, that the system switches its stability around the positive equilibrium point multiple times (see Figure 9). This indicates that the system's stability relies on the absence or minimal degree of dependency between the two species. If the dependency exceeds a certain initial threshold, both species begin oscillating around unstable steady states. However, if the dependency attains another threshold, the system regains stability. This provides a clear demonstration of how changes in the degree of dependency impact the ecological dynamics in the predator-prey system.

Additionally, we explored the Hopf-bifurcation with respect to the fear parameter $k$ and obtained the threshold value $k^* = 2.58$, where the system destabilizes. This shows that the system is stable if the degree of fear is low, but if the fear is increased beyond a level, then both the species start oscillating around their unstable steady state. From the ecological point of view, this nature can be justified. Initially, when $k$ is very small, both populations coexist with stability. Now, we start increasing the cost of fear, which leads to a reduction in the reproduction rate of prey. Due to this, the prey population starts decreasing and consequently, the predators also start decreasing due to the scarcity of food. At the low density of predators, the fear effect gets weak and the prey population starts growing. Now, the abundance of food predators also increases. After a certain level of cost of fear, both the populations oscillate perpetually. This leads to limit cycles in the system at the higher values of cost of fear.

In conclusion, we find that predator-induced fear and the predator's dependency on prey both are effective in controlling the oscillations in the system, thereby playing a crucial role in stabilizing or balancing the ecosystem. These factors contribute significantly to the regulation and persistence of predator-prey dynamics, highlighting their importance in ecological stability. The outcome of our study has an important effect on the conservation of prey populations, as they are supported by a system with additional food for predators. We believe that our research will deepen the understanding of how ecological dynamics change in the presence of fear effect, additional food resources, and degree of dependency are added.

Use of AI tools declaration

The authors declare that they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

Author's are thankful to anonymous reviewers for careful reading and constructive suggestions that improved the quality and presentation of the paper.

Conflict of interest

The authors declare that they have no competing interests regarding the present study.