1. Introduction

The study of prey-predator interactions has been an important issue in mathematical modeling and hence received considerable attention from the researchers for the past few decades. After the classical work of Lotka [11] and Volterra[19] significant progress has been made both in modeling approach and their mathematical study. According to several previous works it is now an established fact that applicability of the seminal theory based on the Lotka-Volterra formulation is limited to a certain class of ecological problems as it mostly hints that if the population density is small then the intraspecific competition between the species will be less. As discussed in [16] we see that these classical modeling approach can not capture phenomena such as predator saturation, group defense etc. In fact the Lotka-Volterra system fails to address many important processes and one of them is Allee-effect, which was first reported by W. C. Allee in 1931 [3]. Allee effect mainly signifies a decrease in per-capita growth rate at low population density which in turn refers that individual fitness is directly proportionate to population density. Allee effect can be caused by difficulties in mate finding, social dysfunction at small population size, inbreeding depression, food exploitation, and predator avoidance or defense etc [7,8,16]. During the last decades, several works can be found in the literature addressing the issue of Allee effect [1,2,5,6,12,14,15]. Before going into further details it is important to note that there can be two types of Allee effects namely the strong and the weak Allee effects respectively. The strong Allee effect is subjected to a threshold below which the population growth becomes negative [8,10]. But in the case of weak Allee effect the population growth decreases but remains positive even in small population size [4,7,16].

In [6] and [14] the authors have studied models with Allee effect in the prey growth where the functional responses were Holling Type-Ⅱ and ratio-dependent respectively. The authors have reported rich dynamics around the non trivial equilibrium points, such as Hopf-bifurcation, Bogdanove-Takens bifurcation, disappearance of limit cycle through homoclinic bifurcations etc. for both the systems. Specially the model considered in [14] exhibits very rich dynamics around the trivial equilibrium point. Gonzalez et al have shown that incorporation of different types of Allee effect can contribute to a significant change in the system dynamics [1,2]. In fact they studied the existence of two and three limit cycles around the nontrivial equilibria in [2] and [1] respectively. During 2015, in [15] Sen et al considered a prey predator model where the predator was subjected to intraspecific competition. Through a complete bifurcation analysis the authors have reported the existence of cusp-bifurcation, homoclinic, heteroclinic bifurcations, generalized Hopf bifurcation, along with other usual bifurcations such as transcritical, saddle-node and Hopf bifurcations due to the presence of Allee effect in the prey growth function.

Thus from the literature it is well understood that the incorporation of Allee effect in the modeling approach reflects and can counter the mechanisms such as animal aggregation/grouping, predator invasion, intraspecific competition etc and hence indeed increases the scope to capture the dynamics of a wider range of species. Also going through the above works one important conclusion can be drawn is that the presence of Allee effect in the prey growth for two-dimensional prey-predator model prevents the appearance of large amplitude limit cycles due to the enrichment at basic tropic level. Since in almost all the works reported till date on Allee effect the researchers have considered two dimensional systems how the nature of the chaotic dynamics would be affected was beyond the scope of the study unless we consider diffusion or stochasticity in the system. But to capture the dynamics of complex natural systems two dimensional ecological systems are merely adequate. Hence the most important and ecologically relevant question is whether the inclusion of Allee effect to one or more tropic levels for three and higher dimensional food-chain and interacting population models can alter the scenario of chaotic oscillations or not. This question is relevant for the reason that the occurrence of chaos in three and higher dimensional population models are theoretical outcomes of the dynamics for the concerned model but not agreed by every ecologist due to the lack of support from field data.

Also from the mathematical point of view many works in literature on three or higher dimensional prey-predator systems can be found, exploring very rich dynamics such as local and global bifurcations, different types of chaos etc. Therefore it is very important to address if the introduction of Allee effect in any way affects the dynamical complexities for three or higher dimensional prey-predator systems.

In this work especially we are interested to see how Allee effect affects the chaotic dynamics of a three dimensional prey-predator system where the predator is generalized and both the prey species are subjected to Allee effect in their growth functions. Here we present a complete study of the system through stability and bifurcation analysis. To emphasize on how the incorporation of Allee effect suppresses or enhances the chaotic behavior of the corresponding analogous model [18] without Allee effect, we focus on the two parametric bifurcation diagrams taking the two Allee effects as the bifurcation parameters. The paper is organized as follows. In the next section the model has been proposed with a small discussion on the motivation. The subsequent section addresses the study of the equilibria and their stability behaviour. Next we present a complete study on the bifurcations that the system undergoes. Finally we focus on the bifurcation diagrams and their significance in the context of the underlying system parameters and finally the conclusion in the

---

2. Model formulation

Takeuchi et al. [18] have studied a three dimensional two prey one predator model in [18], where the predator is, generalist i.e. survived on the two prey populations. The authors have worked with the model given by,

$\frac{dN_1}{dt} = N_1(a_{11}-N_1-a_{12} N_2-m_1 P),$

(1a)

$\frac{dN_2}{dt} = N_2(a_{22}-a_{21} N_1- N_2-m_2 P),$

(1b)

$\frac{dP}{dt} = P(-d_3+e m_1 N_1+em_2 N_2).$

(1c)

with, $N_1(0)>0, N_2(0)>0, P(0)>0$. Here, $a_{ii}(i = 1, 2)$ are the intrinsic growth rate of the two preys $N_1, N_2$ respectively. $a_{12}\, \, \&\, \, a_{21}$ represent the coefficient of competition of $N_1,$ and $N_2$, while the parameters $m_1$ and $m_2$ signify the decrease of $N_1$ and $N_1$ due to predation by the generalist predator $P$. The predator death rate and an equal transformation rate of predator to the remaining $N_1$ and $N_2$ are respectively denoted by $d_3$ and $e$. In the above mentioned work the authors have considered the two prey's subjected to logistic growth rates and the exponential type functional response. Even with such simplest modeling approach the system exhibits rich dynamical behaviour such as bistability of equilibira, Hopf-bifurcation, period doubling chaos etc. Hence it is a very natural question that how the introduction of Allee effect will change the system dynamics. Influenced by these ideas we pose a generalist prey-predator model where the prey growth rates are subjected to Allee effects. But before we introduce the final model we would like to give a brief discussion on the motivation for the incorporation of particular type of Allee effects in the present work.

After the work of W. C. Allee[3] the literature has been enriched by several works where the Allee effect has been addressed and incorporated in the modelling. Evidently in most of the works the system with Allee effect exhibits rich dynamics than that without the Allee effect [1,2,5,6,12,14,15]. Surprisingly in those earlier works authors have mainly considered the Allee effect in the prey growth function in two-dimensional prey-predator system. Recently there are few article available where the authors have discussed the effects of Allee effect in the predator growth rate. But hardly any work in literature can be found where generalist predator-prey system with Allee effect in both the preys is considered. One more thing comes up, when we go through the literature, is that the most common way of incorporating Allee effect is multiplicative Allee effect. But many researchers have shown that the additive Allee effect also shows complex dynamics [1,2,21]. In [21], the authors have incorporated the additive Allee effect in a two dimensional population model in a different but interesting way. Firstly they have considered the single species population growth with logistic type growth as

$\frac{dN}{dt} = N(b-d-\alpha N).$

(2)

where, $b$ is the per capita maximum fertility rate of the population, $d$ is the per capita death rate and $\alpha$ denotes the strength of intra-competition. Then the authors have incorporated an additive Allee effect in the prey-growth as follows.

$\frac{dN}{dt} = N\left(\frac{bN}{A+N}-d-\alpha N\right),$

(3)

where $A$ is the strength of the Allee effect. A similar formalism was also adopted in [20].

Based upon these aforesaid literature, we now propose the following model:

$\frac{dN_1}{dt} = N_1\left(\frac{b_1N_1}{N_1+A_1}-d_1-k_1N_1\right)-a_{12} N_1N_2-m_1 N_1P,$

(4a)

$\frac{dN_2}{dt} = N_2\left(\frac{b_2N_2}{N_2+A_2}-d_2-k_2N_2\right)-a_{21} N_1N_2-m_2 N_2P,$

(4b)

$\frac{dP}{dt} = P(-d_3+em_1N_1+em_2N_2).$

(4c)

After the transformation $N_1 = \frac{d_1}{k_1} x, N_2 = \frac{d_1}{k_2} y$, $P = \frac{d_1}{m_2} z, t = \frac{1}{d_1} T$, the dimensionless model is given by,

$\frac{dx}{dt} = x\left[\frac{\beta_1 x}{x+\alpha_1}-1-x-\alpha y-\epsilon z\right],$

(5a)

$\frac{dy}{dt} = y\left[\frac{\beta_2 y}{y+\alpha_2}-\gamma-\beta x-y-z\right],$

(5b)

$\frac{dz}{dt} = z[-\beta_3+d \epsilon x+d \mu y],$

(5c)

where $\alpha_1 = \frac{k_1}{d_1}A_1$, $\beta_1 = \frac{b_1}{d_1}$, $\alpha = \frac{a_{12}}{k_2}$, $\epsilon = \frac{m_1}{m_2}$, $\alpha_2 = \frac{k_2}{d_1}A_2$, $\beta_2 = \frac{b_2}{d_1}$, $\beta = \frac{a_{21}}{k_1}$, $\gamma = \frac{d_2}{d_1}$, $\beta_3 = \frac{d_3}{d_1}$, $d = \frac{e m_2}{k_1}$, $\mu = \frac{k_1}{k_2}$.

---

2.1. Positivity and boundedness:

Firstly we discuss the boundedness and positivity of solutions of the system (5) starting from positive initial conditions. From (5a) it is clear that the solution $x(t)\rightarrow 0$ as $t\rightarrow \infty$ for any $x(0)> 0$ when $\beta_1\leq1$. Similarly, $y(t)\rightarrow 0$ as $t\rightarrow \infty$ when $\beta_2\leq\gamma$. Hence we assume that $\beta_1>1$ and $\beta_2>\gamma$ hereafter.

The equation (5) can be written as,

$x(t) = x(0)exp\Big[\int_0^t\Big(\frac{\beta_1 x(s)}{x(s)+\alpha_1}-1-x(s)-\alpha y(s)-\epsilon z(s)\Big)ds\Big], \\ y(t) = y(0)exp\Big[\int_0^t\Big(\frac{\beta_2 y(s)}{y(s)+\alpha_2}-\gamma-\beta x(s)-y(s)-z(s)\Big)ds\Big], \\ z(t) = z(0)exp\Big[\int_0^t\Big(-\beta_3+d \epsilon x(s)+d \mu y(s)\Big)ds\Big],$

showing that $x(t)\geq0, \, \, y(t)\geq0$ and $z(t)\geq0$ whenever the initial conditions are all positive. Hence all solutions remain within the first quadrant of the xyz-plane starting from an interior point.

Next let us set, $W = d x+d \mu y+z$ and choose $0<\lambda<\beta_3$. Then after some simple calculations from (5) we have,

$\dot{W}+\lambda W\leq d x(\beta_1+\lambda-x)+d\mu y(\beta_2+\lambda-y)+(\lambda-\beta_3)z.$

Both $d x(\beta_1+\lambda-x)$ and $d\mu y(\beta_2+\lambda-y)$ are quadratic functions, $x$ and $y$ are non-negative variables, and bounded by, $\frac{d(\beta_1+\lambda)^2}{4}$ and $\frac{d\mu(\beta_2+\lambda)^2}{4}$ respectively. Hence,

$\dot{W}+\lambda W \leq \frac{d(\beta_1+\lambda)^2}{4}+\frac{d\mu(\beta_2+\lambda)^2}{4}+(\lambda-\beta_3)z.\\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\leq \frac{d(\beta_1+\lambda)^2}{4}+\frac{d\mu(\beta_2+\lambda)^2}{4} = M (\textrm{say}),$

where we have used the positivity of $z$ and the parametric restriction $0<\lambda<\beta_3$. The above inequality implies $W$ is bounded above and hence the boundedness of $x$, $y$ & $z$.

---

3. Equilibria and their stability:

In this section we focus on the existence of various equilibrium points and their local stability by analyzing the eigenvalues of the Jacobian matrix evaluated around the equilibrium points.

Proposition 1.   The model (5) admits trivial equilibrium points denoted by $E_0(0, 0, 0)$ which is always locally asymptotically stable.

Proof.   Trivial.

Proposition 2.   Assume $\beta_1\geq(1+\sqrt{\alpha_1})^2$ then,

(a) the model (5) can have at most two axial equilibria of the form $E_1^+(x_1^+, 0, 0)$ and $E_1^-(x_1^-, 0, 0)$ respectively for which the second prey and the predator population are absent.

(b) $E_1^-$ is always a saddle point whenever exists. $E_1^+$ is locally asymptotically stable if $x^+<\frac{\beta_3}{d \epsilon}$ otherwise it is a saddle point. In particular if we consider $\alpha_1$ to be the varying parameter then $E_1^+$ is always locally asymptotically stable if $\beta_1<\left(1+\frac{\beta_3}{d\epsilon}\right)$ or we will always have a range of values of $\alpha_1$ for which $E_1^+$ is locally asymptotically stable if $\left(1+\frac{\beta_3}{d\epsilon}\right)<\beta_1<\left(1+\frac{\beta_3}{d\epsilon}\right)^2$.

Proof.  We consider the equilibrium point for which second prey and predator population are absent. Substituting $y = 0$ and $z = 0$ in (5) we find the following quadratic equation,

$\frac{\beta_1 x}{x+\alpha_1}\, = \, 1+x\, \, \Rightarrow\, \, x^2+(\alpha_1+1-\beta_1)x+\alpha_1\, = \, 0.$

(6)

The constant term is positive and hence we can find at most two positive real roots leading to at most two axial equilibrium points, on $x$-axis, which will be denoted by $E_1^+(x_1^+, 0, 0)$ and $E_1^-(x_1^-, 0, 0)$ respectively with $x_1^-<x_1^+$. As the constant term is positive, both the roots will be of same sign and hence the possibility of the case $x_1^-<0<x_1^+$ is excluded. Two roots of the equation are given by,

$x^{\pm}_1\, = \, \frac{-(1+\alpha_1-\beta_1)\pm \sqrt{(1+\alpha_1-\beta_1)^2-4\alpha_1}}{2}.$

Thus two roots are real and positive provided $(1+\alpha_1-\beta_1)^2-4\alpha_1>0$ and $1+\alpha_1-\beta_1<0$. Note that $(1-\alpha_1-\beta_1)^2-4\alpha_1 = (\beta_1-(\sqrt{\alpha_1}+1)^2) (\beta_1-(\sqrt{\alpha_1}-1)^2)>0$ implies that $\beta_1>(\sqrt{\alpha_1}+1)^2$ or $\beta_1<(\sqrt{\alpha_1}-1)^2$, which lead to the following conclusions,

$\bullet$  Case1. $\beta_1<(1+\sqrt{\alpha_1})^2$.

In this case (6) has no positive root and the $E_1$ does not exist.

$\bullet$  Case2. $\beta_1>(1+\sqrt{\alpha_1})^2$.

In this case (6) has two distinct positive roots $x^+_1$ and $x^-_1$. Hence (5) has two axial equilibria $E^+_1$ and $E^-_1$.

$\bullet$  Case3. $\beta_1 = (1+\sqrt{\alpha_1})^2$.

In this case (6) has a unique axial equilibrium point $E_1$.

Next we evaluate the Jacobian at $E_1$ and which is given by,

$J(E_1) = \begin{pmatrix} x_1\left[\frac{\alpha_1\beta_1}{(x_1+\alpha_1)^2}-1\right]&-\alpha x_1&-\epsilon x_1 \\ 0&-\gamma-\beta x_1&0 \\ 0&0&-\beta_3+d\epsilon x_1 \end{pmatrix}.$

After some algebraic calculations we have $\left[\frac{\alpha_1\beta_1}{(x_1^-+\alpha_1)^2}-1\right]>0, \, \, \&\, \, \left[\frac{\alpha_1\beta_1}{(x_1^++\alpha_1)^2}-1\right]<0.$ Hence $E^-_1$ is always a saddle point and $E^+_1$ is locally asymptotically stable if $x^+_1<\frac{\beta_3}{d \epsilon}$ and a saddle point if $x^+_1>\frac{\beta_3}{d \epsilon}$. Further let us consider $\alpha_1$ to be the varying parameter and

$\lambda_1^{E_1^+} = x^+_1-\frac{\beta_3}{d \epsilon} = \frac{-(1+\alpha_1-\beta_1)+ \sqrt{(1+\alpha_1-\beta_1)^2-4\alpha_1}}{2}-\frac{\beta_3}{d \epsilon}, \\ \frac{d\lambda_1^{E_1^+}}{d\alpha_1} = \frac{1}{2}\left(-1+\frac{\alpha_1-\beta_1-1}{\sqrt{(1+\alpha_1-\beta_1)^2-4\alpha_1}}\right)<0,$

when $E_1^+$ exists since $\alpha_1\leq(\sqrt{\beta_1}-1)^2$. Thus, $E_1^+$ will be always stable if $\beta_1<\left(1+\frac{\beta_3}{d\epsilon}\right)$ or it can be stable if $\left(1+\frac{\beta_3}{d\epsilon}\right)<\beta_1<\left(1+\frac{\beta_3}{d\epsilon}\right)^2$. Otherwise $E_1^+$ is a saddle point.

Proposition 3. Assume $\beta_2\geq(\sqrt{\gamma}+\sqrt{\alpha_2})^2$ then,

(a) the model (5) can have at most two axial equilibria given by $E_2^+(0, y_2^+, 0)$ and $E_2^-(0, y_2^-, 0)$ respectively.

(b) $E_2^-$ is always a saddle point whenever exists. $E^+_2$ is locally asymptotically stable if $y^+_2<\frac{\beta_3}{d \mu}$ and a saddle point if $y^+_2>\frac{\beta_3}{d \mu}$. If we consider $\alpha_2$ to be the varying parameter then $E_2^+$ will be always stable if $\beta_2<\left(\gamma+\frac{\beta_3}{d\mu}\right)$ or there exists a range of values of $\alpha_2$ such that $E_2^+$ is stable if $\left(\gamma+\frac{\beta_3}{d\mu}\right)<\beta_2<\frac{1}{\gamma}\left(\gamma+\frac{\beta_3}{d\mu}\right)^2$.

Proof.  Similar as Proposition 2.

Proposition 4. The model (5) has a boundary equilibrium point $E_3(x_3, 0, z_3)$ where the second prey population is absent if $d\beta_1\beta_3\epsilon>(\beta_3+d\epsilon)(\beta_3+d\alpha_1\epsilon)$. $E_3$ is locally asymptotically stable if $\alpha_1\beta_1<\Big(\frac{\beta_3}{d \epsilon}+\alpha_1\Big)^2$. Otherwise $E_3$ is a saddle point.

Proof.  Substituting $y = 0$ in (5) we have,

$x_3 = \frac{\beta_3}{d\epsilon}, \, \, y_3 = 0, \, \, z_3 = \frac{d\beta_1\beta_3\epsilon-(\beta_3+d\epsilon)(\beta_3+d\alpha_1\epsilon)}{(\beta_3+d\alpha_1\epsilon)d\epsilon}$

(7)

It is easily seen that the in this case a unique equilibrium, denoted by $E_3(x_3, 0, z_3)$, is feasible according to the following conditions.

$\bullet$  Case1. $d\beta_1\beta_3\epsilon<(\beta_3+d\epsilon)(\beta_3+d\alpha_1\epsilon)$.In this case $E_3$ does not exist.

$\bullet$  Case2. $d\beta_1\beta_3\epsilon>(\beta_3+d\epsilon)(\beta_3+d\alpha_1\epsilon)$.In this case $E_3$ is feasible.

The Jacobian at $E_3$ is given by,

$J(E_3) = \begin{pmatrix} x_3\left[\frac{\alpha_1\beta_1}{(x_3+\alpha_1)^2}-1\right]&-\alpha x_3&-\epsilon x_3 \\ 0&-\gamma-\beta x_3-z_3&0 \\ d\epsilon z_3&d\mu z_3&0 \end{pmatrix}.$

It is easy to see that $E_3$ is locally asymptotically stable if $\alpha_1\beta_1<(x_3+\alpha_1)^2$ as all the three eigenvalues of the matrix $J(E_3)$ are negative in this case. Otherwise $E_3$ is a saddle point.

Proposition 5.   The model (5) has a boundary equilibrium point $E_4(0, y_4, z_3)$ where the first prey population is absent if $d\beta_2\beta_3\mu>(\beta_3+d\mu\gamma)(\beta_3+d\alpha_2\mu)$. $E_4$ is locally asymptotically stable if $\alpha_2\beta_2<\Big(\frac{\beta_3}{d \mu}+\alpha_2\Big)^2$. Otherwise $E_4$ is a saddle point.

Proof.  Similar as Proposition 4.

Proposition 6. If $\beta_1>(1+\sqrt{\alpha_1})^2$ & $\beta_2>(\sqrt{\gamma}+\sqrt{\alpha_2})^2$, the model (5) can exhibit at most four boundary equilibrium point of the form $E_5(x_5, y_5, 0)$ where the predator population is absent of which at most one equilibrium point can be locally asymptotically stable. Otherwise all are unstable.

Proof.  Let us consider the case of the existence of the predator free equilibrium point. Thus if we plug $z = 0$ in (5) we have system of two quadratic equations given by,

$\frac{\beta_1 x}{x+\alpha_1} = 1+x+\alpha y,$

(8)

$\frac{\beta_2 y}{y+\alpha_2} = \gamma+\beta x+y.$

(9)

For each solution with positive components of the above equations the system (5) will possess an equilibrium point given by $E_5(x_5, y_5, 0)$. To analyze the existence of positive solutions of the above system of equations we consider the following cases.

If we assume $\beta_1\leq(1+\sqrt{\alpha_1})^2$ or $\beta_2\leq(\sqrt{\gamma}+\sqrt{\alpha_2})^2$, then from the result of $E_1$ & $E_2$, we know that,

$\beta_1 x-(1+x)(x+\alpha_1)<0, \,\,\,\,\,\forall x\geq0, \\ \textrm{or, }\, \, \beta_2 y-(\gamma+y)(y+\alpha_2)<0, \,\,\,\,\,\forall y\geq0.$

Hence the $E_5$ does not exist in this case.

Now if both $\beta_1>(1+\sqrt{\alpha_1})^2$ & $\beta_2>(\sqrt{\gamma}+\sqrt{\alpha_2})^2$ hold, (8) and (9) can be written as

$y = -\frac{(x-x^+_1)(x-x^-_1)}{\alpha(x+\alpha_1)} \, \, \, \, \& \, \, \, \, x = -\frac{(y-y^+_2)(y-y^-_2)}{\beta(y+\alpha_2)},$

respectively. (8) has a local maxima at $x = \sqrt{\alpha_1}(\sqrt{\beta_1}-\sqrt{\alpha_1})>0$ and (9) has a local extrema at $y = \sqrt{\alpha_2}(\sqrt{\beta_2}-\sqrt{\alpha_2})>0$. Now if we solve for the value of $y$ from the above equations we will have a fourth degree polynomial in $x$ for which at most four feasible solutions are possible. It is difficult to find out the exact parametric restrictions for which that polynomial will have four positive roots. But graphically we can assert the possible feasibility and positions of those roots. Hence, we may have at most four possible $E_5$ as shown in Fig. 1.

Now the Jacobian at $E_5$ is given by,

$J(E_5) = \begin{pmatrix} x_5\left[\frac{\alpha_1\beta_1}{(x_5+\alpha_1)^2}-1\right]&-\alpha x_5&-\epsilon x_5 \\ -\beta y_5&y_5\left[\frac{\alpha_2\beta_2}{(y_5+\alpha_2)^2}-1\right]&-y_5 \\ 0&0&-\beta_3+d\epsilon x_5+d\mu y_5 \end{pmatrix}\\ =\left( {\frac{{\begin{array}{*{20}{c}} {}&{}&{\left| {} \right.}&{ - {\rm{e}}x5}\\ {}&{{J_{{x_5}{y_5}}}}&{\left| {} \right.}&{\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; - {y_5}} \end{array}}}{{\begin{array}{*{20}{c}} 0&0&{\left| {} \right.}&{ - {\beta _3} + d{\rm{e}}{x_5} + d\mu {y_5}} \end{array}}}} \right)$

Thus $E_5$ is locally asymptotically stable iff,

$-\beta_3+d\epsilon x_5+d\mu y_5<0, \, \, \displaystyle \textrm{Tr}(J_{x_5 y_5})<0\, \, \&\, \, \displaystyle \textrm{Det}(J_{x_5 y_5})>0$

Also $E_5$ is unstable if one of the following conditions satisfy.

$-\beta_3+d\epsilon x_5+d\mu y_5>0,$

(10)

$\left(\frac{\alpha_1\beta_1}{(x_5+\alpha_1)^2}-1\right)\left(\frac{\alpha_2\beta_2}{(y_5+\alpha_2)^2}-1\right)<0,$

(11)

$\frac{\alpha_1\beta_1}{(x_5+\alpha_1)^2}>1\, \, \&\, \, \frac{\alpha_2\beta_2}{(y_5+\alpha_2)^2}>1.$

(12)

Although we are unable to describe the complete analytical expression for which $E_5$'s are asymptotically stable or unstable, it can be easily verified that all the equilibria with red dot in Fig. 1 are unstable since either (11) or (12) is satisfied. Only the equilibrium point pictured as green dot in Fig. 1 may be stable. Hence the proposition is shown.

Proposition 7.  Let $\beta_1>(1+\sqrt{\alpha_1})^2$ & $\beta_2>(\sqrt{\gamma}+\sqrt{\alpha_2})^2$ hold.

(a) The model (5) can have at most three interior equilibira denoted by $E_{1*}(x_1, y_1, z_1)$, $E_{2*}(x_2, y_2, z_2)$ and $E_{3*}(x_3, y_3, z_3)$ such that, $x_1>x_2>x_3$.

(b) $E_{1*}$ and $E_{3*}$ always remain unstable. $E_{2*}$ can be asymptotically stable or unstable depending on suitable parametric restrictions.

Proof.  Lastly we concentrate on the interior equilibrium point i.e. the equilibrium point of (5) for which all the components are strictly positive. The nullclines representing the interior equilibrium point are given by,

$\frac{\beta_1 x}{x+\alpha_1}-1-x-\alpha y = \epsilon z,$

(13)

$\frac{\beta_2 y}{y+\alpha_2}-\gamma-\beta x-y = z,$

(14)

$d \epsilon x+d \mu y = \beta_3.$

(15)

Similar to Proposition 6. if $\beta_1\leq(1+\sqrt{\alpha_1})^2$ or $\beta_2\leq(\sqrt{\gamma}+\sqrt{\alpha_2})^2$, it is clear that,

$\frac{\beta_1 x}{x+\alpha_1}-(1+x)-\alpha y<0, \,\,\,\,\,\,\forall x, y\geq0, \\ \textrm{or, }\, \, \frac{\beta_2 y}{y+\alpha_2}-(\gamma+y)-\beta x<0, \,\,\,\,\,\,\forall x, y\geq0.$

Hence no interior equilibrium point is feasible in this case.

Now if $\beta_1>(1+\sqrt{\alpha_1})^2$ and $\beta_2>(\sqrt{\gamma}+\sqrt{\alpha_2})^2$, then from the above three equations if we solve for $x$ by substituting the values of $y \, \&\, z$ we get,

$G(x)\equiv A_1x^3+A_2x^2+A_3x+A_4 = 0$

(16)

where,

$A_1 = -\alpha d^2 \epsilon^2 + d^2 \epsilon^3 + d^2 \epsilon \mu-\beta d^2 \epsilon^2 \mu, \\ A_2 = 2 \alpha \beta_3 d \epsilon - 2 \beta_3 d \epsilon^2-\alpha \alpha_1 d^2 \epsilon^2 +\alpha_1 d^2 \epsilon^3-\beta_3 d \mu + \beta \beta_3 d \epsilon \mu + d^2 \epsilon\mu +\alpha_1 d^2 \epsilon\mu \\ \, \, \, \, \, \, + \alpha \alpha_2 d^2 \epsilon\mu-\beta_1d^2\epsilon\mu-\alpha_2 d^2\epsilon^2\mu-\alpha_1\beta d^2 \epsilon^2 \mu + \beta_2 d^2 \epsilon^2 \mu-d^2 \epsilon^2 \gamma \mu - \alpha_2 d^2 \mu^2 \\ \, \, \, \, \, \, + \alpha_2 \beta d^2 \epsilon \mu^2, \\ A_3 = -\alpha \beta_3^2 + \beta_3^2 \epsilon + 2 \alpha \alpha_1 \beta_3 d \epsilon-2 \alpha_1\beta_3 d \epsilon^2-\beta_3 d \mu - \alpha_1 \beta_3 d \mu - \alpha \alpha_2 \beta_3 d \mu \\ \, \, \, \, \, \, + \beta_1 \beta_3 d \mu + \alpha_2 \beta_3 d \epsilon \mu + \alpha_1 \beta \beta_3 d \epsilon \mu -\beta_2 \beta_3 d \epsilon \mu + \alpha_1 d^2 \epsilon \mu + \alpha \alpha_1 \alpha_2 d^2 \epsilon \mu$

$\, \, \, \, \, - \alpha_1 \alpha_2 d^2 \epsilon^2 \mu + \alpha_1 \beta_2 d^2 \epsilon^2 \mu + \beta_3 d \epsilon \gamma \mu -\alpha_1 d^2 \epsilon^2 \gamma \mu - \alpha_2 d^2 \mu^2\\ \, \, \, \, \, - \alpha_1 \alpha_2 d^2 \mu^2 + \alpha_2 \beta_1 d^2 \mu^2 + \alpha_1 \alpha_2 \beta d^2 \epsilon \mu^2 + \alpha_2 d^2 \epsilon \gamma \mu^2, \\ A_4 = -\alpha \alpha_1 \beta_3^2 + \alpha_1 \beta_3^2 \epsilon - \alpha_1 \beta_3 d \mu - \alpha \alpha_1 \alpha_2 \beta_3 d \mu + \alpha_1 \alpha_2 \beta_3 d \epsilon \mu - \alpha_1 \beta_2 \beta_3 d \epsilon \mu \\ \, \, \, \, \, + \alpha_1 \beta_3 d \epsilon \gamma \mu -\alpha_1 \alpha_2 d^2 \mu^2 + \alpha_1 \alpha_2 d^2 \epsilon \gamma \mu^2.$

The system (5) will have interior equilibrium points only if (16) has positive root $x_*$ such that $y_* = \frac{\beta_3-d\epsilon x_*}{d\mu}>0$ and $z_* = -\gamma-\beta x_*-y_*+\frac{\beta_2 y_*}{\alpha_2+y_*}>0$.

The Jacobian at $E_*$ is given by,

$J(E_*) = \begin{pmatrix} x_*\left[\frac{\alpha_1\beta_1}{(x_*+\alpha_1)^2}-1\right]&-\alpha x_*&-\epsilon x_* \\ -\beta y_*&y_*\left[\frac{\alpha_2\beta_2}{(y_*+\alpha_2)^2}-1\right]&-y_* \\ d\epsilon z_*&d\mu z_*&0 \end{pmatrix}.$

The corresponding characteristic equation is,

$\lambda^3+B_1 \lambda^2+B_2\lambda+B_3 = 0$

where,

$B_1 = -\left(x_*\left[\frac{\alpha_1\beta_1}{(x_*+\alpha_1)^2}-1\right]+y_*\left[\frac{\alpha_2\beta_2}{(y_*+\alpha_2)^2}-1\right]\right), \\ B_2 = \left(x_*y_*\left[\frac{\alpha_1\beta_1}{(x_*+\alpha_1)^2}-1\right]\left[\frac{\alpha_2\beta_2}{(y_*+\alpha_2)^2}-1\right]+d\mu y_*z_*-\alpha\beta x_*y_*\right), \\ B_3 = -\displaystyle \textrm{Det}(J(E_*)).$

The interior equilibrium point will be locally asymptotically stable if $B_i>0, \, \, i = 1, 2, 3$, and $B_1B_2>B_3$.

It is evident from the equation (16) that the system (5) can have at most three interior equilibrium points. Although it is quite difficult to find out the analytical conditions for their existence and stability, numerically it can be easily verified that the system can possess three feasible interior equilibira. Also extensive numerical results confirm that out of the three interior equilibrium points only one can be stable while the other will remain unstable whenever they exist. To validate our claim we give some numerical results as follows.

Let us consider the parameter set $\alpha = 1, \alpha_1 = .001, \alpha_2 = .00001, \gamma = 1, \beta = 1.5, \beta_3 = 1, \epsilon = 4, d = .5, \mu = 1.$

$\bullet$  For $\beta_1 = 2.5, \beta_2 = 2.6$ the system (5) has only one feasible interior equilibrium points given by, $E_*(x_*, y_*, z_*) = (0.499992, 0.0000332445, 0.248746).$ In this case $E_*$ is an unstable point with two stable and one unstable manifolds.

$\bullet$  For $\beta_1 = 3.5, \beta_2 = 2.6$ the system possess two interior equilibria $E_{1*}$ and $E_{2*}$ whose components are given by $E_{1*}(0.499984, 0.0000639227, 0.498241)$ and $E_{2*}(0.298348, 0.806608, 0.345838)$. Of the two $E_{1*}$ is an unstable equilibrium point with two stable and one unstable manifolds and $E_{2*}$ is a locally asymptotically stable equilibrium point.

$\bullet$  If $\beta_1 = 5, \beta_2 = 3$ the system (5) exhibits three feasible interior equilibrium points namely, $E_{1*}(0.499983, 0.0000694776, 0.872492)$, $E_{2*}(0.283221, 0.867116, 0.708018)$ and $E_{3*}(0.00151324, 1.99395, 0.00376804)$. Here $E_{1*}$ and $E_{3*}$ are unstable equilibria both having two stable and one unstable manifolds. $E_{2*}$ is locally asymptotically stable.

---

4. Local bifurcations:

The present section mainly reflects on how the equilibrium points appear or disappear from one another and how the stability of the equilibria changes through different types of local or global bifurcations.

In proposition 8 we present different conditions how the two branches of different equilibria appear or disappear through several saddle node bifurcations.

Proposition 8.  (a) The system (5) undergoes a saddle-node bifurcation at $\alpha_1 = (\sqrt{\beta_1}-1)^2$, when $E_1^+$ and $E_1^-$ coincide.

(b) The system (5) undergoes a saddle-node bifurcation at $\alpha_2 = (\sqrt{\beta_2}-\sqrt{\gamma})^2$, when $E_2^+$ and $E_2^-$ coincide.

(c) Suppose that at $\alpha_1 = \alpha_1^{5*}$, $\displaystyle \textrm{Det}(J_{x_5 y_5})|_{\alpha_1^{5*}} = 0$. Then the system (5) undergoes a saddle-node bifurcation at $\alpha_1 = \alpha_1^{5*}$ when two $E_5$ coincide.

(d) If for $\alpha_1 = \alpha_1^{SN}$, $G(x) = 0$ has a double root then the system (5) exhibits another saddle node bifurcating where two interior equilibrium points coincide.

Proof.  (a) Let us assume $\alpha_1 = \alpha_1^* = (\sqrt{\beta_1}-1)^2$. It is easy to see that $E_1^+$ and $E_1^-$ coincide at $\alpha_1 = \alpha_1^*$. Let $v_1 = \begin{pmatrix} 1\\ 0\\ 0 \end{pmatrix}$ and $w_1 = \begin{pmatrix} 1\\ -\frac{\alpha(\sqrt{\beta_1}-1)}{\gamma+\beta(\sqrt{\beta_1}-1)}\\ -\frac{\epsilon(\sqrt{\beta_1}-1)}{-\beta_3+d\epsilon(\sqrt{\beta_1}-1)} \end{pmatrix}$ be the eigenvectors corresponding to the zero eigenvalue of the matrices $J(E_1)$ and $J(E_1)^T$ respectively at $\alpha_1 = (\sqrt{\beta_1}-1)^2$.

Let us rewrite the system (5) as $\frac{dX}{dt} = F$. Then we have the following

$w_1^T F_{\alpha_1}|_{\alpha^*_1} = -1\, \, \neq 0, \\ w_1^T D^2F|_{\alpha^*_1}(v_1, v_1) = -\frac{2}{\sqrt{\beta_1}}\, \, \neq 0 .$

Thus the system undergoes saddle-node bifurcation at $\alpha_1 = \alpha_1^*$.

(b) The proof is similar as given in (a).

(c) The slopes of the curves (8) and (9) at any point $(x, y)$ are respectively given by $\frac{1}{\alpha}\left[\frac{\alpha_1\beta_1}{(x+\alpha_1)^2}-1\right]$ and $\frac{\beta}{\left[\frac{\alpha_2\beta_2}{(y+\alpha_2)^2}-1\right]}$. Now $\displaystyle \textrm{Det}(J_{x_5 y_5})|_{\alpha_1^{5*}} = 0$ implies that the curves (8) and (9) touches each other and consequently two $E_5$ coinside. Now if we proceed as above then it is easy to prove that the system undergoes a saddle node bifurcation at $\alpha = \alpha_1^{5*}$.

(d) Similar to (a) & (c).

In the following we discuss how the equilibira $E_3$ and $E_4$ exchange their stability with that of $E_1$ and $E_2$ respectively through transcritical bifurcations.

Proposition 9.(a) The system (5) undergoes Transcritical bifurcation at $\beta_3 = \beta^*_3 \equiv \frac{d\epsilon\left[-(1+\alpha_1-\beta_1)\pm \sqrt{(1+\alpha_1-\beta_1)^2-4\alpha_1}\right]}{2}$ and $E_3$ exchanges stability with $E_1^+$ if $\beta_1<\left(1+\frac{\beta_3}{d*\epsilon}\right)^2$ and with $E_1^-$ if $\beta_1>\left(1+\frac{\beta_3}{d*\epsilon}\right)^2$ respectively.

(b) The system (5) undergoes another Transcritical bifurcation at $\beta_3 = \beta^*_3 \equiv \frac{d\mu\left[-(\gamma+\alpha_2-\beta_2)\pm \sqrt{(\gamma+\alpha_2-\beta_2)^2-4\alpha_2\gamma}\right]}{2}$ and $E_4$ exchanges stability with $E_2^+$ if $\beta_2<\frac{1}{\gamma}\left(\gamma+\frac{\beta_3}{d*\mu}\right)^2$ and with $E_2^-$ if $\beta_2>\frac{1}{\gamma}\left(\gamma+\frac{\beta_3}{d*\mu}\right)^2$ respectively.

Proof.   (a) Let $v_3 = \begin{pmatrix} \frac{\epsilon}{\frac{\alpha_1\beta_1}{(x_3+\alpha_1)^2}-1}\\ 0\\ 1 \end{pmatrix}$ and $w_3 = \begin{pmatrix} 0\\ 0\\ 1 \end{pmatrix}$ be the eigenvectors corresponding to the zero eigenvalue of the matrices $J(E_3)$ and $J(E_3)^T$ respectively.

Then we have the following,

$w_3^T F_{\beta_3}|_{\beta^*_3} = 0, \\ w_3^T DF_{\beta_3}v_3|_{\beta^*_3} = -1, \\ w_3^T D^2F(v_3, v_3)|_{\beta^*_3} = \frac{2d\epsilon^2}{\frac{\alpha_1\beta_1}{(x_3+\alpha_1)^2}-1}.$

(b) Similarly as above.

In the next proposition we see how under different parametric restrictions the system exhibits several Hopf-bifurcations of different equilibria. Specially we focus on how the interior equilibrium point looses its stability through Hopf-bifurcation. subsequently in proposition 11 we discuss the possible existence of Bogdanov-Takens and generalized Hopf bifurcations of the interior equilibrium point.

Proposition 10.   (a) The equilibrium point $E_3$ undergoes an Hopf-bifurcation at $\alpha_1 = \alpha_1^{H_3}$ where $\alpha_1^{H_3}$ is a root of the equation $\beta_1-\frac{(\beta_3+d\epsilon\alpha_1)^2}{d^2\epsilon^2\alpha_1} = 0$.

(b) The equilibrium point $E_4$ undergoes an Hopf-bifurcation at $\alpha_2 = \alpha_2^{H_4}$ where $\alpha_2^{H_4}$ is a root of the equation $\beta_2-\frac{(\beta_3+d\mu\alpha_2)^2}{d^2\mu^2\alpha_2} = 0$.

(c) The interior equilibrium point $E_{2*}$ looses its stability at $\alpha_1 = \alpha_1^{H_*}$ such that $\displaystyle J_{E_{2*}}$ has two purely imaginary eigenvalues.

Proof.   (a) At the threshold value $\beta_1 = \beta_1^{H_3} = \frac{(\beta_3+d\epsilon\alpha_1)^2}{d^2\epsilon^2\alpha_1}$ the eigenvalues of the Jacobian at $E_3$ for the system (5) are given by, $\lambda_1 = -\gamma-\beta x_3-z_3, \, \, \lambda_2 = +i\sqrt{\epsilon\beta_3z_3}, \, \, \lambda_3 = -i\sqrt{\epsilon\beta_3z_3}$. Also $\frac{d}{d\beta_1}\displaystyle (Tr(J_{E_3}))|_{\beta_1 = \beta_1^{H_3}} = \frac{d\alpha_1\mu}{(\beta_3+d\alpha_1\mu)^2}\neq0$.

Hence, $E_3$ undergoes a Hopf-bifurcation at $\beta_1 = \beta_1^{H_3}$.At the threshold value $\beta_1 = \beta_1^{H_3} = \frac{(\beta_3+d\epsilon\alpha_1)^2}{d^2\epsilon^2\alpha_1}$ the eigenvalues of the Jacobian at $E_3$ for the system (5) are given by, $\lambda_1 = -\gamma-\beta x_3-z_3, \, \, \lambda_2 = +i\sqrt{\epsilon\beta_3z_3}, \, \, \lambda_3 = -i\sqrt{\epsilon\beta_3z_3}$. Also $\frac{d}{d\beta_1}\displaystyle (Tr(J_{E_3}))|_{\beta_1 = \beta_1^{H_3}} = \frac{d\alpha_1\mu}{(\beta_3+d\alpha_1\mu)^2}\neq0$.

Hence, $E_3$ undergoes a Hopf-bifurcation at $\beta_1 = \beta_1^{H_3}$.

(b) By similar arguments as above it can be shown that $E_4$ undergoes a Hopf-bifurcation at $\beta_2 = \frac{(\beta_3+d\mu\alpha_2)^2}{d^2\mu^2\alpha_2}$.

(c) It is quite difficult to give an analytical proof for the Hopf-bifurcation of the interior equilibrium point $E_{2*}$ but extensive numerical simulations show that the under certain conditions $E_{2*}$ undergoes Hopf bifurcations which can be either supercritical or subcritical depending upon the choice of different parametric combinations.

Proposition 11.   (a) If there exists a set of parameters such that the Jacobian $\displaystyle J_{E_{2*}}$ has zero eigenvalue of multiplicity two, the system undergoes a Bogdanov-Tackens bifurcation at the point $E_{2*}$.

(b) If there exists a set of parameters such that the Lyapunov coefficient of the Hopf-bifurcating limit cycle around $E_{2*}$ is zero then the system also exhibits a Bautin or Generalized Hopf bifurcation.

Proof.  Analytical conditions for these two bifurcations are difficult to produce. But numerically it can be easily verified that the system exhibits both these bifurcations. Here we present two such parameter sets in the following.

$\bullet$  The system undergoes BT bifurcation at $\alpha = 1, \alpha_1 = .001, \alpha_2 = .00001, \gamma = 1, \beta = 1.5, \beta_1 = 4.290306, \beta_2 = 3.2376946, \beta_3 = 1, \epsilon = 4, d = .5, \mu = 1.$

$\bullet$  $E_{2*}$ undergoes Bautin or generalized Hopf bifurcation at $\alpha = 1, \alpha_1 = .001, \alpha_2 = .00001, \gamma = 1, \beta = 1.5, \beta_1 = 2.977584, \beta_2 = 2.7448137, \beta_3 = 1, \epsilon = 4, d = .5, \mu = 1.$

---

5. Local and global bifurcations: Numerical simulation results

In this section we mainly focus on how the introduction of Allee effect influences the dynamics of the underlying system. We construct two dimensional bifurcation diagrams taking $\alpha_1$ and $\alpha_2$ as the bifurcation parameters. We will also try to make comparison of our results with the results for the model considered in [17] pp. 62-72. Here we present two bifurcation diagrams in the $\alpha_1\alpha_2$ parametric plane shown in Fig. 2 and Fig. 5.

In the first schematic bifurcation diagram, presented in Fig. 2 we have considered the parameter set for which the system without Allee effect described in [17] pp. 62-72 possesses locally stable coexisting equilibrium point. Now we consider the results for the model(5). The trivial equilibrium point $E_0(0, 0, 0)$ is always locally asymptotically stable, for the system (5), irrespective of any parametric restrictions and hence we are not going to mention its stability at any parametric domain of the bifurcation diagrams. It is worthy to mention that for any choice of parameter values there exists a nonempty basin of attraction for $E_0$. In first bifurcation diagram presented at Fig. 2 the two vertical lines given in violet and black colours, whose expressions can be obtained from Proposition 9 and Proposition 8 respectively, represent the transcritical bifurcation curve for $E_3$ and saddle-node bifurcation curve for the equilibrium point $E_1$. On the left of the violet line there is one more yellow coloured vertical line which is the Hopf-bifurcation curve for $E_3$ as described in the Proposition 10 On the left of this line i.e. in the regions $R_1$ to $R_7$, $E_3$ is locally asymptotically stable and in the regions in between the yellow and violet curves $E_3$ is unstable and it looses stability through a Hopf-bifurcation on the yellow line. Thus as we increase the value of $\alpha_1$, keeping $\alpha_2$ fixed, the parameters move through the domains $R_1\, \rightarrow\, R_{15}\, \rightarrow\, R_{16}\, \rightarrow\, R_{17}$. $E_3$ is stable in $R_1$ and looses its stability through Hopf-bifurcation as parameter $\alpha_1$ moves from $R_1$ to $R_{15}$ and then gradually disappears through the appearance of $E_1$ and later the two branches of $E_1$ disappear through a saddle-node bifurcation. The same phenomena for $E_1$ and $E_3$ are observed as $\alpha_1$ moves through $R_2/R_3\, \rightarrow\, R_{14}\, \rightarrow\, R_{13}\, \rightarrow\, R_{12}$. The dark cyan coloured horizontal line signifies the saddle-node bifurcation curve for $E_2$, below the line we find two equilibira $E_2^+$ and $E_2^-$ among which $E_2^+$, is locally asymptotically stable but $E_2^-$ is a saddle point and as $\alpha_2$ crosses the line from below they disappear through saddle node bifurcation.

The light cyan and the magenta curves are the saddle-node bifurcation curve for $E_5$ and the transcritical bifurcation curve for $E_*$ respectively. The two curves touch each other at a point given by the blue dot. Below the cyan curve $E_5^1$ and $E_5^2$ exist but none of which is stable as mentioned in section 3. Above this curve both the equilibria disappear through the saddle-node bifurcation, the threshold for this saddle-node bifurcation is discussed in Proposition 8. The magenta curve is the transcritical bifurcation curve of the equilibrium point $E_*$ whose equation can not be obtained explicitly. As $\alpha_1$ moves from the left of this curve to the region $R_7$ ($R_6\, \rightarrow\, R_7$), first a transcritical bifurcation takes place through which $E_{1*}$ disappears, and then one more transcritical bifurcation occurs as $\alpha_1$ crosses from $R_8\, \rightarrow\, R_9$ where $E_{2*}$ disappears, there is an exchange of stability with $E_5^2$. Now on the left of the magenta curve and below the light cyan curve there is the dotted blue curve which represents the saddle-node bifurcation curve for the interior equilibrium points. This curve touches the magenta curve at the red dot. The magenta and the blue dotted curve together divide the parametric plane into different domains signifying the existence of the interior equilibria. The domains $R_7, R_8$ bounded above by the magenta curve where one interior equilibrium point is feasible which is locally unstable; no interior point exists in the domains lying above the magenta and the blue dotted curve. The domains $R_4, R_5, R_6$ bounded by the magenta and the blue dotted curve contain two interior equilibria of which $E_{1*}$ is always locally unstable and stability of $E_{2*}$ depends on different parametric restrictions. Now within the domain bounded by the magenta and the blue dotted curve there are three curves including the blue dotted curve coinciding at a point marked as black dot. This black dot represents the Bogdanov-Taken's(BT) bifurcation point, which is the point of intersection of the saddle-node bifurcation curve and Hopf-bifurcation curve for the interior equilibrium point $E_2*$. The green curve arising from the BT point represents the Hopf-bifurcation curve of $E_2*$. On the left of this curve i.e in the region $R_4$, $E_2*$ is locally asymptotically stable and in the region bounded by the green and red curves $E_2*$ looses its stability through a supercritical Hopf-bifurcation and thus a locally stable limit cycle appears around $E_2*$. We have observed that as $\alpha_1$ increases beyond the Hopf-bifurcation threshold, two periodic orbit exists which appears through a period doubling at $\alpha_1 = 0.0118$ (see Fig. 3). Interstinly if we further increase the value of $\alpha_1$ the system exhibits peak adding bifurcation, the first peak appears at $\alpha_1 = 0.0121$ which is reflected at the bifurcation diagram. Appearance of peaks are shown in Fig. 4 and ultimately we observe chaotic dynamics and we find chaotic attractor around $E_2*$. The red curve in Fig. 2 represents the first period doubling bifurcation curve through which the limit cycle first changes its period from 1 to 2, whenever $\alpha_1$ crosses it. Clearly the route to chaos is not period doubling rather it is a combination of period doubling and peak adding bifurcations.

An illustration of peak adding bifurcation is shown in Fig. 4. In this figure we have plotted the time series for first prey population ($x_{11}$) against time, after deleting significant amount of initial transients, as $\alpha_1$ increases from 0.0119 to 0.0123 and the other parameters are fixed at $\alpha = 1$, $\alpha_2 = 0.01$, $\beta = 1.5$, $\beta_1 = 2$, $\beta_2 = 2$, $\beta_3 = 1$, $\gamma = 1$, $d = 0.5$, $\mu = 1$, $\epsilon = 5$. In peak adding bifurcation supplementary local maxima and minima emerges successively but the period of the limit cycle does not change due to this bifurcation. Supplementary local maxima and minima appear through the appearance of point of inflexions in the time series plot. Another characteristic feature is the difference in the heights of supplementary maxima and minima with the change of parameter value. Detailed description of peak-adding bifurcation and its application in the context of single species population model are available at [9,13].

We have presented another bifurcation diagram at Fig. 3 for a better understanding of the dynamics as $\alpha_1$ moves through the domains $R_4\, \rightarrow\, R_5\, \rightarrow\, R_6$. In this figure we have presented $x-$ components of various equilibrium points against a range of values for $\alpha_1$ varies. Here blue, green, red, magenta, black lines represent the $x-$ components of $E_0,$ $E_3,$ $E_5,$ $E_{1*},$ $E_{2*}$ respectively where the dotted lines represent the components of unstable equilibia and the solid lines represent components of the stable equilibia. In case of periodic and chaotic dynamics, local maxima and minima of concerned $x-$ component is plotted after deleting the initial transients.Stability of various equilibia for parameter values lying in different domains of Fig. 2 are summarized in Table. 2.

The second schematic diagram is presented in Fig. 5. In this case the parameters are the same as discussed in [17] pp. 62-72 and the system is in the chaotic regime in the absence of the Allee effects. This bifurcation diagram is quite similar to that of Fig. 2. The only visible difference is due to the change in positions of the yellow vertical line representing the Hopf-bifurcation curve for the equilibrium point $E_3$ or the saddle-node bifurcation curve of $E_1$(black vertical line). As a result three new qualitatively different regions $R_{4A}$, $R_{5A}$ and $R_{6A}$ have come up and one region $R_{11}$ disappears. The Figure-5 is prepared to have more insight to understand how the dynamics is changing if $\alpha_2$ moves through $R_6\, \rightarrow\, R_5\, \rightarrow\, R_4$. In this diagram the blue, green, red, magenta, black lines represent the $x-$ components of $E_0,$ $E_3,$ $E_5,$ $E_{1*},$ $E_{1*}$. Stability of various equilibria for parameter values lying in different domains of Fig. 5 are summarized in Table. 3.

---

6. Conclusion

Prey-predator models with Allee effect in prey growth recently have received significant attention from the researchers [1,2,5,6,12,14,15]. Several ecological species are identified which exhibit Allee effects due to various reasons [7,8,16]. Models with one prey and their specialist predator with various types of functional responses exhibit comparatively rich dynamics compared to the corresponding models without Allee effect. Most common observation for these investigations are the possibility of system's collapse due to the extinction of both the species depending upon their initial population densities and also due to some global bifurcations when grazing pressure on prey species is significantly high. That both the prey and predator species become extinct, depending upon the initial population densities, is a common feature for the models with strong Allee effect. Here we have made an attempt to understand the influence of Allee effect on a three dimensional prey-predator model consisting with two prey and one predator. In some sense the model can be considered as a prey-predator model with generalist predator also as the predator can survive on any one of two prey populations. We have introduced Allee effect in the growth equations for both the prey species and the Allee effects are known to be additive in nature [1,2,21].

Firstly we admit that the inclusion of Allee effects in the growth equations of both the prey species makes the mathematical analysis quite difficult and in most of the cases we are unable to find explicit conditions for stability of equilibrium point(s) and thresholds for various local bifurcations. However, with the help of numerical simulations we have explored the rich dynamics exhibited by the model by considering the Allee effect parameters as bifurcation parameters. In case of the three dimensional model we have considered here, the trivial equilibrium point is always stable as the basin of attraction of the extinction steady-state is a non-empty set under any choice of parameter values (this is clear from Table-1 and Table-2). All possible local and global bifurcation scenarios are presented in two schematic bifurcation diagrams, where we have used schematic diagrams as some of the bifurcation curves are very close to each other when we plot them against actual parameter values. Another important observation is the suppression of chaos due to the Allee effects in prey growths as we have observed chaotic dynamics for a short range of values for the Allee effect parameter. However the appearance and disappearance of chaos is not only due to period-doubling and reverse period-doubling bifurcations rather we have observed the appearance of peak adding bifurcation also. Hence we can say that Allee effects in prey growths can suppress the chaotic dynamics and the route to chaos is different from the model without Allee effect.

Negative growth rate of prey population at their low population density results in the extinction of one or more species depending upon the strengths of various interactions as well as the initial population densities. However with the increased strength of Allee effect on any one or both the prey population always drives the system towards total extinction. Our claim is based upon the stability of extinction steady-state in the regions $R_{15}$, $R_{16}$, $R_{17}$ in Fig. 2 and in the regions $R_{14}$, $R_{15}$, $R_{16}$ in Fig. 4. All the populations coexist at their steady-state or exhibit oscillatory aperiodic coexistence when the strengths of Allee effects are not very high and of course depending upon the initial population densities. Another interesting feature is the appearance of tri-stability for a range of parameter values. Fig. 3 shows that $E_0$, $E_3$ and $E_{2*}$ are stable for $\alpha_1<\alpha_{1_H}$, we see the stability of $E_0$ and $E_3$ and oscillatory or aperiodic coexistence of three species for $\alpha_{1_H}<\alpha_1<\overline{\alpha}_1$. Here $\alpha_1 = \alpha_{1_H}$ is the Hopf-bifurcation threshold and chaotic dynamics disappears through crisis at $\overline{\alpha}_1$. In order to visualize the existence of chaotic regime we have plotted the bifurcation diagram for a short range of values of $\alpha_1$ but for large $\alpha_1$ we find extinction of one or more species. Similar argument holds for Fig. 5. In summary, the introduction of Allee effects in both the prey population induces rich dynamics due to appearance of various equilibrium points and their change in stability behaviors due to number of local and global bifurcations. Survival of three species is solely dependent upon the strengths of inter-and intra-specific interactions as well as the initial population densities. It is important to mention here that the conclusions are based upon a relatively simple model, as the consumption of prey by the predator is assumed to follow the law of mass action. Our future goal will be to examine dynamics of similar or other type of models with a saturating functional response as well as predator dependent functional responses.

---

Acknowledgments

We are grateful to the anonymous referees for their valuable comments towards improving our manuscript.

---