Combined impact of fear and Allee effect in predator-prey interaction models on their growth

1.   Introduction and motivation

Ecologists have observed the prey-predator interaction cycle with great interest. In past, they believed that the prey population declined due to direct killing by predators. In recent years, they realized through field studies that the mere presence of predators affects the growth of the prey population. The effect of stress in the reduction of the prey population is more powerful than direct killing. It has been pointed out that environmental stress affects the physiology and behavior of the prey to such an extent that it causes delayed mammal reproduction. Fear of the presence of predators, social life, competition for food, and similar factors are all causes that can lead to a reduction in mammal reproduction ^[1,2]. For small size prey such as sparrows and bank voles, the anti-predation response influence on the prey population is much higher ^[2,3,4].

Prey species where the physical prey size is small make coordinated feeding groups and avoid being detected by predators. They are vulnerable to a greater range of predators and cannot defend themselves or outrun predators. On the other hand, prey species where the physical prey size is large tend to form a large feeding group. They can defend themselves against predators and are more likely to depend upon group defense, group alertness, self-defense within a group, and speed to avoid being killed by a predator. All prey of either small or large physical size respond to predation risk and take a variety of anti-predation measures. They change habitat, practice vigilance whilst foraging, change their behavior whilst mating, and also undergo physiological changes ^[5,6]. It has been observed biologically, based on experimental studies, that mammals suppress breeding in response to strong predation pressure ^[4,7,8,9]. This has a long-term effect on the population of prey. The nonbreeding individuals have a better chance of avoiding predation than those who are breeding. Predation induced breeding suppression affects the prey-predator dynamics. Individuals who do not breed experience less predation because during pregnancy individuals are less active and can be captured more easily.

There have been quite a few previous attempts to incorporate the effect of fear of the predator into predator-prey systems. Kumar and Dubey ^[10] studied a mathematical model to investigate the fear effect and prey refuge in a predator-prey system with gestation time delay. Their analysis demonstrated that the fear effect in the prey population induces Hopf-bifurcation. The combined effects of fear and group defense in a fractional order prey-predator system are investigated by Das and Samanta ^[11]. Wang et al. ^[12] examined a predator-prey model by incorporating the cost of fear of the predator on prey. Their analysis found that a rise in fear stabilizes the system and a drop in fear leads to bistability. Wang and Zou ^[13] incorporated maturity after a delay and adaptive defense in the predator-prey fear model. They found that a higher predator population leads to robust anti-predator defense and higher predation risk implies a weak anti-predatory defense. Mondal et al. ^[14] studied a predator-prey system with a fear of prey and alternative food available for predators. According to their findings, alternative food to predators lowers the attack on prey and thus helps in the growth of both populations. Zhang et al. ^[15] and Wang et al. ^[16] incorporated prey refuge with fear effect. Duan et al. ^[17] investigated that maturation delay in a fear model leads to Hopf-bifurcation. The fear amongst the prey can stabilize as well as destabilize the population model as investigated by Pal et al. ^[18]. Pal et al. ^[19] observed that the fear effect induced by cooperative hunting destabilizes the system. Das and Samanta ^[20] examined a stochastic model where fear is in a prey population and an alternative food source is available for predators.

A lot of papers build on the previous work of Wang et al. ^[12] to introduce fear of the predator into the predator-prey system. If $N$ represents the total population size and $P$ the total predator population size, their fear function reduced the growth rate of the prey population by a multiplicative factor

so in their model the growth rate of the prey population is given by

where $r$ represents the per capita birth rate of the prey population in the absence of the predator. This fear function was used by many authors ^{[21,22,23,24,25]}. However, it has the drawback that as $P$ becomes very large (in other words, the predator population gets very large), the growth rate of the prey population tends to zero. This is not ecologically true since the prey will become habituated to coexist with the predator, and as the size of the predator population increases the growth rate will reduce, but not to zero.

Sarkar and Khajanchi ^[26] corrected this by introducing a fear function

where $\alpha$ represents level of fear and $\eta \in [0, 1]$ represents the minimum cost of fear. They incorrectly state that as the level of fear, $\alpha$, increases the fear function decreases. Although ecologically this should be true, with their fear function it is mathematically not true, as the fear function increases as the level of fear increases. We make a minor redefinition to the model of Sarkar and Khajanchi and define our fear function to be

where $\alpha_1$ is a parameter and $s_{\alpha_1}$ is the level of fear. This fear function decreases as the level of fear increases. The fear function which we are using is more ecologically sound than the fear function of of Wang et al. ^[12], because after a certain level of fear, the prey becomes habituated to live with the predators. Thus, it reduces its reproduction to a certain extent but does not stop it completely.

The paper is structured as follows. In Section 2 we first outline the full eco-epidemiological model that we will be studying. The model combines an Allee effect and the modified fear function. There is a disease spreading amongst the prey, consequently there are three population classes, susceptible prey, infected prey, and predators. We also introduce the full eco-epidemiological model with a time delay corresponding to the delayed effect of predators consuming prey on predator reproduction. In Section 3 we will focus on and analyze the model with no disease present so that it is a predator-prey model with no disease amongst the prey. There is no time delay in the model discussed in this section. Section 4 analyzes the model with no disease present and a time delay. Section 5 will return to the full eco-epidemiological model with disease amongst the prey but without a time delay, and Section 6 analyzes this model with a time delay. This includes finding the length of the time delay to preserve the stability. In all of these models the reproduction of the prey is reduced due to our modified fear function as well as Allee effects. A final conclusion section summarizes and discusses the paper.

2.   Formulation of the eco-epidemiological model

In this section we will outline the full eco-epidemiological model that we will be studying. All prey are born susceptible. The reproduction rate of the prey depends on the combined effects of fear of the predator together with the strong Allee effect. The function $g(\alpha_1, \eta, y) = \eta+\frac{\alpha_1 (1-\eta)}{\alpha_1 + y}$ represents the cost of fear of predators on the reproduction of prey ^[26]. $\eta \in [0, 1]$ represents the minimum cost of fear and $\alpha_1$ is a parameter for which $\frac{1}{\alpha_1} = s_{\alpha_1}$ is the fear factor. In the case of a high predator population, the prey survives under the minimum fear $\eta$, and after a certain level of fear, the prey becomes habituated to living with the predators.

In the full eco-epidemiological model, we will study a predator-prey model where the prey is infected with some disease. We divide the total prey population into two classes called susceptible prey and infected prey. Let $u_1(t_1), u_2(t_1)$, and $v(t_1)$ be the uninfected prey density, infected prey density, and predator density, respectively at time $t_1$. We will formulate the model using the following assumptions:

(a) The uninfected prey population grows logistically when there is no predator.

(b) The susceptible prey population $u_1(t_1)$ becomes infected when it comes in contact with infected prey $u_2(t_1)$. This infection follows the nonlinear incidence rate. If we assume that the incidence rate is given by the law of mass action $\beta_1 u_1f(u_2)$, then $f(u_2)$, will increase as $u_2$ increases, which is not realistic. Lui et al. ^[27] investigated many reasons for taking a nonlinear incidence rate rather than a bilinear one.

(c) The infected prey population will remain infected and will not recover.

(d) Infected prey will not be able to reproduce because it has the disease. Only uninfected prey species have the capability to reproduce.

(e) The predator consumes infected prey more than susceptible prey because infected prey is easily catchable. For example, wolves attack moose more successfully when they are heavily infected by Echinococcus granulosus^[28]. Hence, we assumed that the consumption rate $p_1$ of infected prey is more than the consumption rate $b_1$ of uninfected prey.

(f) The conversion rates $b_1$ and $p_1$ represent the reproduction of predators for capturing each uninfected and infected prey species, respectively.

(g) All the parameters used in the model are positive and $u_1(0) > 0, u_2(0) > 0$, and $v(0) > 0$.

We will investigate the following ecological issues through our model:

(1) The effect of predator fear on the stability of the system.

(2) The level of infected prey.

Using the above assumptions, we will formulate a mathematical model as below:

where $u_1$ represents the susceptible prey population, $u_2$ represents the infected prey population, $v$ represents the infected predator population, $k$ is the carrying capacity in the absence of infected prey, $\theta_1$ is the Allee threshold, $\eta$ represents the minimum cost of fear, $\alpha_1$ is a parameter for which $\frac{1}{\alpha_1} = s_{\alpha_1}$ is the fear factor, $\beta_1$ represents the the disease transmission coefficient, $c_1$ is the half saturation constant, $\delta$ is the predator preference rate of $u_2$, and $a_1$ represents the predation rate between susceptible prey and predator. $g_1$ represents the predation rate between infected prey and predator, $\mu_1$ is the total per capita death rate of infected prey which includes both natural death and disease induced additional mortality, $b_1$ is the consumption rate of uninfected prey, $p_1$ represents the consumption rate of infected prey, and $m_1$ represents the predator per capita natural mortality rate.

We additionally assume that

(h) $u_1(0)+u_2(0) < k$.

Note that if there are no infected prey, then $u_1(t)$, the total size of the susceptible prey population, will always be less than the carrying capacity. So, as normally the disease will be introduced by only a few infected prey, it makes sense to assume that $u_1(0) + u_2(0) < k$, which ensures that $u_1(t) + u_2(t) < k$ for all time.

We consider in our study the following fear function to measure the cost of fear:

where $\frac{1}{\alpha} = s_\alpha$ is the nondimensional form of the fear factor. Figure 1 demonstrates that the fear function decreases as the level of fear $s_\alpha$ increases. Also, as the predator population increases, the fear function declines. Because the fear factor multiplies with prey growth rate in model ($3.1$), the growth of the prey population will be low when the value of the fear function is low, that is, when the predator population is high.

To make the model simple, we non-dimensionalize model ($2.1$) by taking $x = \frac{u_1}{k}$, $y = \frac{u_2}{k}$, $z = \frac{v}{k}$, $t = t_1kr$, $\theta = \frac{\theta_1}{k}$, $\beta = \frac{\beta_1}{kr}$, $c = \frac{c_1}{k}$, $a = \frac{a_1}{r}$, $g = \frac{g_1}{r}$, $\mu = \frac{\mu_1}{kr}$, $b = \frac{b_1}{r}$, $p = \frac{p_1}{r}$, $m = \frac{m_1}{kr}$, and $\alpha = \frac{\alpha_1}{k}$ as the parameter for which $\frac{1}{\alpha} = s_{\alpha}$ is the fear factor. So, model ($2.1$) in nondimensional form becomes:

In Eq (2.2), $axz$ is the total rate at which susceptible prey are consumed by predators, and $bxz$ is the total rate at which predators increase due to the consumption of susceptible prey. So $a$ is the per capita rate at which a single predator attacks and consumes each susceptible prey, and $b$ is the additional per capita rate at which each predator reproduces due to this consumption. As predators are typically much larger than prey, one predator consuming one prey is unlikely to give the predator sufficient energy to reproduce an entire predator. Therefore, $b \le a$, and similarly $p \le g$. Another way to see this is that the conversion rate of prey into predators will always be less than the rate at which predators attack prey because some of the intake energy gained by eating prey will go to waste since it will not all be converted into predator growth. Similar assumptions were made in ^[29] and ^[30].

We also introduce a time delay to make the model (2.1) more realistic. There is a time delay $\tau$ between predators consuming susceptible or infected prey and the susceptible or infected prey which has been consumed being converted into predators through reproduction. The population of predators will grow with interactions with uninfected prey and infected prey which take place at time $t-\tau$. Hence, the population dynamics for the interaction of uninfected prey, infected prey, and predators can be described by the following system of differential equations:

All the variables and parameters have the same meaning as given in ($2.1$). After non-dimensionalizing model ($2.3$) as we did in model ($2.1$), we get:

We now proceed to analyze this model in four cases: the model with no disease present and no time delay, then the model with no disease present and a time delay, next the full eco-epidemiological model but with no time delay, and finally the full eco-epidemiological model with a time delay.

3.   The model with no disease present and no time delay

In this section we focus on the model with no disease present and no time delay. Putting $y = 0$ in Eq (2.2), we obtain

For the model with no disease present, we make the further transformations to simplify the model $N = x$, $P = az$, $\alpha' = \frac{\alpha}{a}$, $T = t \frac{b}{a}$, $r_0 = \frac{b}{a}$, and $d = \frac{ma}{b}$. Then, these equations become

these are the equations which we shall work with in the remainder of this section. For notational simplicity, we write $\alpha$ for $\alpha'$ and $t$ for $T$.

3.1.   Boundedness

Theorem 1. All the solutions $(N(t), P(t))$ of the system ($3.1$) with initial conditions $N(0)\geq 0$, $P(0)\geq 0$ are bounded for all $t > 0$.

Proof. We define $w = N+\frac{1}{r_0} P$. The derivative of $w$ along with $(3.1)$ is

For any positive constant $s$, we have

Here, $N < 1$, and if $s = min\left(\frac{\theta}{r_0}, d \right)$, then

Let $\frac{\theta+1}{r_0} = L$. Then, the above differential inequality can be expressed in the form

Using the theory of differential inequality for $w(t)$, we obtain

where $w_0 = w(N_0, P_0)$. We can deduce that for $t\rightarrow \infty$, $0 < w\leq \frac{L}{s}.$

For any $\epsilon > 0$, define

This shows that the solution of the system represented by ($3.1$) is bounded, i.e., for every trajectory of ($3.1$), there is a time $t_0$ such that $(N(t), P(t))\in \Delta$ for all $t > t_0$. Hence, the theorem is proved. □

3.2.   Equilibrium points of the system

It can be verified that model ($3.1$) has four possible equilibria:

$(1)$ The extinction equilibrium $E_0 = (0, 0)$.

$(2)$ The axial equilibrium $E_\theta = (\theta, 0)$.

$(3)$ The axial equilibrium $E_1 = (1, 0)$.

$(4)$ The coexistence equilibrium $E^* = (N^*, P^*)$, where $N^* = d$ and $P^*$ is a nonnegative root of the quadratic equation: $P^{*2}+ \sigma _1 P^*+\sigma _2 = 0$. Define $\sigma _1 = \alpha - \eta(1-d)(d-\theta)$ and $\sigma _2 = -\alpha (1-d)(d-\theta)$. By Descartes' rule of signs, the quadratic equation has a unique nonnegative root if $\theta < d < 1$, as shown in Figure 2.

3.3.   Stability analysis

The stability matrix ($J$) for system ($3.1$) is

where $S_1 = \frac{1}{r_0}\left(\eta+\frac{ 1-\eta}{1 +s_\alpha P}\right)[(1-2N)(N-\theta)+N(1-N)]-\frac{P}{r_0}$ and $S_2 = -\frac{1}{r_0}\left[N(1-N)(N-\theta)\frac{\alpha (1-\eta)}{(\alpha + P)^2}+N\right]$.

$(1) \; E_0$ is locally asymptotically stable because the characteristic roots are –$\frac{\theta}{r_0} < 0$ and $-d < 0$.

$(2)$ The characteristic roots for $E_\theta$ are $\frac{\theta (1-\theta)}{r_0}$ and $\theta-d$, $\frac{\theta (1-\theta)}{r_0} > 0$ because $0 < \theta < 1$ from system ($3.1$). Thus, $E_\theta$ is unstable.

$(3) \; E_1$ is stable if $d > 1$ because the characteristic roots are –$\frac{(1-\theta)}{r_0} < 0$ and $1-d$.

$(4)$ The stability matrix for the coexistence equilibrium $E^* = (N^*, P^*) = (d, P^*)$ is:

The corresponding characteristic equation is $\lambda ^2+A\lambda+B = 0$, where $A = -\frac{1}{r_0}\left(\eta+\frac{ 1-\eta}{1 +s_\alpha P^*}\right)[-2d^2+(\theta +1)d]$ and $B = \frac{P^*}{r_0}\bigr[d(1-d)(d-\theta)\bigr(\frac{\alpha (1-\eta)}{(\alpha + P^*)^2}\bigl)+d\bigl]$. The Routh-Hurwitz criteria for the second order system is given by: $A > 0$ and $B > 0$.

$\quad A = -\frac{1}{r_0}\left(\eta+\frac{ 1-\eta}{1 +s_\alpha P^*}\right)[-2d^2+(\theta +1)d] > 0$ if $d > \frac{\theta +1}{2}$ and because $(N^*, P^*)$ exists for $\theta < d < 1$, $B = \frac{P^*}{r_0}\bigr[d(1-d)(d-\theta)\bigr(\frac{\alpha (1-\eta)}{(\alpha + P^*)^2}\bigl)+d\bigl] > 0$. Therefore, $(N^*, P^*)$ is locally asymptotically stable if $d > \frac{\theta +1}{2}$, i.e., twice the death rate of predator is one unit higher than the Allee threshold value as shown in Figures 3 and 4.

3.4.   Bifurcation analysis

We observed a transcritical bifurcation and a Hopf-bifurcation that occurs around the equilibrium points when the parameter passes through some critical values.

3.4.1.   Transcritical bifurcation

According to Sotomayor's theorem ^[31], we investigate the transcritical bifurcation for the system ($3.1$) taking $d = d_0$ as the bifurcation parameter. If the following conditions are satisfied:

$1) \; w^T f_d(E_1, d_0) = 0$,

$2) \; w^T (Df_d(E_1, d_0). \upsilon)\neq 0$,

$3) \; w^TD^2f(E_1, d_0)(\upsilon, \upsilon)\neq 0,$

then the system ($3.1$) experiences a transcritical bifurcation at the equilibrium point $E_1 = (1, 0)$ as the parameter $d$ passes through the bifurcation value $d = d_0 = 1$, where $f_d$ denotes the vector of partial derivatives of the components of $f$ with respect to the scalar $d$ and $Df$ represents the matrix of partial derivatives of the components of $f$ with respect to the components of $X = (N, P)$. Here,

$Df(X) = \begin{pmatrix} \frac{\partial f_1}{\partial N} & \frac{\partial f_1}{\partial P}\\ \frac{\partial f_2}{\partial N} & \frac{\partial f_2}{\partial P}\\ \end{pmatrix}$, ${D^2f(X) = \begin{pmatrix} \frac{\partial^2 f_1}{\partial N\partial N } & \frac{\partial^2 f_1}{\partial N\partial P } & \frac{\partial^2 f_1}{\partial P\partial N } & \frac{\partial^2 f_1}{\partial P\partial P }\\ \frac{\partial^2 f_2}{\partial N\partial N } & \frac{\partial^2 f_2}{\partial N\partial P } & \frac{\partial^2 f_2}{\partial P\partial N } & \frac{\partial^2 f_2}{\partial P\partial P } \end{pmatrix}}.$

$\upsilon$ is an eigenvector of $A = Df(\overline{E}_1, d_0)$ corresponding to the eigenvalue $\lambda_2 = 0$ and $w$ is an eigenvector of $A^T$ corresponding to the eigenvalue $\lambda_2 = 0$.

Theorem 2. The system $(3.1)$ undergoes a transcritical bifurcation at the positive equilibrium $E_1 = (1, 0)$, as the parameter $d$ passes through the bifurcation value $d = d_0 = 1$.

Proof. The linearized system around the equilibrium point $E_1$ is:

Let us define $\upsilon = (\upsilon_1, \upsilon_2)^T$ and $w = (w_1, w_2)^T$ to be the right and left eigenvectors of $\lambda_2 = 0$. Now solving $J(E_1, d_0)\upsilon = 0$ implies that $\upsilon = \bigr(\frac{-\upsilon_2}{1-\theta}, \upsilon_2\bigl)$, and by solving $J^T(E_1, d_0)w = 0$, we get $w = (0, w_2)$ where $\upsilon_2, w_2$ are any nonzero real numbers. Now, system $(3.1)$ can be rewritten as in the following vector form:

Taking the derivative of $f(X)$ with respect to $d$, we get $f_d(X) = \begin{pmatrix} 0\\ -P \end{pmatrix}$, then $f_{(E_1, d_0)}(X) = \begin{pmatrix} 0\\ 0 \end{pmatrix}$. Hence $w^T f_{(E_1, d_0)}(X) = 0$. Next, taking the derivative of $f_d(X)$ with respect to $X = (N, P)^T$, we get

We have $w^T (Df_{(E_1, d_0)}(X). \upsilon) = -w_2\upsilon_2\neq 0$. Furthermore,

Hence,

where $(\upsilon, \upsilon)$ is a tensor product of $(\upsilon_1, \upsilon_2)^T$, so $w^TD^2f_{E_1, d_0}(X)(\upsilon, \upsilon) = (\upsilon_1\upsilon_2+\upsilon_2\upsilon_1)w_2 = 2\upsilon_1\upsilon_2 w_2\neq 0$. Therefore, according to Sotomayor's theorem ^[31] for local bifurcation, system $(3.1)$ has a transcritical bifurcation at steady-state $E_1$ when the parameter $d$ passes through the bifurcation value $d_0 = 1$.

□

3.4.2.   Hopf-bifurcation

Theorem 3. Any system of the form

that has an equilibrium $(N^*, P^*)$ for the parameter $d$ whose linearization has eigenvalues $\lambda_{1, 2} = \gamma(d) \pm i\beta(d)$ such that $\gamma(\overline{d}) = 0$, $\beta(\overline{d}) = \overline{\beta} > 0$ and satisfying the following condition:

experiences Andronov-Hopf bifurcation.

Proof. From the previous analysis, we know that the two equilibrium points $(0, 0)$ and $(1, 0)$ have no complex eigenvalues. Therefore, we use the coexistence fixed point $(N^*, P^*)$ to check if the conditions stated in the above theorem apply to the model.

The possibility of the Hopf-bifurcation at $E^* = (N^*, P^*)$ has been analyzed by taking $d$ as a bifurcation parameter and keeping the rest of the parameters as constants. To begin, we investigate the linearization condition, which is

The eigenvalues of $(3.3)$ in terms of the trace and determinant are given by

Let $\gamma(d) = \frac{trJ}{2}$. For the condition of the theorem to be satisfied, we try to find a value $d = \overline{d}$ such that $\gamma(\overline{d}) = 0$.

Thus, there exists a value $d = \overline{d}$ at which $\gamma(\overline{d}) = 0$, or similarly the real part of the eigenvalues is $0$.

Next, we consider $\frac{d\gamma(d)}{d(d)}\Big\vert_{d = \overline{d}}$. We have $\gamma(d) = \frac{1}{2r_0}\left(\eta+\frac{ 1-\eta}{1 +s_\alpha P^*}\right)[-2d^2+(\theta +1)d]$.

Therefore, the condition $\frac{d\gamma(d)}{d(d)}\Big\vert_{d = \overline{d}}\neq 0$ is satisfied.

Now we show that the complex part of the eigenvalues exist at $d = \overline{d}$ and they are not zero. Let $i\beta (d) = \pm \frac{i\sqrt{4detJ-tr^2J}}{2}$. We know that at $\overline{d}$, $trJ = \gamma (\overline{d}) = 0$, which implies that $tr^2J = 0$. Therefore

Next, we find $detJ(d)$:

Hence, the conditions for the existence of the Hopf-bifurcation are satisfied. □

We run the bifurcation diagram using XPPAUT software and using the parameter $d$ as the bifurcation parameter. Figure 5 shows the Hopf-bifurcation at the coexistence equilibrium $E^*$ when $d = \frac{1+\theta}{2}$ and the transcritical bifurcation at the equilibrium $E_1$ when $d = 1$.

3.5.   Effect of cost of fear $s_\alpha$ and the Allee threshold effect $\theta$

It is clear from Figure 6 that as the level of fear $s_\alpha$ increases, the initial short-term decrease of the prey population becomes larger. However, after a long period of time, the level of fear has no effect on the long-term level of the prey population as the prey become habituated to the presence of the predator. However, as the level of fear increases, the long-term level of the predator population decreases.

The existence of the coexistence equilibrium point $E^* = (N^*, P^*)$ implies that $\theta < d < 1$. If the Allee threshold $\theta$ is small, the coexistence equilibrium is stable when it exists and the Allee effect has no impact on the long-term prey density, as shown in Figure 7. Moreover, for small values of the Allee threshold, the long-term level of the predator population decreases. However, as the Allee threshold increases further, both the long-term predator and long-term prey populations go extinct.

4.   Model with no disease and a time delay

Next, we consider the model with no disease present and a time delay. Putting $y = 0$ in Eq ($2.4$) and transforming as in Eq ($3.1$), we get:

At equilibrium $\overline{N} = d$. Let $N = \overline{N}+u$ and $P = \overline{P}+v$ where $u$ and $v$ are very small perturbations. Then substituting into the second equation of (4.1), we deduce that

Let $v = Be^{\lambda t}$ and $u = Ae^{\lambda t}$. Substituting into the equations and dividing by $e^{\lambda t}$,

The stability matrix is $\begin{pmatrix} L-\lambda & M\\ \overline{P}e^{-\lambda \tau} & -d+\overline{N}de^{-\lambda \tau}-\lambda \end{pmatrix}$, where $\overline{N} = d$. The corresponding characteristic equation is:

Let $\lambda = \sigma+i\rho$. Substituting into the characteristic equation and expanding,

Equating real and imaginary parts,

At $\tau = \tau_1^*$, $\sigma(\tau_1^*) = 0$,

Squaring and adding the two parts of ($4.4$), we obtain

Let $\xi = \rho_1^{*2}$, then ($4.5$) reduces to

If $\rho_1^{*2} = \xi$ is the last positive single root of ($4.6$), then

Using ($4.4$), after removing $sin(\tau_1^* \rho_1^*)$ terms, we get

Using the analytic version of the implicit function theorem ^[32] to establish Hopf-bifurcation at $\tau = \tau_1^*$, we need to show that

Differentiate ($4.3$) with respect to $\tau$ for $\rho(\tau_1^*) = \rho_1^*$ and $\sigma(\tau_1^*) = 0$, and we get

Let $R = a_1+d\tau_1^*(b_1 \; cos (\tau_1^* \rho_1^*)+\rho_1^* \; sin (\tau_1^* \rho_1^*)) - d\; cos (\tau_1^* \rho_1^*)$

and $S = -2\rho_1^* + d\tau_1^*b_1\; sin(\tau_1^* \rho_1^*) - d\; sin (\tau_1^* \rho_1^*) - d\tau_1^* \rho_1^* \; cos(\tau_1^* \rho_1^*)$.

Equation ($4.10$) can be expressed as

Solving ($4.11$) for $\frac{d\sigma}{d\tau}$, we get

Using inequality $(4.7)$, we obtain

Figure 8 shows the Hopf-bifurcation at $\tau = \tau_1^*$.

5.   Eco-epidemiological model with no time delay

We now return to the full eco-epidemiological model with no time delay. Recall that the non-dimensionalized form of this is Eq (2.2):

5.1.   Boundedness

Propsition: The trajectory of system ($2.2$) is bounded.

Proof. Let $s = x+y+z$. Differentiating along the solution of model ($2.2$), we get

For any positive constant $P$, we have

since $a\geq b$ and $g\geq p$.

as $u_1+u_2 < k$ for all time implies that $x+y < 1$.

Now if $m > P$, then

Let $(1+P) = M$, therefore we have

Using the theory of differential inequality for $s(t)$, we get

where $s_0 = (x_0, y_0, z_0)$. It can be deduced that $l i m_{t\to\infty} sup\; s\leq \frac{M}{P}$, independently of the initial conditions. □

Corollary: If $\epsilon > 0$, then the region

is an invariant region for model ($2.2$).

5.2.   Positivity of Solutions

Since model ($2.2$) deals with animals populations, it is necessary to show that the prey and predator population will remain positive for all time.

Hence, all solutions $(x(t), y(t), z(t))$ will remain positive for all time with positive initial values, that is, $(x(0), y(0), z(0))\in \mathbb{R}^3_+$.

Before analyzing the model ($2.2$), we define two threshold quantities:

$(1)$ The infection reproduction number $r_0 = \frac{\beta}{c \mu}$. This is the product of the average infectious duration $\frac{1}{\mu}$ and $\frac{\beta}{c}$, the number of secondary infections per unit time at the disease-free equilibrium in a prey population with no predators. If $r_0 < 1$, disease will not invade a disease-free equilibrium in the absence of predators. If $r_0 > 1$, disease will invade such a population.

$(2)$ The predator reproduction number is $r_d = \frac{b}{m}$. This is the product of the predator life expectancy $\frac{1}{m}$ and the rate $b$ at which each predator reproduces in a prey population at equilibrium with no disease present. If $r_d < 1$, the predator cannot invade such a population, whereas if $r_d > 1$, the predator will invade such a population.

5.3.   Equilibrium points of the system

The possible equilibria for the system are:

(a) $E_0 = (0, 0, 0)$ is the trivial equilibrium and it always exists. Here, the populations of susceptible prey, infected prey, and predator species go to extinction.

(b) The two axial equilibria, $E_1 = (1, 0, 0)$ and $E_\theta = (\theta, 0, 0)$.

(c) The planar equilibrium $E_z = \left(\frac{1}{r_d}, 0, \overline{z}\right)$ where the infected prey is extinct. This exists when $1 < r_d < \frac{1}{\theta}$. In this case, $\overline{z}$ is the unique nonnegative root of the quadratic equation

By Descartes' rule of signs, Eq ($5.1$) has a unique nonnegative root if, and only if,

in which case $1 < r_d < \frac{1}{\theta}$.

(d) The planar equilibrium $E_y = (\overline{x}, \overline{y}, 0)$, where $\overline{x} = \frac{1}{r_0} \left(c+\delta \overline{y}\right)$ and $\overline{y}$ is the nonnegative root of the cubic equation

This equation can be written as $\gamma_1 y^3+\gamma_2y^2+\gamma_3y+\gamma_4 = 0$, where

Note that $\gamma_4 < 0$ if $1 < r_0 < \frac{1}{\theta}$.

(e) The interior equilibrium $E^* = (x^*, y^*, z^*)$, where

Note that

For this equilibrium to exist, we need

$y^*$ is a positive root of the equation $\xi _1 y^4+ \xi _2 y^3+\xi _3 y^2+\xi _4 y+ \xi _5 = 0$, where $\xi_1, \ \xi_2, \ \xi_3, \ \xi_4$, and $\xi_5$ are given in Appendix A.

5.4.   Stability analysis

It is simple to show that the stability matrix for the system ($2.1$) linearized about the equilibrium point is

5.4.1.   Behavior of the system around the equilibrium points $E_0, E_1$, and $E_\theta$

$E_0$ is locally asymptotically stable because all the corresponding eigenvalues are negative. The stability matrix for equilibrium $E_1 = (1, 0, 0)$ is:

The eigenvalues are $\lambda _1 = -1+\theta (<0)$, $\lambda _2 = \mu (r_0-1)$ and $\lambda _3 = m(r_d-1)$, so $E_1$ is locally asymptotically stable if $r_0 < 1$ and $r_d < 1$, as shown in Figure 9.

The stability matrix for equilibrium $E_\theta = (\theta, 0, 0)$ is:

The eigenvalues are $\lambda _1 = \theta(1-\theta)$, $\lambda _2 = \mu \theta \left(r_0-\frac{1}{\theta}\right)$ and $\lambda _3 = m\theta\left(r_d-\frac{1}{\theta}\right)$. $E_\theta$ is unstable since $\lambda _1 = \theta(1-\theta) > 0$.

5.4.2.   Behavior of the system around the equilibrium point $E_z=\left(\frac{1}{r_d}, 0, \overline{z}\right)$

The stability matrix for $E_z$ is:

One eigenvalue is given by $\lambda_2 = \frac{\beta}{cr_d}-g\overline{z}-\mu$, which is negative if

The other eigenvalues, $\lambda_1$ and $\lambda_3$, are roots of the characteristic equation given by

Clearly, the roots of above equation have negative real parts if

Therefore, the equilibrium $E_z$ is locally asymptotically stable if

and $r_d < \frac{2}{1+\theta}$. Figure 10 shows the local stability for the equilibrium point $E_z$.

5.4.3.   Behavior of the system around the equilibrium point $E_y=(\overline{x}, \overline{y}, 0)$.

The stability matrix for $E_y$ is:

where $\tau_1 = \overline{x} (1-2 \overline{x} - \overline{y} - \theta)$, $\tau_2 = -\overline{x}(\overline{x} - \theta) - \frac{\beta \overline{x} c}{(c + \delta \overline{y})^2}$ and $\tau_3 = \frac{\overline{x} (1-\overline{x}-\overline{y})(\overline{x}-\theta)(1-\eta)}{\alpha} - a \overline{x}$. One eigenvalue is given by $\lambda_3 = b\overline{x}-p \overline{y} -m$. The other eigenvalues, $\lambda_1$ and $\lambda_2$, are roots of the characteristic equation given by:

Here, $\ K = \tau_1, \ L = \tau_2 < 0, \ M = \frac{\beta \overline{y}}{c+\delta \overline{y}} > 0, \ \text{and} \ N = - \frac{\delta \beta \overline{x} \overline{y}}{(c + \delta \overline{y})^2} < 0.$

The roots of the characteristic equation have negative real parts if $K+N < 0$ and $(KN-ML) > 0$. Figure 11 shows the solution converges to the equilibrium point $E_y = (\overline{x}, \overline{y}, 0)$.

5.4.4.   Behavior of the system around the coexistence equilibrium point $E^*=(x^*, y^*, z^*)$

The stability matrix at the interior attractor is given by

where $A = \left(\eta +\frac{ 1-\eta}{1+s_\alpha z^*}\right)[x^*(1-2x^*-y^*+\theta)], \ B = -x^*(x^*-\theta)\left(\eta+\frac{ 1-\eta}{1 +s_\alpha z^*}\right)+\frac{\delta \beta x^*y^*}{(c+\delta y^*)^2}-\frac{\beta x^*}{c+\delta y^*} < 0$, $C = \frac{-x^*(1-x^*-y^*)(x^*-\theta)\alpha (1-\eta)}{(\alpha+z^*)^2}-ax^* < 0, \ \ D = \frac{\beta y^*}{c+\delta y^*} > 0$, $E = -\frac{\delta \beta x^* y^*}{(c+\delta y^*)^2} < 0$, $F = -gy^* < 0$, $G = bz^* > 0$, and $H = p z^* > 0.$

The characteristic equation associated with this stability matrix is given by

where $a_1 = tr[J_E^*] = -(A+E), a_2 = AE-BD-FH-CG$, and $a_3 = -det[J_E^*] = AHF-BGF-DHC+GEC$. From the Routh-Hurwitz criteria, the equilibrium $E^*$ is locally asymptotically stable if, and only if, $a_1 > 0$, $a_3 > 0$, and $a_1a_2 > a_3$. Figure 12 demonstrates that the solution converges to the equilibrium point $E^* = (x^*, y^*, z^*)$.

5.5.   Bifurcation analysis

When the parameters go through specific critical values, a transcritical bifurcation and a Hopf-bifurcation occur around the equilibrium points.

5.5.1.   Transcritical bifurcation

According to Sotomayor's theorem for local bifurcation ^[31], if the following conditions are satisfied:

$(1) \; w^T f_m(E_1, m_0) = 0$,

$(2) \; w^T (Df_m(E_1, m_0). \upsilon)\neq 0$,

$(3) \; w^TD^2f(E_1, m_0)(\upsilon, \upsilon)\neq 0,$

where $w$ and $\upsilon$ are the left and right eigenvectors of $f_m(E_1, m_0)$, then the system ($2.2$) experiences a transcritical bifurcation at the equilibrium point $E_1 = (1, 0, 0)$, when $r_d = 1$ (i.e $m = b$).

Theorem 4. If $r_0 \ne c$, then the system $(2.2)$ undergoes a transcritical bifurcation at the positive equilibrium $E_1 = (1, 0, 0)$ when $m = b$.

Proof:

The linearized system around the equilibrium point $E_1$ is:

Now, $\lambda_3 = 0$ at $r_d = 1, \ $(i.e, at $m_0 = b)$. So

Let us define $\upsilon = (\upsilon_1, \upsilon_2, \upsilon_3)^T$ and $w = (w_1, w_2, w_3)^T$ to be the right and left eigenvectors of $\lambda_3 = 0$. Now solving $J(E_1, m_0)\upsilon = 0$ implies that $\upsilon = \bigl(\frac{a\upsilon_3}{-1+\theta}, 0, \upsilon_3\bigr)$, and by solving $J^T(E_1, m_0)w = 0$, we get $w = (0, 0, w_3)$ where $\upsilon_3, w_3$ are any nonzero real numbers.

Now, if $X = (x, y, z)$, system $(2.2)$ can be rewritten in the following vector form:

Taking the derivative of $f(X)$ with respect to $m$, we get $f_{m}(E_1, m_0) = \begin{pmatrix} 0\\ 0\\ 0 \end{pmatrix}$. Hence, $w^T f_m(E_1, m_0) = 0$.

Next, taking the derivative of $f_m(X)$ with respect to $X = (x, y, z)^T$, we get

Therefore, we have $w^T (Df_m(E_1, m_0). \upsilon) = -w_3\upsilon_3\neq 0$. Furthermore,

where $A_1 = \left(\eta+\frac{ 1-\eta}{1 +s_\alpha z}\right)[-2(x-\theta)-2x]$, $A_2 = \left(\eta+\frac{ 1-\eta}{1 +s_\alpha z}\right)[-(x-\theta)-x]-\frac{\beta c}{(c+\delta y)^2}$, $A_3 = - \frac{\alpha(1-\eta)}{(\alpha+z)^2}[(1-x-y)(2x-\theta)-x(x-\theta)]-a$, $A_4 = \left(\eta+\frac{ 1-\eta}{1 +s_\alpha z}\right)[-(x-\theta)-x]-\frac{\beta c}{(c+\delta y)^2}$, $A_5 = \frac{2 \beta c \delta y}{(c+\delta y)^3}$, $A_6 = \frac{x(x-\theta)\alpha(1-\eta)}{(\alpha+z)^2}$, $A_7 = \frac{\alpha(1-\eta)}{(\alpha+z)^2}[-(1-2x-y)(x-\theta)-x(1-x-y)]-a$, $A_8 = \frac{x(x-\theta)\alpha(1-\eta)}{(\alpha+z)^2}$, $A_9 = \frac{2x(1-x-y)(x-\theta)\alpha(1-\eta)}{(\alpha+z)^3}$, $A_{10} = 0$, $A_{11} = \frac{ \beta c}{(c+\delta y)^2}$, $A_{12} = 0$, $A_{13} = \frac{\beta c y}{(c+\delta y)^2}$, $A_{14} = -\frac{2\delta c \beta x}{(c+\delta y)^3}$, $A_{15} = -g$, $A_{16} = 0$, $A_{17} = -g$, $A_{18} = 0$, $A_{19} = 0$, $A_{20} = 0$, $A_{21} = b$, $A_{22} = 0$, $A_{23} = 0$, $A_{24} = p$, $A_{25} = b$, $A_{26} = p$, $A_{27} = 0$.

So, $w^TD^2f(E_1, m_0)(\upsilon, \upsilon) = 2b\upsilon_1\upsilon_3 w_3+2p\upsilon_2\upsilon_3 w_3 \neq 0$. Therefore, according to Sotomayor's theorem ^[31] for local bifurcation, system $(2.2)$ has a transcritical bifurcation at steady-state $E_1 = (1, 0, 0)$ when the parameter $r_d = 1$, (i.e, $m = b$) (Figure 13).

5.5.2.   Hopf-bifurcation analysis

We observe that the system shows the Hopf-bifurcation for the equilibrium point $E_z$ with respect to the parameter $\theta$ when $\theta = 0.333335$, $(E_z = (0.666667, 0, 0.046177))$, as shown in Figure 14 by using Matcont software for bifurcation. The blue color shows that the equilibrium is stable, while the red color shows the unstable equilibrium.

In addition, Figure 15 shows another Hopf-bifurcation for the equilibrium point $E_z$ with respect to the parameter $m$ when $m = 0.900001$, $(E_z = (0.600000, 0, 0.062795))$.

Also, we observed the Hopf-bifurcation for the equilibrium point $E_z$ with respect to the parameter $b$ when $b = 1.666664$, $(E_z = (0.600001, 0, 0.062795))$, as seen in Figure 16.

Figure 17 shows the Hopf-bifurcation for the equilibrium point $E_y$ with respect to the parameter $\mu$ when $\mu = 0.845893 \; (E_y = (0.502786, 0.075836, 0))$.

Also, we found the Hopf-bifurcation for the equilibrium point $E_y$ with respect to the parameter $\delta$ when $\delta = 0.101810 \; (E_y = (0.556617, 0.073434, 0))$, as shown in Figure 18.

5.6.   Effect of cost of fear $s_\alpha$ and the Allee threshold effects $\theta$

We observe from Figure 19 that the population of uninfected prey declines as fear intensity rises, and after a certain level of fear, population of uninfected prey will not change due to habituation to stress. Increasing the fear parameter stabilizes the system. For an intermediate value of fear, the system becomes unstable, and after a certain level of fear, the system shows instability and the predator population goes extinct. As a result, the infected prey and the uninfected prey interact and this increases disease transmission, and, consequently, the infected prey increases.

Figure 20 shows the effect of the Allee threshold effect on uninfected prey, infected prey, and predator. It is clear that when the Allee threshold increases, the population of uninfected prey and predator decrease while the population of infected prey increases. When the Allee threshold exceeds $0.3$, all populations go to extinction.

6.   Full eco-epidemiological model with a time delay

Recall the non-dimensionalized full eco-epidemiological model ($2.4$) with a time delay included:

Also define $b\overline{x}+p\overline{y} = m$. The characteristic equation for model ($2.4$) is

If $\tau = 0$, then $E^*$ is asymptotically stable under certain conditions. When $\tau\neq 0$, we take $\xi = \iota+i\psi$ where $\iota, \psi\in R$. Substituting the value of $\xi$ into $(6.1)$ and separating in terms of real and imaginary parts:

Let $\tau_1 = \tau^*_1$ be such that $\iota(\tau^*_1) = 0$ and $\psi(\tau^*_1) = \psi^*_1$, then ($6.2$) can be written as

where $A_1 = Q_2\psi^{*2}_1-Q_0$ and $B_1 = Q_1\psi^*_1$. By squaring and adding, on simplification, it follows that

Let $\psi^{*2}_1 = w$. Also

where $P^*_0 = P^2_0-Q^2_0, \ P^*_1 = -2P_2P_0+P^2_1+2Q_0Q_2-Q^2_1$, and $P^*_2 = P^2_2-2P_1-Q^2_2$. Note that $\varphi\prime (w) = 3w^2+2P^*_2w+P^*_1$. The two roots of the equation $\varphi\prime (w) = 0$ are: $w_1 = \frac{-P^*_2+\sqrt{P^{*2}_2-3P^*_1}}{3}$ and $w_1 = \frac{-P^*_2-\sqrt{P^{*2}_2-3P^*_1}}{3}$. The conditions for ($6.5$) to have positive roots are as follows:

(a) If $P^*_0 < 0$, then ($6.5$) has at least one positive root.

(b) If $P^*_0\geq0$ and $P^{*2}_2\leq3P^*_1$, then ($6.5$) has no positive root.

(c) If $P^*_0\geq 0$ and $P^{*2}_2 > 3P^*_1$, then ($6.5$) has a positive root if, and only if, $w_1 > 0$ and $\phi(w_1)\leq 0$.

Using the above conditions, we can say that if condition $(b)$ is satisfied then the stability of $E^*$ will not change on increasing $\tau$. If $(6.5)$ has a positive root, then the stability of $E^*$ may change as $\tau$ changes. Let $w_1, w_2$, and $w_3$ be three positive roots of $(6.5)$, then $(6.4)$ has positive roots $\psi^*_1 = \sqrt{w_k}$, $k = 1, 2, 3$. Using $(6.3)$, we get:

where $\psi^*_1$ is the last positive root of ($6.4$). To establish Hopf-bifurcation at $\tau = \tau^*_1$, we need to show that

Differentiating ($6.2$) with respect to $\tau$ and setting $\iota = 0$ and $\psi = \psi^*_1$, we get

where $U = Q_2 \psi^{*2}_1\tau^*_1-Q_0\tau^*_1+Q_1, \ V = 2Q_2\psi^*_1-Q_1\psi^*_1\tau^*_1, \ M = -3\psi^{*2}_1+P_1, \ N = 2 P_2 \psi_1^*, \ X = Q_2\psi^{*3}_1-Q_0\psi^*_1$, and $Y = Q_1(\psi^*_1)^2$. Finally, by solving $\frac{d\iota}{d\tau}\Big\vert_{\tau = \tau^*_1}$ and $\frac{d\psi}{d\tau}\Big\vert_{\tau = \tau^*_1}$, we have

where $S_1 = Ucos(\tau^*_1 \psi^*_1)+V sin(\tau^*_1 \psi^*_1)+M, \ S_2 = Usin(\tau^*_1 \psi^*_1)-V cos(\tau^*_1 \psi^*_1)-N$, and $S^2_1+S^2_2 > 0$.

If $\psi^*_1$ is the last positive root of ($6.4$), then using ($6.5$),

Figure 21 shows the Hopf-bifurcation of model $(2.4)$ with respect to $\tau$.

6.1.   Maximum delay to maintain stability

Let $X(t) = \overline{X}+u(t)$, $Y(t) = \overline{Y}+v(t)$, and $Z(t) = \overline{Z}+w(t)$. Linearizing the system ($2.4$) at the coexistence equilibrium $E^*$, we get:

Taking Laplace transforms for ($6.11$), we obtain

where $L_2 = \int_{-\tau}^0 e^{-s(t_1)}v(t_1)dt_1 = \int_{-\tau}^0 e^{-s(t)}v(t)dt$ and $L_3 = \int_{-\tau}^0 e^{-s(t_1)}w(t_1)dt_1 = \int_{-\tau}^0 e^{-s(t)}w(t)dt$.

If $L\left\lbrace u(t) \right\rbrace, L\left\lbrace v(t) \right\rbrace$ and $L\left\lbrace w(t) \right\rbrace$ have poles with positive real parts, then the inverse Laplace transformation of $L\left\lbrace u(t) \right\rbrace, L\left\lbrace v(t) \right\rbrace$, and $L\left\lbrace w(t) \right\rbrace$ will have terms which exponentially increase with time. Thus, $E^*$ will be locally asymptotically stable if, and only if, all poles of $L\left\lbrace u(t) \right\rbrace, L\left\lbrace v(t) \right\rbrace$, and $L\left\lbrace w(t) \right\rbrace$ have negative real parts. Using the Nyquist criterion $E^*$ will be asymptotically stable if the following two conditions are satisfied ^[33,34]:

where $G(s) = s^3+P_2s^2+P_1s+P_0+e^{-s\tau}(Q_2s^2+Q_1s+Q_0) = 0$ and $w_0$ is the smallest positive root of ($6.13$). Equating the real and imaginary parts, we get

To find the length of maximum delay so that stability is preserved, we write $(6.15)$ and ($6.16$) as

Equation $(6.17)$ and inequality $(6.18)$ are the conditions for stability and both should be satisfied simultaneously. We first find an upper bound $w^+$, independent of $\tau$, such that ($6.18)$ is valid for $w_0 \le w \le w^+$. We will later use this inequality to determine a range of feasible values for $\tau^*$.

Since $\mid cos(w\tau)\mid\leq 1$ and $\mid sin(w\tau)\mid\leq 1$, using these inequalities in $(6.17)$,

Let $w_+$ denote the maximum value of $w$ satisfying (6.18),

where $w_+\geq w_0$. From ($6.18$), we get

Since $E^*$ is locally asymptotically stable for $\tau = 0$, by assumption, the inequality ($6.22$) will continue to hold for sufficiently small $\tau$ and $w = w_0$. Using $(6.17)$ and $(6.22)$, we get

where $sin(w\tau)\leq w\tau$ and $1-cos(w\tau)\leq 2 sin^2(\frac{w\tau}{2})\leq \frac{w^2 \tau^2}{2}$. So

since $w_0\leq w\leq w_+$.

Let $\frac{\mid Q_2w^2_+ + P_2Q_1-Q_0\mid w_+^2}{2} = L_1$, $\mid Q_1-P_2 Q_1\mid w^2_++\mid P_2Q_0\mid = L_2$, and $P_2P_1-P_0+P_2Q_1-Q_0+Q_2w_+^2 = L_3$. So, $L_1 \tau^2+L_2\tau \leq L_3$. Suppose $\tau^*$ is the positive root of $L_1 \tau^2+L_2\tau = L_3$, that is,

Theorem 5. If $0 < \tau < \tau^*$ and the Nyquist criterion to the local asymptotic stability is satisfied, then $\tau^*$ estimates the maximum delay length to preserve stability.

7.   Conclusions

This paper has examined prey-predator models with the Allee effects and the cost of fear in prey reproduction, where the predator predates the prey with a Holling Type I functional response. We have investigated the dynamic behavior of the model mathematically, including the existence and stability of equilibria, the occurrence of Hopf-bifurcation around the coexistence equilibrium point, and the existence of a limit cycle that emerges from Hopf-bifurcation. We have described equilibria, the stabilities of the system, and the occurrence of Hopf-bifurcation. The coexistence equilibrium is locally asymptotically stable if twice the death rate of the predator is one unit higher than the Allee threshold value.

We have shown that as the fear level rises, the population of prey remains unaffected because, at $N^* = d$, the fear level has no effect. However, when the level of fear increases, the population of predators decreases. The existence of the coexistence equilibrium point $E^* = (N^*, P^*)$ implies that $\theta < d < 1$. Since $N^* = d$, the Allee effect has no affect on the prey population, and the coexistence equilibrium $E^*$ is stable when it occurs for a small Allee effect. Furthermore, we have demonstrated that as the Allee threshold rises, the predator population rises. However, after certain values of the Allee threshold, both populations become unstable and then go extinct. If we take the death rate of the predator $d$ as a bifurcation parameter, the system experiences Andronov-Hopf bifurcation. When the parameter $d$ passes through the bifurcation value $d_0 = 1$, our system experiences a transcritical bifurcation at steady state $E_1$ according to Sotomayor's theorem ^[31]. The dynamics of the delayed system have been investigated and we have discussed the Hopf-bifurcation.

Finally, we have looked at an eco-epidemiological model that takes the cost of fear and the Allee effect into account in the prey population. We have examined a predator-prey model in which the prey is afflicted with a disease. The overall prey population is divided into two categories, susceptible prey and infected prey. When susceptible prey comes into contact with infected prey, it becomes infected. Infected prey is unable to procreate because it has the disease. The predator consumes both uninfected and infected prey, but infected prey is easy to capture whilst it takes more time to catch the uninfected prey. As a result, we have assumed that the consumption rate of infected prey is higher than that of uninfected prey. There are five biologically significant equilibria. We have discussed the stabilities of these equilibria. We have performed a numerical study to examine how the Allee effect and fear effect affect our system, and we have discovered that the population of uninfected prey and predator drops as the Allee threshold rises, whereas the population of infected prey rises. All populations become extinct when the Allee threshold reaches a specified value. We have demonstrated that when the fear intensity rises, the population of uninfected prey declines, and after a certain level of fear, there is no effect on the population of uninfected prey as the fear increases, and this is due to physiological effects when they get acclimatised to fear. The system becomes stable for an intermediate level of fear, and after a certain level of fear, the system exhibits instability because the predator population goes extinct, leaving only uninfected and infected prey, which increases disease transmission and leads to an increase in infected prey. We have found that the equilibria $E_z$ and $E_y$ have a Hopf-bifurcation under specific parametric conditions.

To make the model ($2.1$) more realistic, we have incorporated a time delay and Hopf-bifurcation occurs by taking the delay as the bifurcation parameter. Hopf-bifurcation will occur when the delay $\tau$ passes through a series of critical values. If $0 < \tau < \tau^*$ and the Nyquist criterion to the local asymptotic stability is true then $\tau^*$ estimates the maximum delay length to preserve stability.

Hopf-bifurcation has helped us in finding the existence of a region of instability in the neighborhood of non-zero equilibrium, where prey and predator species will persist undergoing regular fluctuation. For the transcritical bifurcations a stable coexistence equilibrium $E^*$ combines with another equilibrium and disappears and the positive equilibrium $E_1$ switches from instability to stability.

Use of AI tools declaration

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

Data availability

Data sharing is not applicable to this article as no datasets were generated or analysed during the current study.

Conflict of interest

The authors declare that there is no conflict of interest.

Appendix A