Allee effect-driven complexity in a spatiotemporal predator-prey system with fear factor

1.   Introduction

The Allee effect is primarily associated with the relationship between per capita growth rate and population density ^[1,2,3,4]. The prevailing situation of the species being threatened due to the Allee effect has drawn major attention to researchers who are studying both the theoretical and practical ecological conservation in exploited environments. Biologists now perceive a diverse range of factors which include the benefit of reproduction, cooperative breeding ^[5,6,7], anti-predator character among the prey ^[8,9], foraging ability and environmental conditioning as well as the purpose of the influence of the Allee effect on the growth of prey ^[10,11,12]. A few reasons behind the emergence of the Allee effect may be stated as difficulties in mate finding, reproductive facilitation, predation, environmental conditioning and inbreeding depression ^[4,13]. In recent years, many eminent researchers have paid considerable attention to investigate the impact of the Allee effect in the growth of prey and the dynamical complexity of the predator-prey model ^{[14,15,16,17]}. The Allee effect is usually categorized by the manifestation of density dependence in the event of its relatively lower range. The population experiences a strong Allee effect resulting in the existence of a critical density where the population growth becomes negative. However, a weak Allee effect causes the population to have a reduced per capita growth rate at low population density, therefore exhibiting a positive per capita growth rate.

A single species logistic model with adjusting for mating encounters was considered, which takes the following form ^[13]:

where $r$ and $K$ represent the intrinsic growth rate and the carrying capacity of the environment, respectively. $h$ is a non-negative constant, which depicts the population density at which the probability of mating is half. The term $\frac{h}{u+h}$ is the probability of not mating, and the term $\frac{\mu h u}{u+h}$ represents the reduction of reproduction due to a mating shortage. Moreover, if $0 < \mu < r$, it is called a weak Allee effect, and if $\mu > r$ then it has the strong Allee effect. Keeping the special aspects in mind, we rewrote the classical single species model (1.1) with additive Allee effect by using the following notation:

where $a, b$ and $c$ are positive constants that represent the birth rate, natural death rate and the death rate induced by intra-species competition, respectively. To obtain (1.1), we only need to set $r = a-b$, $K = \frac{a-b}{c}$ and $m = \mu h$.

Recently, Zanette et al. ^[8] found that the fear of predators could reduce about $40%$ of offspring that song sparrows could produce, the predation risk had an important effect on both the birth and survival rate of offspring. Based on the field observations in ^[8], Wang et al. ^[9] proposed a predator-prey model incorporating the fear effect in prey species as follows:

where the parameter $k$ reflects the level of fear, $f:\mathbb{R}_+\rightarrow \mathbb{R}_+$ is the functional response of predators where $\mathbb{R}_+$ represents nonnegative real number set, $d$ is the natural death rate of the predators, $e$ is the conversion rate of prey's biomass to predator's biomass and the function $f_0$ accounts for the cost of anti-predator defense due to fear and satisfies the following properties:

Moreover, their analyses showed that high levels of fear can stabilize the predator-prey system. However, relatively low levels of fear can induce multiple limit cycles through subcritical Hopf bifurcations, which is different from the classical predator-prey models, ignoring the fear factor where Hopf bifurcations are typically supercritical.

Also, field observations show that fear can create an Allee effect ^[18]. Combining (1.2) and (1.3), we obtained a temporal predator-prey model with Allee and fear effects in prey species and density-dependent death factor of the predator species as follows:

where the parameters $a, b, c, d, e, h, m$ and $k$ have the same meanings as in (1.2) and (1.3). $p$ represents predation factor of the predator species, $s$ is the density-dependent death rate of the predators and $g(u, v)$ denotes the functional response between predator and prey, which usually satisfies the following hypothesis:

$\mathrm{(H_0)}$ $g\in{C^1(\mathbb{R}^2, \mathbb{R})}$ satisfies $g$ and is non-negative, and $g(0, v) = 0$ for all $v \geq 0$. Furthermore, there exists a non-negative continuous function $\phi(u)$ and a polynomial function $\psi(v)$ (of arbitrary degree) such that $g(u, v)\leq \phi(u)\psi(v)$ for any $(u, v)\in \mathbb{R}_+^2$, where $\mathbb{R}_+ = [0, \infty).$

About the hypothesis $\mathrm{(H_0)}$, the typical forms of the functional response $g(u, v)$ are as follows:

${\bf{(i)}}$ Holling type Ⅰ $g(u, v) = u$, type Ⅱ $g(u, v) = \frac{u}{u+\kappa}$, type Ⅲ $g(u, v) = \frac{u^2}{u^2+\kappa}$;

${\bf{(ii)}}$ Beddington-DeAngelis type $g(u, v) = \frac{u}{1+\kappa_1u+\kappa_2v}$;

${\bf{(iii)}}$ ratio-dependent type $g(u, v) = \frac{u}{\kappa_1u+\kappa_2 v}$;

${\bf{(iv)}}$ Crowley-Martin type $g(u, v) = \frac{u}{(1+\kappa_1u)(1+\kappa_2 v)}$.

Throughout the manuscript, attention is focused on the particular form of fear factor defined by $f_0(k, v) = \frac{1}{1+kv}$ ^[9]. By incorporating this function into the system of equations of (1.4), one obtains

Considering the effect of spatial diffusion, one gets the corresponding spatiotemporal model which is given as follows:

Here, $u({\bf x}, t)$ and $v({\bf x}, t)$ are densities of the species $u$ and $v$ on spatial location ${\bf x}$ and time $t$, respectively. $\Omega\subset\mathbb{R}^2$ is a bounded domain with a smooth boundary having no-flux boundary conditions, where $\mathbb{R}$ represents the set of real number. $d_u$ and $d_v$ are diffusion coefficients that represent the respective diffusive velocity of the two species.

In fact, the effect of the fear factor and multiple Allee effects in a temporal or spatiotemporal predator-prey system have been investigated in some literatures (cf. ^[19,20,21]). The research investigations conducted by Li et al. ^[20] and Shi et al. ^[21] are primarily based on the original model established in ^[19], where it was assumed that the fear factor impacts the entire growth rate of the prey species. In ^[20], a one-parameter saddle-node bifurcation and degenerate Hopf bifurcation, as well as two-parameter cusp bifurcation and Bogdanov-Takens bifurcation, were studied in a predator-prey system with the Leslie-Gower functional response. In ^[21], the existence and non-existence of non-constant positive steady states and codimension one Hopf bifurcation were investigated in a diffusive predator-prey model with multiple Allee effects and the fear factor. However, they didn't consider the combination effect of Allee and fear effects on Turing instability and pattern selection.

Recently, one-parameter bifurcation phenomena, including those within a predator-prey model with the Allee effect and the fear effect, were investigated by Lai et al. ^[22]. In this model, the Allee effect was suggested to be an additive, and the fear factor was assumed to influence the entire growth rate of the prey species. Their findings revealed that the density of predator species decreased as the fear factor increased. Notably, the manner in which the fear factor affected prey differed from previously documented approaches in the existing literature ^[19,22]. In actuality, field experiments have demonstrated that the fear factor leads to a reduction in prey species production ^[8]. Considering the long-term shift in prey populations, several experimental observations ^[9,23,24] indicated that the fear factor significantly impacted interactions between prey and predator species. Frightened song sparrow nestlings, for instance, were unable to survive due to a scarcity of food supplies caused by fear. Long-term evolution resulted in substantial changes in both behaviour and physiology driven by fearful predation rather than actual predation. Prey species, in an effort to minimize predation rates, would migrate from high-risk to low-risk zones. Examples of significantly reduced fear-driven reproduction rates include mule deer versus mountain lions, elk versus wolves, snowshoe hares versus dogs ^[26] and dugongs versus sharks ^[25], among others. Wang et al. ^[9] developed a mathematical model to investigate the reduction of the reproduction rate in the presence of fear and to explore the stability mechanism of the system. Subsequently, the influence of fear in predator-prey models, along with the inclusion of prey refuge, was studied by Wang et al. ^[27] and Samaddar et al. ^[28]. In ^[29], authors explored the combined effects of the Allee effect, fear factor and prey refuge on the dynamics of a fractional order prey-predator system. Liu, in the context of a spatiotemporal prey-predator model that incorporated fear and Allee effects, examined the existence and non-existence of non-constant positive solutions and bifurcations, including Hopf and steady-state bifurcations ^[30]. Most recently, Pal et al. investigated pattern formation in a cross-diffusive Leslie-Gower predator-prey model with fear and Allee effects ^[31], although they did not delve into pattern selection and codimension two bifurcations.

The novelty of this study lies in the comprehensive exploration of the bifurcation phenomena, encompassing transcritical bifurcation, saddle-node bifurcation, and Hopf bifurcation, each of codimension one, extending to cusp and Bogdanov-Takens bifurcations of codimension two. These bifurcations are induced by both the Allee effect and the fear factor (refer to Theorems 3.3, 3.5, 3.6 and 3.9 below). Additionally, we investigated the Turing instability mechanism for the spatiotemporal system and interpreted the transition of patterns among three classical types: Hexagonal patterns, stripe patterns, and their combinations. It has been observed that the Allee effect plays a crucial role in inducing Turing instability. Furthermore, one may find that a small Allee effect factor stabilizes the considered system, whereas a large Allee effect factor destabilizes the unique positive coexistence steady state by inducing Turing instability (see Theorem 4.1 below), steady-state bifurcation (see Theorem 4.2 below), and spatiotemporal Hopf bifurcation (see Theorem 4.3 below). Given the significance of these factors, the primary objectives of this investigation revolve around ecological concerns:

● How does the Allee effect and fear factor influence the dynamics of the temporal/spatiotemporal system?

● What are the parametric criteria for stability and instability for both the spatiotemporal interacting species?

● Does intra-specific competition among predators encourage reduced biodiversity?

● What effects do the codimension two cusp and Bogdanov-Takens bifurcations around the coexistence equilibrium have on the environment?

● Is the Allee effect a key mechanism for inducing Turing instability, and if so, how does it impact the spatially heterogeneous distribution of species?

A thorough analysis of stability and bifurcation is successfully carried out in the proposed model system in order to address all the issues mentioned above. The present article is organized as follows: Section $2$ deals with the fundamental properties regarding global classical solvability, positivity, and uniformly boundedness of the spatiotemporal model entities for a very generic situation. The qualitative behavior of the temporal model with linear functional response is highlighted in Section $3$ in terms of the occurrence of various bifurcations. Subsequently, the spatiotemporal model is analyzed and the corresponding numerical simulations are performed to illustrate the effectiveness of the theoretical findings in Section $4$. Finally, the concluding remarks of the entire outcomes of the present model system are included in the concluding Section $5$.

2.   Global classical solvability of the spatiotemporal model

This section consists of some basic properties such as positivity, uniqueness, and boundedness of solutions of the spatiotemporal model in order to verify that the considered model is well-posedness.

Let $\rho\in(n, \infty)$. It is known that $W^{1, \rho}(\Omega, \mathbb{R}^2)$ is continuously embedded in $C(\bar\Omega, \mathbb{R}^2)$. One may then consider the solution of system (1.6) in $S = \{\omega\in W^{1, \rho}(\Omega, \mathbb{R}^2)\big{|}\frac{\partial w}{\partial \nu} = 0\quad \partial\Omega\}$. In view of Amann's theory (see Theorems 1, 14.4, 14.6, 14.7 and 17.1) ^[32,33,34]) and comparison principle of the parabolic equation, one can easily have the following:

Theorem 2.1. Assume that all system parameters are positive and $\mathrm{(H_0)}$ is satisfied. Furthermore, suppose that the initial data $(u_0({\bf x}), v_0({\bf x}))\in S$ for some $\rho > n$ and is nonnegative. Then there exists a maximal time $T_{max}$ such that the considered system (1.6) has a unique bounded nonnegative solution defined on $\Omega\times(0, T_{max})$ fulfilling $(u, v)\in C([0, T_{max}), S)\bigcap C^{2, 1}(\bar\Omega\times(0, T_{max}), \mathbb{R}^2)$. Furthermore, it exists globally in time.

Proof. Since the diffusion coefficient matrix $D = diag\{d_u, d_v\}$ is a positive definite and diagonal matrix, (1.6) is normally elliptic ^[32]. Since $g\in C^1$, then the local existence, uniqueness, and regularity of the solution follows directly according to Theorem 14.4 and Theorem 14.6 ^[34]. Since $u_0({\bf x})\geq 0$ and $v_0({\bf x})\geq 0$, by the comparison principle of parabolic equation, $u({\bf x}, t)\geq 0$ and $v({\bf x}, t)\geq 0$ for all $({\bf x}, t)\in\Omega\times[0, T_{max})$. To show the global existence, one needs to prove the boundedness of $\|u(\cdot, t)\|_{L^\infty(\Omega)}$ and $\|v(\cdot, t)\|_{L^\infty(\Omega)}$ for $t\in[0, T_{max})$ according to Amann's theory (see Theorems 5.2 and 17.1) ^[33,34]. Since the non-negativity of $u$, $v$ and $g(u, v)$, one must conclude from the first equation of (1.6) that

which implies that $u({\bf x}, t)\leq \max\{\sup_\Omega(u_0({\bf x}), \frac{a-b}{c})\}\triangleq M_u$ for all $({\bf x}, t)\in \Omega\times[0, \infty)$, according to the comparison principle of the parabolic equation. On the other hand, since $g$ is non-negative continuous and $u$ is non-negative and bounded, we conclude that there always exists a constant $C_0 > 0$ such that

for $u, v\geq 0$. Since $u$ is a priori bounded, it follows from Theorem 1 ^[35] that there exists a constant $M_v > 0$ such that $v({\bf x}, t)\leq M_v$ for all $({\bf x}, t)\in \Omega\times[0, T_{max})$. Therefore, the solution of (1.6) exists globally in time. This completes the proof of the theorem.

3.   Qualitative analysis of the corresponding non-dimensionalized temporal model: Holling type-Ⅰ functional response

In this section, for convenience of analysis, it is chosen $g(u, v)\equiv u$, then temporal model (1.5) turns into

By using the following non-dimensional transformations

one gets the corresponding non-dimensionalized form (after dropping tildes) given as follows:

where $U_0 = \frac{cu_0}{a}$, $V_0 = \frac{cv_0}{ea}$, and

3.1.   Existence of equilibria of the non-dimensionalized temporal system

Equilibrium points are the of intersection of the zero growth prey curve and zero growth predator curve in the non-negative quadrant of the $uv$-plane. The biologically significant equilibrium points of (3.2) must satisfy the following equations

By solving the second equation of (3.3), we obtain the non-negative roots $v = 0$, $v = \frac{1}{\sigma}(-\delta+\eta u)$, provided $u > \frac{\delta}{\eta}$. Clearly, (3.2) has the trivial equilibrium point $E_{0} = (0, 0)$.

3.1.1.   Semi-trivial equilibrium point

To obtain the existence of arbitrary semi-trivial equilibrium point $(u, 0)$, substituting $v = 0$ into the first equation of (3.3) yields the following equation with respect to $u$

where $A_1 = \beta+\rho-1$ and $A_2 = \gamma-(1-\beta)\rho$. By analyzing the distribution of positive root of (3.4), one may obtain that the existence of the semi-trivial equilibrium point $(u, 0)$ has the following cases:

Case-ⅰ: If $A_{2} < 0$, that is, $\gamma < (1-\beta)\rho$, then (3.4) has a unique positive real root, $u_{1} = \frac{-A_{1}+\sqrt{A_{1}^{2}-4A_{2}}}{2}$, which means that (3.2) has only one semi-trivial equilibrium point $E_{1} = (u_{1}, 0)$;

Case-ⅱ: If $A_{2} = 0$, and $A_{1} < 0$, that is, $\gamma = (1-\beta)\rho$, and $\beta + \rho < 1$, then (3.4) has a unique positive root $u_2 = -A_{1}$. Therefore, (3.2) has only one semi-trivial equilibrium point $E_{2} = (u_{2}, 0)$;

Case-ⅲ: If $A_{2} > 0$, $A_{1} < 0$ and $A_{1}^{2}-4A_{2} > 0$, that is, $\gamma > (1-\beta)\rho$, $\beta+\rho < 1$ and $(\beta-\rho)^{2}+1 > 2(\beta-\rho+2\gamma)$, then (3.4) has two positive roots as $u_{3, 4} = \frac{-A_{1}\pm \sqrt{A_{1}^{2}-4A_{2}}}{2}$, and thus (3.2) has two semi-trivial equilibrium points $E_{3} = (u_{3}, 0)$ and $E_{4} = (u_{4}, 0)$;

Case-ⅳ: If $A_{2} > 0$, $A_{1} < 0$ and $A_{1}^{2}-4A_{2} = 0$, that is, $\gamma > (1-\beta)\rho$, $\beta+\rho < 1$ and $(\beta-\rho)^{2} = 2(\beta-\rho+2\gamma)-1$, then (3.4) has only one positive real root, $u_{5} = -\frac{A_{1}}{2}$. Therefore, (3.2) has a unique semi-trivial equilibrium point $E_{5} = (u_{5}, 0)$.

3.1.2.   Coexistence equilibrium point

To show the existence of the coexistence equilibrium point, substituting $v = \frac{1}{\sigma}(\eta u-\delta)$ into the first equation of (3.3) and simplifying it, we obtain the following cubic equation of one variable:

where,

Let us choose $z = u+a_{1}$, then (3.5) takes the form as follows:

with, $b_{1} = 3(a_{2}-a_{1}^{2}), b_{2} = 2a_{1}^{3}-3a_{1}a_{2}+a_{3}$. Now, by analyzing the distribution of the positive roots of $f(z)$, we obtain the distribution of coexistence equilibrium points of (3.2), which can be described as the following lemma.

Lemma 3.1. For the proposed model system (3.2) , the number of coexistence equilibrium points can be described as follows:

(i) If $\Delta = \frac{b_1^3}{27}+\frac{b_2^2}{4} > 0$ and $a_{3} < 0$, then (3.2) has a unique coexistence equilibrium point given by $E_{*} = (u_{*}, v_{*})$ provided $u^{*} < \frac{\delta}{\eta}$.

(ii) If $\Delta = \frac{b_1^3}{27}+\frac{b_2^2}{4} < 0$ and $b_1 < 0$, then (3.2) has at most three coexistence equilibrium points denoted by, $E_{*}^{1} = (u_{*}^1, v_{*}^{1})$, $E_{*}^{2} = (u_{*}^{2}, v_{*}^{2})$ and $E_{*}^{3} = (u_{*}^{3}, v_{*}^{3})$, more precisely,

$(\boldsymbol{b}\mathit{\pmb{1:}})$ $E_{*}^{1}, E_{*}^{2}$ and $E_{*}^{3}$ coexist if $a_{3} < 0$, $a_{1}+\sqrt{a_{1}^{2}-a_{2}} < 0$ holds.

$(\boldsymbol{b}\mathit{\pmb{2:}})$ Only two of $E_{*}^{1}, E_{*}^{2}$ and $E_{*}^{3}$ coexist if one of the following conditions attached here follows:

$(\boldsymbol{b}\mathit{\pmb{21}})$ $a_{3} > 0$, $a_{1}-\sqrt{a_{1}^{2}-a_{2}} < 0$ and $(\boldsymbol{b}\mathit{\pmb{22}})$ $a_{3} = 0$, $a_{1}+\sqrt{a_{1}^{2}-a_{2}} < 0$.

$(\boldsymbol{b}\mathit{\pmb{3:}})$ Only one of $E_{*}^{1}, E_{*}^{2}$ and $E_{*}^{3}$ exists if one of the following conditions holds: $(\boldsymbol{b}\mathit{\pmb{31}})$ $a_{3} < 0$, $-a_{1}+\sqrt{a_{1}^{2}-a_{2}} < 0$ and $(\boldsymbol{b}\mathit{\pmb{32}})$ $a_{3} = 0$, $a_1 > 0$, $a_{1}-\sqrt{a_{1}^{2}-a_{2}} < 0$.

(iii) If $\Delta = \frac{b_1^3}{27}+\frac{b_2^2}{4} = 0$ and $b_1 < 0$, then the model system (3.2) has at most two positive equilibrium points, namely, $E_{*}^{4} = (u_{*}^{4}, v_{*}^{4}), E_{*}^{5} = (u_{*}^{5}, v_{*}^{5})$ or $\hat{E}_{*}^{4} = (\hat{u}_{*}^{4}, \hat{v}_{*}^{4}), \hat{E}_{*}^{5} = (\hat{u}_{*}^{5}, \hat{v}_{*}^{5})$. More specifically,

$(\boldsymbol{c}\mathit{\pmb{1:}})$ Let us suppose that $b_{2}-2(a_{1}^{2}-a_{2})^{\frac{3}{2}} = 0$, then both $E_{*}^{4}: = E_{*}^1 = E_{*}^2$ and $E_{*}^{5}: = E_{*}^3$ coexists if $a_{2} > 0$ and $-\frac{2}{\sqrt{3}}\sqrt{a_{2}} < a_{1} < -\sqrt{a_{2}}$ hold, where $E_{*}^{4}$ has the multiplicity two, and only $E_{*}^{4}$ exists if one of the following conditions holds: $(\boldsymbol{c}\mathit{\pmb{11}})$ $a_{2} > 0$, $a_{1} \leq -\frac{2}{\sqrt{3}}\sqrt{a_{2}}$, $(\boldsymbol{c}\mathit{\pmb{12}})$ $a_{2} < 0$.

$(\boldsymbol{c}\mathit{\pmb{2:}})$ Let us consider $b_{2}+2(a_{1}^{2}-a_{2})^{\frac{2}{3}} = 0$, then both $\hat{E}_{*}^{4}$ and $\hat{E}_{*}^{5}$ coexists if $a_{2} > 0$ and $a_{1}+\sqrt{a_{2}} < 0$ hold, where $\hat{E}_{*}^{4}$ has the multiplicity two, and only $\hat{E}_{*}^{5}$ exists if one of the following conditions holds: $(\boldsymbol{c}\mathit{\pmb{21}})$ $a_{2} \leq 0$, $a_{1}^{2}+a_{2}^{2} \neq 0$, $(\boldsymbol{c}\mathit{\pmb{22}})$ $a_{1} > \frac{2}{\sqrt{3}}\sqrt{a_{2}}$.

Proof. Denote by $\Delta = (\frac{b_1}{3})^3+(\frac{b_2}{2})^2 = \frac{b_1^3}{27}+\frac{b_2^2}{4}$ the discriminant of cubic equation (3.6). By employing Cardano's formula, and Descartes' rule of signs ^[36], we have three cases:

$(i)$ If $\Delta > 0$, i.e., $\frac{b_1^3}{27}+\frac{b_2^2}{4} > 0$, then the cubic equation (3.6) has a unique real root denoted by $z_{*}$. For this case (3.5) has a unique real root $u_{*}$, which is positive if $a_{3} < 0$ follows. Furthermore, according to the Cardano's formula, the unique positive real root $u_\ast$ is given by

$(ii)$ If $\Delta < 0$ and $b_1 < 0$, i.e., $\frac{b_1^3}{27}+\frac{b_2^2}{4} < 0$ and $a_1^2-a_2 > 0$, according to the Cardano's formula, then (3.5) has three distinct real roots given by, $u_{*}^{1}, u_{*}^{2}$ and $u_{*}^{3}$, respectively, where

More precisely, since $\theta\in[0, \pi]$, one obtains $(u_{*}^1)_{min} = -a_1-\sqrt{a_1^2-a_2}$, $(u_{*}^2)_{min} = -a_1+\sqrt{a_1^2-a_2}$, and $(u_{*}^3)_{min} = -a_1-2\sqrt{a_1^2-a_2}$. $(u_{*}^1)_{max} = -a_1+\sqrt{a_1^2-a_2}$, $(u_{*}^2)_{max} = -a_1+2\sqrt{a_1^2-a_2}$, and $(u_{*}^3)_{max} = -a_1-\sqrt{a_1^2-a_2}$. Thus, one can conclude the following facts.

If $a_1+2\sqrt{a_1^2-a_2} < 0$ then $(u_{*}^1)_{min} > 0$, $(u_{*}^2)_{min} > 0$, and $(u_{*}^3)_{min} > 0$, which implies that $({{\bf{b1}}})$ is satisfied;

If $a_1-\sqrt{a_1^2-a_2} < 0$, then $(u_{*}^2)_{min} > 0$, thus $u_{*}^2 > 0$. Since $a_3 > 0$, it follows from Descartes' rule of signs that (3.5) has two positive real roots, which implies (b21) holds. On the other hand, if $a_3 = 0$, then (3.5) has at most two positive real roots. Since $a_1+\sqrt{a_1^2-a_2} < 0$, it follows that $(u_{*}^1)_{min} > 0$ and $(u_{*}^2)_{min} > 0$, which manifests that (b22) holds.

If $a_3 < 0$, then (3.5) has at least one positive real roots. Since $-a_1+\sqrt{a_1^2-a_2} < 0$, one obtains $a_1 > 0$ and $a_2 > 0$. Thus, (3.5) has exactly one positive real root. If $a_3 = 0$, $a_1 > 0$ and $a_1-\sqrt{a_1^2-a_2} < 0$, one obtains $a_2 < 0$, and thus, (3.5) has only one positive real root, which implies that the conclusion (b3) holds.

$(iii)$ If $\Delta = 0$ and $b_1 < 0$, i.e., $\frac{b_1^3}{27}+\frac{b_2^2}{4} = 0$ and $b_{1} < 0$, then (3.5) has three real roots, one of them a double roots denoted by $E_{*}^{4}$ and the other denoted by $E_{*}^{5}$. More specifically, according to Cardano's formula, we have $u_{*}^{4}: = u_{*}^{1} = u_{*}^{2} = \sqrt{a_{1}^{2}-a_{2}}-a_{1}$, $u_{*}^{5}: = u_{*}^{3} = -2\sqrt{a_{1}^{2}-a_{2}}-a_{1}$ when we set $b_{2}-2(a_{1}^{2}-a_{2})^{\frac{3}{2}} = 0$ in formula (3.7). Also, it is to be noted that $u_{*}^{4} > u_{*}^{5}$, so $u_{*}^{4}$ and $u_{*}^{5}$ are positive if $a_{1}+2\sqrt{a_{1}^{2}-a_{2}} < 0$ and $a_{1}^{2}-a_{2} > 0$ hold. These two conditions are identical to $a_{2} > 0$ and $-\frac{2}{\sqrt{3}}\sqrt{a_{2}} < a_{1} < -\sqrt{a_{2}}$, whereas only $u_{*}^{4}$ is positive when $a_{2} > 0$, $a_{1}+\frac{2}{\sqrt{3}}\sqrt{a_{2}} \leq 0$, or $a_{2} < 0$, which implies that (c1) holds.

When $b_{2}+2(a_{1}^{2}-a_{2})^{\frac{3}{2}} = 0$ in (3.7), one has $\hat{u}_{*}^{4}: = u_{*}^{1} = u_{*}^{3} = -\sqrt{a_{1}^{2}-a_{2}}-a_{1}$, $\hat{u}_{*}^{5}: = u_{*}^{2} = 2\sqrt{a_{1}^{2}-a_{2}}-a_{1}$. Since, $\hat{u}_{*}^{5} > \hat{u}_{*}^{4}$, it is to be noted that if $a_{2} > 0$, $a_{1}+\sqrt{a_{2}} < 0$, then one has that both $\hat{u}_{*}^{5}$ and $\hat{u}_{*}^{4}$ are positive. Moreover, there are two sets of conditions, $a_{2} \leq 0, a_{1}^{2}+a_{2}^{2} \neq 0$ and $a_{1}-\frac{2}{\sqrt{3}}\sqrt{a_{2}} > 0$, that can result in $\hat{u}_{*}^{5} > 0$, and thus, the conclusion (c2) holds.

Based on the analysis above, one may easily prove Lemma 3.1, which ends the proof of the theorem.

Remark 3.2. Lemma 3.1 provides an upper bound estimation on the number of coexistence equilibrium points $E_{*}(u_{*}, v_{*})$. One has to impose the condition $u_\ast\leq\frac{\delta}{\eta}$ to ensure the existence of the exact number of coexistence equilibrium points. Additionally, one may obtain the range of original model parameters, which ensures the existence of coexistence equilibrium points by substituting the coefficients $a_1, a_2$ and $a_3$ of (3.5) into the expressions $b_1$ and $b_2$.

In Figure 1, graphical representation of the system equilibria are provided for various system parameter values through Lemma 3.1. Especially, this implies that the Allee effect and fear factor can affect the number of interior equilibrium points, as well as the results in a complex distribution of the population in the habitat.

3.2.   Stability and bifurcation analysis of trivial and semi-trivial equilibrium points of the non-dimensionalized temporal system

This subsection consists of various dynamical behavior of the temporal model system (3.2) around the trivial equilibrium point $E_0$ as well as the semi-trivial equilibrium points $E_i$ for $i = 1, 2, 3, 4, 5$, respectively. The Jacobian matrix of (3.2) at any points $E(u, v)$ in the $uv$-plane is given by

Theorem 3.3. The trivial equilibrium point $E_{0} = (0, 0)$ is

(i) a hyperbolic attractor if ${\gamma} > \rho(1-\beta)$,

(ii) a hyperbolic saddle point if ${\gamma} < \rho(1-\beta)$.

Proof. The proof is trivial, so it is omitted here.

By direct calculation, we have the following:

Theorem 3.4. The model system (3.2) undergoes a transcritical bifurcation around the trivial equilibrium point $E_{0} = (0, 0)$ as the system parameter $\gamma$ crosses the threshold value $\gamma = \gamma^{[TC]} = \rho(1-\beta)$, provided $\beta+\rho < 1$.

Proof. Since the occurrence of transcritical bifurcation around the trivial equilibrium point $E_{0} = (0, 0)$, one of the eigenvalues of the corresponding Jacobian matrix must be zero, which implies that $\det J_{E_{0}} = 0$. One may easily verify that at the critical system parameter value $\gamma^{[TC]} = \rho(1-\beta)$, we get $\det J_{E_{0}} = 0$.

Let $\mathcal{V}$ and $\mathcal{W}$ be the eigenvectors of the Jacobian matrix $J_{E_{0}}$ and its transpose matrix $J_{E_{0}}^{T}$, respectively, corresponding to the simple zero eigenvalue. Then, at the critical parameter value $\gamma^{[TC]} = \rho(1-\beta)$, we obtain the eigenvectors as $\mathcal{V} = (v_{1}, v_{2}) = (1, 0)^{T}$ and $\mathcal{W} = (w_{1}, w_{2}) = (1, 0)^{T}$. Now, we compute the transversality conditions $\Delta_{1}, \Delta_{2}$ and $\Delta_{3}$ as defined below:

Thus, if $\beta+\rho < 1$ then one obtains $\Delta_{3}: \mathcal{W}^{T}[D^{2}F(0, 0; \gamma^{[TC]})(\mathcal{V}, \mathcal{V})] \neq 0$. Therefore, (3.2) undergoes a transcritical bifurcation around the trivial equilibrium point $E_{0} = (0, 0)$ whenever the parameter $\gamma$ attains the threshold value $\gamma^{[TC]} = \rho(1-\beta)$, according to Sotomayor's theorem ^[37].

Direct applications of Routh-Hurwitz criteria give the following:

Theorem 3.5. The semi-trivial equilibrium points $E_{i} = (u_{i}, 0)$ with $i = 1, 2, 3, 4, 5$ are locally asymptotically stable if ${\gamma} < (u_i+\rho)^2$ and $u_i < \frac{\delta}{\eta}$. Additionally, each $E_{i} = (u_{i}, 0)$ for $i = 1, 2, 3, 4, 5$ is saddle point if $\gamma > (u_i+\rho)^2$ and $u_i < \frac{\delta}{\eta}$ or $\gamma < (u_i+\rho)^2$ and $u_i > \frac{\delta}{\eta}$.

Theorem 3.6. The model system (3.2) undergoes a transcritical bifurcation around the semi-trivial equilibrium point $\bar{E} = (\bar{u}, 0)$ for the critical system parameter value $\gamma = \gamma^{[TC]}$ when $\gamma$ satisfies $\delta-\eta\bar{u} = 0$ under the two parametric restrictions

where $\bar{u}$ is determined by (3.4).

Proof. To obtain the transcritical bifurcation requires that one of the eigenvalues of the corresponding Jacobian matrix at the semi-trivial equilibrium point $\bar{E}$ must be zero, which implies that $\det(J_{\bar{E}}) = 0$. We accordingly obtained the critical system parameter value $\delta = \delta^{[TC]} = \eta \bar{u}$.

Next we checked Sotomayor's condition ^[37]. Let $\mathcal{A}$ and $\mathcal{B}$ be the eigenvectors of the Jacobian matrix $J_{\bar{E}}$ and its transpose matrix $J_{\bar{E}}^{T}$, respectively, corresponding to the simple zero eigenvalue. At the critical system parameter value $\gamma = \gamma^{[TC]}$, we obtain the eigenvectors $\mathcal{A}$ and $\mathcal{B}$ as $\mathcal{A} = (A_{1}, A_{2})^{T}\triangleq(\bar{J}_{12}, -\bar{J}_{11})^{T}$ and $\mathcal{B} = (B_{1}, B_{2})^{T}\triangleq(0, -\bar{J}_{12})^{T}$. Now, we compute the three numerical values $\Delta_{1}, \Delta_{2}$ and $\Delta_{3}$ as follows:

So, if ${\gamma}\neq (\bar{u}+\rho)^{2}$, then we have $\Delta_{2}: \mathcal{B}^{T}[DF_{\delta}(\bar{u}, 0; \delta^{[TC]})(\mathcal{A})] \neq 0$; if $\sigma \neq \frac{\eta A_{1}}{A_{2}}$, that is, $\gamma\neq \frac{(\eta(\eta+\alpha)+\sigma)(\bar{u}+\rho)^2}{\sigma}$, then we have $\Delta_{3}: \mathcal{B}^{T}[D^{2}F(\bar{u}, 0; \delta^{[TC]})(\mathcal{A}, \mathcal{A})] \neq 0.$ Hence, a transcritical bifurcation takes place into the system surrounding the semi-trivial equilibrium point $\bar{E} = (\bar{u}, 0)$ whenever the system parameter $\delta$ crosses the threshold value $\delta^{[TC]} = \eta \bar{u}$ under the pre-assumptions as defined in (3.8).

Remark 3.7. Biological interpretation: From Theorem 3.3 to Theorem 3.6, we obtained that large Allee effects are prone to causing both the prey and predator population extinction, and small Allee effects can only cause the predator species' extinction. However, the prey species $u$ persists even when the prey species density is low.

3.3.   Bifurcation analysis of the non-dimensionalized temporal system around the coexistence equilibrium point

This section serves various types of bifurcation phenomena and takes place in (3.2) around the coexistence equilibrium point $E_{*}(u_{*}, v_{*})$ for some specific system parameter values. We showed that the temporal system undergoes codimension one bifurcations including saddle-node bifurcation, Hopf-bifurcation and codimension two cusp and Bogdanov-Takens bifurcations. As we know, at the coexistence equilibrium point $E_{*} = (u_{*}, v_{*})$, the characteristic equation of linearized system is given by

where $T = \mathrm{tr}(J_{E_{*}}) = J_{11}^{\ast}+J_{22}^{\ast}$ and $D = \det(J_{E_{*}}) = J_{11}^{\ast}J_{22}^{\ast}-J_{12}^{\ast}J_{21}^{\ast}$ with

Remark 3.8. System (3.2) may exhibit bi-stability between the trivial equilibrium $E_0$ and the coexistence equilibrium $E_\ast$. For example, if one takes a set of system parameters as $\alpha = 0.03, \beta = 0.01, \gamma = 0.21, \rho = 0.08, \eta = 0.3, \delta = 0.06$ and $\sigma = 0.2$, then it is easy to check that $\gamma-\rho(1-\beta) = 0.1292 > 0$, $\mathrm{tr}(J_{E_\ast(0.3453, 0.4579)}) = -0.036 < 0$ and $\mathrm{det}(J_{E_\ast(0.3453, 0.4579)}) = 0.0105 > 0$ are satisfied. So, the trivial equilibrium $E_0$ and the coexistence equilibrium $E_\ast$ are all stable in this case.

3.3.1.   Existence of Hopf-bifurcation

The following theorem follows directly from the Hopf bifurcation theorem ^[38].

Theorem 3.9. The model system (3.2) exhibits a Hopf-bifurcation and gives birth of a limit cycle around the coexistence equilibrium point $E_{*} = (u_{*}, v_{*})$ at the critical system parameter value $\gamma = \gamma^{[HB]}$ provided by $T(\gamma^{[HB]}) = 0$, $D(\gamma^{[HB]}) > 0$ and $\dot{T}(\gamma^{[HB]}) = \frac{d T}{d \gamma}\rvert_{\gamma = \gamma^{[HB]}} \neq 0$.

3.3.2.   Direction and stability of Hopf-bifurcation

This section deals with the determination of the stability and direction of Hopf-bifurcation, which takes place in (3.2) around the coexistence equilibrium point $E_{*} = (u_{*}, v_{*})$, through the computation of the first Lyapunov number. For this purpose, here we will follow the normal form technique ^[36]. Let us choose the following transformation $Q_{1} = u-u_{*}$ and $Q_{2} = v-v_{*}$, which shift the coexistence equilibrium point $E_{*} = (u_{*}, v_{*})$ to the origin. Using this transformation, (3.2) may be written as:

By Taylor's series expansion of the right hand side of (3.11) about the origin, and we obtain:

where $Q = (Q_{1}, Q_{2})^{T}$ and the coefficients $c_{ij}, d_{ij}$; $i, j = 0, 1, 2, 3, \ldots$ are given by

System (3.12) may be rewritten as

where $\Phi = \begin{pmatrix} c_{20}Q_{1}^{2}+c_{11}Q_{1}Q_{2}+c_{02}Q_{2}^{2}+c_{30}Q_{1}^{3}+c_{12}Q_{1}Q_{2}^{2}+c_{03}Q_{2}^{3} \\ d_{11}Q_{1}Q_{2}+d_{02}Q_{2}^{2} \end{pmatrix} = \begin{pmatrix} \Phi_{1} \\ \Phi_{2} \end{pmatrix}$. According to the argument in ^[7], after transformation, we get that (3.13) is equivalent to the following

where, the terms $\mathcal{N}^{1}$ and $\mathcal{N}^{2}$ are given by;

with,

Now, to observe the stability behavior of the bifurcating limit cycle, one may calculate the first Lyapunov number $\mathcal{L}_{1}$ by computing the following formula:

where, $\mathcal{N}^{k}_{ij} = \frac{\partial^{2} \mathcal{N}^{k}}{\partial y_{i} \partial y_{j}}\bigg{\rvert}_{(y_{1}, \; y_{2}; \; \sigma) = (0, \; 0;\; \sigma^{[HB]})}$, $\mathcal{N}^{k}_{ijl} = \frac{\partial^{3} \mathcal{N}^{k}}{\partial y_{i} \partial y_{j} \partial y_{l}}\bigg{\rvert}_{(y_{1}, \; y_{2}; \; \sigma) = (0, \; 0; \; \sigma^{[HB]})}$. Therefore, depending upon the positivity of $\mathcal{L}_{1}$, one may consider the following cases in order to determine the stability-instability behavior of the limit cycle:

Case-ⅰ: If $\mathcal{L}_{1} > 0 (<0)$, then a sub-critical (super-critical) Hopf-bifurcation takes place into (3.2) or an unstable (stable) limit cycle will be created surrounding the coexistence equilibrium point $E_{*} = (u_{*}, v_{*})$.

Case-ⅱ: Whenever the first Lyapunov number $\mathcal{L}_{1}$ vanishes, then a local bifurcation of co-dimension two, namely, generalized Hopf bifurcation or Bautin bifurcation appears into (3.2) about the coexistence equilibrium point $E_{*} = (u_{*}, v_{*})$.

3.3.3.   Numerical simulations of Hopf bifurcation

In this subsection, we studied the Hopf bifurcation numerically. The parameter values provided for numerical simulations of (3.2) are: $\alpha = 0.38, \beta = 0.01, \gamma = 0.21, \rho = 0.08, \eta = 0.3, \delta = 0.06$ and $\sigma = 0.2$. To solve the nonlinear system (3.2) numerically we have used MATLAB R2016a by constructing a suitable code. In Figure 2, each system parameter is chosen as a bifurcation parameter to illustrate the occurrence of Hopf bifurcation, respectively. All numerical simulations suggest that Hopf bifurcations are subcritical in nature, which can be interpreted by the positive first Lyapunov numbers (Figure 2(ⅰ) with $\mathcal{L}_{1} = 2.901765$; Figure 2(ⅱ) with $\mathcal{L}_{1} = 2.918016$; Figure 2(ⅲ) with $\mathcal{L}_{1} = 2.848071$; Figure 2(ⅳ) with $\mathcal{L}_{1} = 2.905451$; Figure 2(ⅴ) with $\mathcal{L}_{1} = 2.940986$; Figure 2(ⅵ) with $\mathcal{L}_{1} = 2.803728$; Figure 2(ⅶ) with $\mathcal{L}_{1} = 2.863622$).

3.3.4.   Saddle-node bifurcation

In this section, one may discuss the existence of saddle-node bifurcation in (3.2) around the coexistence equilibrium point $E_{*} = (u_{*}, v_{*})$ for the critical system parameter value $\gamma = \gamma^{[SN]}$. For this purpose, here one may use Sotomayor's theorem ^[37].

Theorem 3.10. The model system (3.2) may exhibit a saddle-node bifurcation around the coexistence equilibrium point $E_{*} = (u_{*}, v_{*})$ whenever the system parameter $\gamma$ attains the critical value $\gamma = \gamma^{[SN]}$ and $\mathcal{T}\neq 0$, where

and $\mathcal{T}$ is specified in the proof of the theorem.

Proof. It is easy to check that at the critical system parameter value $\gamma = \gamma^{[SN]}$ the determinant of the corresponding Jacobian matrix $J_{E_{*}}$ vanishes, that is, $\det(J_{E_{*}}) = 0$, which implies that one of the eigenvalues of $J_{E_{*}}$ is equal to zero. Next, we checked the transversality condition. Let $\mathcal{\hat{V}}$ and $\mathcal{\hat{W}}$ be the eigenvectors of the Jacobian matrix $J_{E_{*}}$ and its transpose matrix $J_{E_{*}}^{T}$ respectively corresponding to the simple zero eigenvalue. Then, at the critical system parameter value $\gamma^{[SN]}$, one may calculate the eigenvectors as

Next, we calculated the transversality conditions $\Gamma_{1}$ and $\Gamma_{2}$ as follows:

Now, to compute the second transversality condition, one needs to compute the second order derivatives at the coexistence equilibrium point $E_{*} = (u_{*}, v_{*})$ as follows:

If $\mathcal{T}\neq 0$ by Sotomayor's theorem, one may conclude that (3.2) undergoes a saddle-node bifurcation around the coexistence equilibrium point $E_{*} = (u_{*}, v_{*})$ whenever the system parameter $\gamma$ attains the critical value $\gamma^{[SN]}$.

Example 3.11. Now, to verify the analytical findings of the occurence of saddle-node bifurcation around the coexistence equilibrium point, here an attempt is made to preform a numerical example. For this purpose, the following set of system parameter values are taken into account; $\alpha = 0.38, \beta = 0.01, \rho = 0.08, \eta = 0.3, \sigma = 0.2$ and $\delta = 0.06$. Then, at the critical system parameter value $\gamma^{[SN]} = 0.227779$, (3.2) has a coexistence equilibrium point $E_{*} = (0.258868, 0.088301)$. Also, at that critical value the eigenvectors $\mathcal{\hat{V}}$ and $\mathcal{\hat{W}}$ are given by $\mathcal{\hat{V}} = (-0.169746, -0.254621)^{T}, \, \, \mathcal{\hat{W}} = (0.024690, -0.254621)^{T}$, where ${T}$ represents the transpose of matrix. Next, we calculated the transversality conditions $\Gamma_{1}$ and $\Gamma_{3}$ as follows:

Hence, both of the transversality conditions for the occurence of saddle-node bifurcation holds good. Thus, one may conclude that (3.2) undergoes a saddle-node bifurcation around the coexistence equilibrium point $E_{*} = (0.258868, 0.088301)$ whenever the system parameter $\gamma$ crosses the threshold value $\gamma^{[SN]} = 0.227779$ (cf. Figure 3).

3.3.5.   Bogdanov-Takens bifurcation

When in a system the saddle-node bifurcation curve and the Hopf bifurcation curve meet tangentially, at that point a local bifurcation of codimension two, namely a Bogdanov-Takens bifurcation takes place into the system ^[38]. Here, it is observed that, (3.2) exhibits Hopf-bifurcation for the critical system parameter value $\sigma = \sigma^{[HB]}$ surrounding the coexistence equilibrium point $E_{*} = (u_{*}, v_{*})$. Also, (3.2) undergoes a saddle-node bifurcation around the coexistence equilibrium point $E_{*} = (u_{*}, v_{*})$ whenever the system parameter $\gamma$ attains the threshold value $\gamma^{[SN]}$. Subsequently, a Bogdanov-Takens bifurcation takes place into (3.2) around the coexistence equilibrium point $E_{*} = (u_{*}, v_{*})$.

Theorem 3.12. The temporal system (3.2) experiences a Bogdanov-Takens bifurcation surrounding the coexistence equilibrium $E_{*} = (u_{*}, v_{*})$ for the critical parameter value $(\sigma, \gamma) = (\sigma^{[BT]}, \gamma^{[BT]})$ provided the following conditions must be followed:

where $c_{20}(\epsilon^*), d_{11}(\epsilon^*)\, \, d_{20}(\epsilon^*)$, $\vartheta_1(\epsilon)$, and $\vartheta_2(\epsilon)$ are given in the proof of the theorem. Then, there exists a parameter dependent nonlinear smooth invertible variable transformation, and a direction preserving time reparametrization, which together reduce (3.2) to the following equivalent norm form

where $\vartheta_{1}$ and $\vartheta_{2}$ are defined in (A.17), and $s = \pm 1$.

Proof. For the detailed proof please see the Appendix.

Without loss of generality, let us start with the negative value of $s$, that is, let $s = -1$, then for (3.2) there exists a neighbourhood of $(\vartheta_{1}, \vartheta_{2}) = (0, 0)\in \mathbb R^{2}$, which divides the bifurcation plane into four regions through the following curves according to ^[38]:

Here, $SN$ represents the saddle-node bifurcation curve with two branches $SN^{+}$ and $SN^{-}$, corresponding to the sign of $\vartheta_{2}$ such as $\vartheta_{2} < 0$ and $\vartheta_{2} > 0$, respectively, $H$ represents the Hopf-bifurcation curve and $HL$ represents the Homoclinic bifurcation curve. In addition, the argument against the case for the positive value $s$, that is, for $s = 1$ is very much identical to the argument in favor of $s = -1$. For this case the local representation of bifurcation curves in a very small neighborhood of $(\vartheta_{1}, \vartheta_{2}) = (0, 0)$ were generated using the following linear transformation given by $(\xi_{1}, \xi_{2}, t, \vartheta_{1}, \vartheta_{2}) = (\xi_{1}, -\xi_{2}, -t, \vartheta_{1}, -\vartheta_{2}).$

3.3.6.   Numerical simulations of Bogdanov-Takens bifurcation

In this subsection, we devote to perform some numerical simulations to illustrate the effectiveness of Theorem 3.12. Without loss of generality, we choose $\gamma$ and $\sigma$ as the bifurcation parameters, and take a set of parameter values as follows:

Under this set of parameters, (3.2) has the interior equilibrium point $E_\ast(0.4030, 0.0459)$ and the Jacobian matrix associated with $E_\ast$ has two zero eigenvalues. We now perturb the parameters $\gamma$ and $\sigma$ in the small neighbourhood of 1.324208 and 0.268561, respectively, that is, let $\sigma = 0.268561+\epsilon_1$ and $\gamma = 1.324208+\epsilon_2$, where $\epsilon_1$ and $\epsilon_2$ are small perturbation parameters. Then, we get the system corresponding to the form (A.9) as follows

where $G_{10} = -0.8345\epsilon_2$, $H_{10} = -0.0021\epsilon_1$, $m_{10} = 0.0609-0.3429\epsilon_2$, $n_{10} = 0.0138$, $m_{01} = -0.2688$, $n_{01} = -0.0609-0.092\epsilon_1$, $m_{20} = -1.6187+1.4198*\epsilon_2$, $m_{11} = -0.6671$, $m_{02} = 0.1105$, $n_{20} = 0$, $n_{11} = 0.3$, and $n_{02} = -0.537122-2\epsilon_1$.

Since $d_{11}(\epsilon^*)+2c_{20}^{*} = -5.0018\neq 0$ and $d_{20}(\epsilon^\ast) = -00892\neq 0$, we used the formulae presented in Section 3.3.5 to calculate $\gamma_1(\epsilon), \gamma_2(\epsilon), \mathcal{D}(\epsilon)$, and $\mathcal{E}(\epsilon)$. Then, we obtained the following norm form of the Bogdanov-Takens bifurcation corresponding to (A.16) by using sign calculation of MATLAB R2016a software:

where

Then, the generic norm form of B-T bifurcation for (3.2) near the origin for small $\|\epsilon\|$ is locally topologically equivalent to the system:

where $\vartheta_1(\epsilon)$ and $\vartheta_2(\epsilon)$ can be obtained with the expression of $\gamma_1(\epsilon), \gamma_2(\epsilon), \mathcal{D}(\epsilon)$, and $\mathcal{E}(\epsilon)$ given in (3.20), $s = sign\left(\frac{\theta}{d_{20}(\epsilon^\ast)}\right) = 1$, since $\det\left(\frac{\partial(\vartheta_1, \vartheta_2)}{\partial(\epsilon_1, \epsilon_2)}\big{|}_{(\epsilon_1, \epsilon_2) = (0, 0)}\right) = 1507084.8078\neq 0,$ which means that the transversality condition BT. 3 is satisfied. In view of Theorem 3.12, (3.2) undergoes B-T bifurcation at the coexistence equilibrium $E^{*}$, and the local representation of bifurcation curves in $\epsilon_1, \epsilon_2$ in the small neighborhood of the origin are given in (3.15)–(3.18). This can be seen from the phase portrait diagrams depicted in Figure 4. Similarly, we obtained different B-T bifurcation phase portraits that occurred in the vicinity $(\epsilon_{1}, \epsilon_{2})$ of the chosen bifurcating parameters (see Figure 5). The detailed explanation in each figure corresponding to Bogdanov-Takens bifurcation are listed in Tables 1 and 2.

3.4.   Exploration of codimension one bifurcations of the non-dimensionalized temporal system

In this present section we have explored the occurrence of various codimension one bifurcations, namely the Hopf-bifurcation, transcritical bifurcation and saddle-node bifurcation for some specific system parameter value surrounding the ecologically meaningful equilibrium point of (3.2). In Figure 6, we have illustrated the one-parameter bifurcation diagrams with respect to the system parameter as specified in the $x$-axis of each figure ^[39].

In Figure 6(ⅰ) we have observed that a Hopf-bifurcation took place in (3.2) around the coexistence equilibrium point $E_{*} = (0.352582, 0.228873)$ at the system parameter value $\gamma^{[HB]} = 0.211421$, and we obtained the corresponding first Lyapunov number as $\mathcal{L}_{1} = 2.901765$. Also, saddle-node bifurcation appeared in (3.2) at $\gamma^{[SN]} = 0.227779$ around $E_{*} = (0.258868, 0.088301)$. Additionally, (3.2) exhibited a transcritical bifurcation at the critical system parameter value $\gamma^{[TC]} = 0.2212$ about the axial equilibrium point $\bar{E} = (0.2, 0)$.

The bifurcation diagram, with respect to the system parameter $\delta$ as depicted through the Figure 6(ⅱ), illustrated that (3.2) exhibits a Hopf-bifurcation around the coexistence equilibrium point $E_{*} = (0.350447, 0.233725)$ for the critical system parameter value $\delta^{[HB]} = 0.058389$. At this critical system parameter value we obtained the first Lyapunov number as $\mathcal{L}_{1} = 2.918016$, which implied that the appeared Hopf-bifurcation is subcritical. Additionally, we have observed that a saddle-node bifurcation takes place at (3.2) at the critical system parameter value $\delta^{[SN]} = 0.043514$ surrounding the coexistence equilibrium point $E_{*} = (0.247418, 0.153555)$. Moreover, we have observed that the system experiences a transcritical bifurcation about the predator free equilibrium point $\bar{E} = (0.178911, 0)$.

From the bifurcation diagram of (3.2), with respect to the system parameter $\beta$ as depicted in Figure 6(ⅲ), we have found that the model system exhibits a Hopf-bifurcation at the critical system parameter value $\beta^{[HB]} = 0.014360$ about the coexistence equilibrium point $E_{*} = (0.351251, 0.226876)$. Moreover, a saddle-node bifurcation took place in the system at $\beta^{[SN]} = 0.063588$ surrounding the equilibrium point $E_{*} = (0.244684, 0.067026)$. Additionally, the system experienced a transcritical bifurcation about the predator free equilibrium point $\bar{E} = (0.2, 0)$ at $\beta^{[TC]} = 0.05$.

3.5.   Detection of codimension two bifurcations of the non-dimensionalized temporal system

In this present section we have explored the occurrence of codimension two bifurcations, namely Bogdanov-Takens bifurcation, and cusp bifurcation for some specific system parameter value surrounding the ecologically meaningful equilibrium point by depicting the two-parameter bifurcation diagram ^[39]. It should be noted that we have dealt with a local bifurcation of codimension two, namely the cusp bifurcation in (3.2) due to the variation of two system parameters. It is a well-known fact that when two saddle-node bifurcation curves meet at a point, that is where a cusp bifurcation may take place. Also, one more fact is that when the saddle-node bifurcation curve vanishes on the transcritical bifurcation curve, a cusp bifurcation may arise. This suggests that three equilibrium points may be merged at the cusp bifurcation point due to the variation of two significant system parameters. For more detailed derivation about cusp bifurcation, please refer to Step-1 and Step-2 in the Appendix.

In Figure 7(ⅰ) we have detected two cusp bifurcation points appearing on the saddle-node bifurcation curve, which indicates that (3.2) exhibits a cusp bifurcation at the critical system parameter values $(\eta^{[CP]}, \gamma^{[CP]}) = (0.170821, 0.275520)$ about the coexistence equilibrium point $E_{*} = (0.352761, 0.000693)$. The normal form coefficient is given by $-0.6960690$ ^[39]. Similarly, at $(\eta^{[CP]}, \gamma^{[CP]}) = (0.829643, 0.139782)$ about the coexistence equilibrium $E_{*} = (0.072383, 0.000254)$, the normal form coefficient is given by $-0.3359687$.

It is evident from Figure 7(ⅱ) that, for the critical system parameter value $(\eta^{[BT]}, \sigma^{[BT]}) = (0.818409, 1.580402)$, (3.2) exhibits a Bogdanov-Takens bifurcation about the coexistence equilibrium point $E_{*} = (0.281675, 0.1079)$, which is an intersection point of Hopf bifurcation and saddle-node bifurcation curves. The normal form coefficients are given by $(-0.1878629, -2.409671)$ ^[39]. Furthermore, a cusp bifurcation took place in the system at the system parameter value $(\eta^{[CP]}, \sigma^{[CP]}) = (0.335815, 0.114164)$ about the equilibrium point $E_{*} = (0.179915, 0.002268)$, and the normal form coefficient is given by $-0.9066523$ ^[39].

Similarly, a Bogdanov-Takens bifurcation point has been detected in Figure 7(ⅲ) at the system parameter value $(\beta^{[BT]}, \sigma^{[BT]}) = (0.142135, 1.559992)$. That is, a Bogdanov-Takens bifurcation takes place in the system for this system parameter value about the coexistence equilibrium point $E_{*} = (0.351249, 0.029086)$, and the normal form coefficient is given by $(-0.0408996, -1.974066)$. Also, a cusp bifurcation point has been detected in the Figure at the critical system parameter value $(\beta^{[CP]}, \sigma^{[CP]}) = (0.05, 0.121539)$ about the predator free equilibrium point $\bar{E} = (0.2, 0)$, and the normal form coefficient is given by $-0.218866$.

4.   Stability analysis of the non-dimensionalized spatiotemporal model: Holling type-Ⅰ functional response

In this subsection, stability and Turing instability mechanisms are studied. The corresponding non-dimensionalized form of the spatiotemporal model (1.6) with Holling type-Ⅰ functional response (i.e., $g(u, v)\equiv 1$) is given as follows:

The linearized system of (4.1) can be written as

Here, $J_{ij}^\ast$ for $i, j = 1, 2$ are the same as those in (3.10). Let us suppose $\Omega = [0, L]\times[0, L]$. Then, $k^2 = (\frac{m\pi}{L})^2+(\frac{n\pi}{L})^2$ with $(m, n)\in\mathbb{N}^2\setminus\{(0, 0)\}$ as the eigenvalues of $-\nabla^2$ on $\Omega$ under homogeneous Neumann boundary conditions. The corresponding eigenfunctions are given by ${\bf{w}}_{m, n}({\bf{x}}) = {\bf{C}}_{m, n}\cos\frac{m\pi x}{L}\cos\frac{n\pi y}{L}$, where ${\bf{x}} = (x, y)$ and ${\bf{C}}_{m, n}$ are Fourier coefficients with $m$ and $n$ as integers. Denote

It follows that the eigenvalues of $L(d)$ are given by the eigenvalues of $L_k(d)$. So, the characteristic equation of (4.2) is

with $A(k^2) = (1+d)k^2-(J_{11}^\ast+J_{22}^\ast), \, \, B(k^2) = dk^4-(d J_{11}^\ast +J_{22}^\ast)k^2+J_{11}^\ast J_{22}^\ast-J_{12}^\ast J_{21}^\ast.$

4.1.   Turing instability induced by large Allee effect for the unique positive constant steady state

Clearly, if $\gamma\leq (u_\ast+\rho)^2$, $E_\ast$ is locally asymptotically stable. In this subsection, one may always assume that $\gamma > (u_\ast+\rho)^2$ (i.e., $J_{11}^{*} > 0$) is satisfied and investigate the Turing instability mechanism for the considered spatial system. To investigate Turing instability, it is assumed that

are satisfied.

It is noted that $A(k^2) > 0$ for all non-negative wavenumber $k$ under assumption (4.5). Hence, Turing instability can only be attained in the case of $B(k^2) < 0$ for some nonzero wavenumber $k$, which implies that the coefficient of $k^2$ must satisfy

It should be noted that $B(k^2) > 0$ for all non-negative wavenumber $k$ when $\gamma = 0$, which implies that there is no Turing instability occurrence in this case although we can require $d > 1$, which is a classical Turing instability mechanism. So, here Turing instability may occur only when the Allee effect and spatial diffusion are present, thus, calling this mechanism Allee-effect-driven Turing instability.

The critical Turing instability condition can be found by imposing the minimum of $B(k^2)$ satisfying $\min\limits_{k^2 > 0}B(k^2) < 0$. Since (4.6), one should $\min_{k^2 > 0}B(k^2) = B(k_m^2)$ with $k_m^2 = \frac{d J_{11}^\ast+J_{22}^\ast}{2d}.$ Then $\min\limits_{k^2 > 0}B(k^2) < 0$ turns into

In fact, the parameter space defined by conditions (4.5)–(4.7) is called Turing instability space, where one can choose an arbitrary system parameter, including the parameter $\gamma$, as a bifurcation parameter to obtain Turing instability. To obtain a explicit condition revealing Turing instability, in what follows, we chose the diffusion coefficient $d$ as a bifurcation parameter to derive it. Since

one may find the following two threshold values

with $\Delta = (2J_{12}^\ast J_{21}^\ast-J_{11}^\ast J_{22}^\ast)^2-J_{11}^{\ast 2} J_{22}^{\ast 2} > 0$. Moreover, one must $d J_{11}^\ast+J_{22}^\ast = 0$ when $d = -\frac{J_{22}^\ast}{J_{11}^\ast}\triangleq d_\ast$. Furthermore, one gets $0 < \chi_- < d_\ast < \chi_+$. It is easy to check that $\min\limits_{k^2 > 0}B(k^2) < 0$ and $d J_{11}^\ast+J_{22}^\ast > 0$ hold simultaneously when $d > \chi_+\triangleq d_c$, which means $E_\ast$ is spatially unstable and Turing instability occurs in this case. Moreover, we have $\min\limits_{k^2 > 0}B(k^2) > 0$ when $d < \chi_+$, which means that $E_\ast$ is stable in this case. With these facts in hand, one may immediately have the following:

Theorem 4.1 (Turing instability induced by Allee effect). Assume that the condition (4.5) and $\gamma > (u_\ast+\rho)^2$ are satisfied. Then, system (4.1) undergoes Turing instability and spatiotemporal forms at $E_\ast$ when $d > d_c: = \chi_+$. Furthermore, Turing bifurcation occurs at $d = d_c$.

In Figure 8, we have plotted the critical $d_c$ values against the fear effect $\alpha$ with different Allee effect threshold values. We have also plotted the temporal Hopf bifurcation curves for the two different Allee effect threshold values, and these two kinds of curves divide the plane into four regions in each case: Stable region (SR), Hopf region (HR), Turing region (TR), and Turing-Hopf instability region (THR). Interestingly, from this figure, we find the $d_c$ and $\alpha^{[HB]}$ values decrease as the Allee effect factor increases, which means the Allee effect can affect Turing instability region or THR region. Next, we shall investigate pattern selection in the TR region and pattern formation in HR and THR regions, respectively.

Another way to obtain Turing bifurcation is through the steady-state solution bifurcation theory ^[40,41]. In what follows, we shall choose the diffusion ratio $d$ as a bifurcation parameter to investigate steady-state bifurcation near the positive steady state $E_\ast$. For this purpose, we rewrite $A(k^2)$ and $B(k^2)$ as $A(k^2, d)$ and $B(k^2, d)$, respectively. In view of the theory frame in Section three ^[41], the sufficient condition for the occurrence of the steady state bifurcation is that there exists a non-zero wave number $k > 0$ and $d_S > 0$ such that

and

In view of the discussion above, we immediately have:

Theorem 4.2. Suppose that (4.5) is satisfied. Furthermore, suppose that

Then, $(d_k, u_\ast, v_\ast)$ is a bifurcation point where a smooth curve $\Gamma$ of the non-constant positive steady-state solutions bifurcates from the line of positive constant steady-state solution $(u_\ast, v_\ast)$, and $\Gamma$ is contained in a connected component $\mathcal{C}$ of the set of non-zero steady-state solutions of (4.1) in $\mathbb{R}_+\times S_+$. Moreover, either $\mathcal{C}$ is unbounded in $\mathbb{R}_+\times S_+$, or $\mathcal{C}\bigcap\, (\partial\mathbb{R}_+\times S_+)\neq \emptyset$, or $\mathcal{C}$ contains another bifurcation point $(d_j, u_\ast, v_\ast)$ with $d_j\neq d_k$. More precisely, $\Gamma$ is locally a curve near $(d_k, u_{*}, v_{*})$ in the form of $\Gamma = \{(d(s), u(s), v(s)):s\in(-\epsilon, \epsilon)\}$, where $u(s) = u_{*}+sa_k\cos\frac{m\pi x}{L}\cos\frac{n\pi y}{L}+s\psi_1(s), v(s) = v_\ast+sb_k\cos\frac{m\pi x}{L}\cos\frac{n\pi y}{L}+s\psi_2(s)$ for $s\in(-\epsilon, \epsilon)$, and $d:\, (-\epsilon, \epsilon)\rightarrow \mathbb{R}$, $(\psi_1, \psi_2):\, (-\epsilon, \epsilon)\rightarrow \mathcal{Z}$ are $C^1$ functions such that $d(0) = d_k, \psi_1(0) = 0, \psi_2(0) = 0$. Here $\mathcal{Z} = \Big{\{}(u, v)\in S_+\Big{|}\int_\Omega uu_{mn}+vv_{mn}dxdy = 0\Big{\}}$ with $u_{mn} = a_k\cos\frac{m\pi x}{L}\cos\frac{n\pi y}{L}, v_{mn} = b_k\cos\frac{m\pi x}{L}\cos\frac{n\pi y}{L}$, and $a_k, b_k$ satisfy $L_k(d_k)(a_k, b_k)^T = (0, 0)^T$.

Proof. From $d^S = d_k = \frac{\mathrm{det}(J_{E_\ast})-J_{22}^{*}k^2}{k^2(J_{11}^{*}-k^2)}$, we immediately have $B(k^2, d^S) = B(k^2, d_k) = 0$, but $B(j^2, d_k)\neq 0$ for $j\neq k$. Clearly, $A(k^2, d_k) = (1+d_k)k^2-(J_{11}^{*}+J_{22}^{*}) > 0$ under the assumption (4.5) for all wavenumbers $k$. Furthermore, $\frac{\partial B(k^2, d)}{\partial d}\Big{|}_{d = d^S} = k^2(k^2-J_{11}^{*})\neq 0$. Therefore, the conditions (4.9) and (4.10) are satisfied, which ends the proof of the theorem.

In view of ^[41], the sufficient condition for the occurrence of Hopf bifurcation is that there exists a non-zero wavenumber $k > 0$ and $d_H > 0$ such that

and

For any Hopf bifurcation point $d^H$, $\alpha(\lambda)\pm \omega(\lambda)$ are the eigenvalues of $L_k(d)$, thus, (4.12) becomes

On the other hand, from $A(k^2, d^H) = 0$, one may have

when $k\in(0, \sqrt{J_{11}^\ast+J_{22}^\ast})$ with $J_{11}^\ast+J_{22}^\ast > 0,$ which requires that $J_{11}^{*} > 0.$ Clearly, $A(k^2, d^H) = 0$ and $A(j^2, d^H)\neq 0$ for $j\neq k$. Now, we only need to verify whether $B(j^2, d^H)\neq 0$ for $k\in(0, \sqrt{J_{11}^\ast+J_{22}^\ast})$, and in particular, $B(k^2, d^H) > 0$. Since

and $\Delta_k = 4J_{11}^{*^2}-4(J_{11}^{*^2}+J_{12}^{*}J_{21}^{*}) = -4J_{12}^{*}J_{21}^{*} > 0$, it follows that $B(k^2, d^H) > 0$ for

Suppose that

then

Set the quadratic function $g(z) = J_{11}^{*}J_{22}^{*}-J_{12}^{*}J_{21}^{*}-\frac{L^2(J_{11}^\ast+J_{22}^\ast)}{\pi^2}J_{11}^{*}z+\Big{(}\frac{J_{11}^{*}+J_{22}^{*}}{k_+^2}-1\Big{)}z^2$, then it is positive for $j\in\mathbb{R}$ if

We summarize the discussion above as following.

Theorem 4.3. Suppose that $J_{11}^{*}+J_{22}^{*} > 0$, $k^+ < \sqrt{J_{11}^{*}+J_{22}^{*}}$ and (4.17) holds. Then, for any $k\in(0, \sqrt{J_{11}^{*}+J_{22}^{*}})\cap(k^-, k^+)$ and $k^{-}$, $k^+$ are defined as in (4.15), the R-D system (4.1) undergoes a Hopf bifurcation at $d = d^H = \frac{J_{11}^{*}+J_{22}^{*}-k^2}{k^2}$ and the bifurcating periodic solutions from $d = d^H$ are spatially inhomogeneous.

4.2.   Turing pattern selection in two-dimensional space domains

To interpret pattern selection among spot pattern, stripe pattern, and the mixture of them, we need to derive amplitude equations near the onset $d = d_c$. Following the method employed in ^[42], we obtain the amplitude equations as follows:

where

with $I_1 = -(u_{11}+u_{00})(\varphi c_{20}+c_{11})-(\varphi c_{11}+c_{02})(v_{11}+v_{00})-\frac{3}{2}c_{12}\varphi-\frac{1}{2}c_{30}\varphi^3-\frac{1}{2}c_{03}, \; I_2 = -(u_\star+u_{00})(\varphi c_{20}+c_{11})-(\varphi c_{11}+c_{02})(v_\star+v_{00})-3c_{12}\varphi-c_{30}\varphi^3-c_{03}, \; J_1 = -(u_{11}+u_{00})d_{11}-(\varphi d_{11}+d_{02})(v_{11}+v_{00})-\frac{3}{2}d_{12}\varphi, \; J_2 = -(u_{\star}+u_{00})d_{11}-(\varphi d_{11}+d_{02})(v_{\star}+v_{00})-3d_{12}\varphi$. Here, $\varphi, \psi, u_{00}, v_{00}, u_{11}, v_{11}, u_\star$ and $v_\star$ possess the similar expressions as those in Section 4.3.1 ^[42]. Denote

Clearly, $\mu_1 < \mu_2 < \mu_3 < \mu_4$ under the assumption of $g_1 > 0$ and $g_2 > 0$. Then, we have some results about the stability of the hexagonal pattern, the stripe pattern, and the mixture of them according to the range of $\mu$ in Theorem 3 ^[42].

4.3.   Numerical Simulations of pattern formation in Turing and non-Turing regions

In this subsection, we performed some numerical simulations to show the pattern formation over two dimensional spatial domain according to the obtained theoretical result (see Theorem 4.1), and indicated some interesting pattern solutions in either the Hopf instability region or Turing-Hopf region besides Turing region (see Figure 8). We solved the reaction-diffusion system (4.1) numerically by using a five-point finite difference scheme for the Laplacian operator and forward Euler scheme for the temporal part with appropriate care at the boundary to accommodate Neumann boundary conditions. All the programs for solving the spatiotemporal system (4.1) were performed by MATLAB R2016a. In addition, to ensure the stability of numerical method, the choice of time step size $\Delta t$, the space step sizes $\Delta x$ and $\Delta y$ and the diffusion coefficient $d$ should satisfy the Courant-Friedrichs-Lewy stability criterion, i.e., $\max\{d, 1\}\Delta t\bigg{(}\frac{1}{(\Delta x)^2}+\frac{1}{(\Delta y)^2}\bigg{)}\leq \frac{1}{2}$. In addition, we have checked that the numerical simulation results for some other smaller values of the time and space step sizes and they remain the same qualitatively.

In Figure 9, we observed pattern selections among the hexagonal pattern, stripe pattern and a mixture of hexagons and stripes for different diffusion coefficients with a set of fixed system parameters, where $\alpha = 0.32, \beta = 0.048, \gamma = 0.0826, \rho = 0.135, \eta = 0.8, \delta = 0.012$ and $\sigma = 0.128$. Here, the initial condition is chosen as a random perturbation around the spatially homogeneous steady state $E_\ast$. Under this set of parameters, we obtained that the one of the Hopf bifurcation threshold values is $\alpha = \alpha^{[HB]} = 0.3254$ and Turing bifurcation threshold value is $d = d_c = 37.9638$. In view of (4.19), we get the coefficients of the amplitude equations (4.18) are $h_0 = -25.0863 < 0$ and $g_2 = 262.3196 > g_1 = 7.7015 > 0$. Moreover, we have $\mu_1 = -0.2955$, $\mu_2 = 0$, $\mu_3 = 0.0748$ and $\mu_4 = 2.6959$, which will not change with the variation of diffusion coefficient $d$. So, when $d = 40.7$, we obtained $\mu = 0.0721\in(\mu_2, \mu_3)$, which means the hexagonal pattern prevails over the whole domain as time goes forward (see Figure 9(ⅰ)). When $d = 150$, we obtained $\mu = 2.9511\in(\mu_4, \infty)$, which implies that the stripe pattern will be formed as time evolves (see Figure 9(ⅲ)). However, when $d = 60$, we obtained $\mu = 0.5805\in(\mu_3, \mu_4)$. This manifests that the mixture of hexagons and stripes forms over the two dimensional spatial domain (see Figure 9(ⅱ)).

In Figure 10, we took another set of system parameters as $\alpha = 0.4, \beta = 0.042, \gamma = 0.1, \rho = 0.12, \eta = 0.86, \delta = 0.06$ and $\sigma = 0.18$ with variational diffusion coefficients $d = 55, 58, 62$ in the TR. The initial condition was chosen as a random perturbation around the coexistence equilibrium. One of the critical Hopf bifurcation values for the temporal system and the critical Turing bifurcation value for the spatial system are $\alpha^{[HB]} = 0.4334$ and $d_c = 53.6568$, respectively. We have plotted the spatiotemporal bifurcation diagram for Turing bifurcation and Hopf bifurcation (see Figure 10(ⅰ)). A direct calculation gives $g_1 = -425.7471 < 0$ and $g_2 = -853.6242 < 0$, which suggests the amplitude equations (4.18) undergoes sub-critical bifurcation. Interestingly, we observed only stable hexagonal pattern solution prevails over the spatial domain as time evolved, although the diffusion coefficients are different (see Figure 10(ⅱ)–(ⅳ)).

Apart from the stationary pattern solution, we are now interested in pattern dynamics in HR and THR regions, which are non-Turing regions. In Figure 11, the system parameters remain consistent with those depicted in Figure 10, with the exception of alterations to the fear factor $\alpha$ and the diffusion coefficient $d$. Under this specific parameter configuration, it is observed that one of the Hopf-bifurcation threshold values is $\alpha = \alpha{^{[HB]}} = 0.4334$. The initial choice is made to set $\alpha = 0.5$, a value situated within the temporal Hopf instability region. Subsequently, the diffusion coefficient $d$ is selected to be $40$, which falls below the Turing bifurcation threshold value ($d_c = 42.6939$). This choice therefore positions the system parameters within the pure Hopf instability region, as indicated in Figure 10(ⅰ). Within Figure 11, the spatiotemporal distribution of the prey species $u$ at various time points is demonstrated. It is noteworthy that the initial condition employed in ^[43] is utilized, elucidating that only a small number of populations from both species are confined to an elliptic domain.

Here, $u_0 = 0.8$ and $v_0 = 0.2$, $x_1 = 153.5, y_1 = 145, x_2 = y_2 = 150, \Delta_{11}^2 = \Delta_{12}^2 = 12.5, \Delta_{21}^2 = 5, \Delta_{22}^2 = 10$. We first observe some circular rings with different densities of the species on it (see Figure 11(ⅰ)) at $t = 200$. As time goes forward, the expanding circular rings hit the domain boundary and break into the stripes (see Figure 11(ⅱ), (ⅲ)). With the further advancement of time, the circular rings gradually lose stability and break into some spots and stripes (see Figure 11(ⅳ)–(ⅵ)), and eventually settle down to coldspot pattern (see Figure 11(ⅶ)).

In Figure 12, we observe a similar spatiotemporal distribution of prey species except the difference of instability mechanism of the circular rings pattern, where the system parameters are chosen the same as those in Figure 11 besides the diffusion coefficient $d = 45$, which above the Turing bifurcation threshold value $d = d_c = 42.6939$. It is easy to check that here the system parameters are located in THR (see Figure 10(ⅰ)) and the initial conditions are in (4.21). After experiencing the rupture of the circular rings, the final pattern transition state is the mixture pattern of spots and stripes (see Figure 12(ⅶ)). It should be noted that we don't observe spatiotemporal chaotic pattern in THR for the certain system parameters and the special choice of initial condition.

5.   Concluding remarks

The present article is dedicated to the derivation of comprehensive theoretical outcomes for various types of bifurcations, followed by their substantiation through numerical simulations. Each numerical simulation is underpinned by the selection of distinct system parameters, treating them as bifurcation parameters, thereby showcasing diverse dynamical responses. A particular focus is placed on the intricate influence of the Allee effect on both temporal and spatiotemporal systems.

For the temporal system, a pivotal finding emerges, indicating that the Allee effect can instigate both codimension-one bifurcations (transcritical, saddle-node and Hopf) and codimension-two bifurcations (cusp and Bogdanov-Takens). This revelation imbues the system with a more intricate dynamical behavior, signifying an enrichment of the dynamics inherent to the model system. Notably, the induced Hopf bifurcations tend to be subcritical, aligning qualitatively with previous observations ^[9]. Additionally, it established that the temporal system can exhibit bistability between the trivial equilibrium $E_0$ and the coexistence equilibrium $E^*$ (see Remark 3.1).

For the spatiotemporal model, we investigated the system's dynamics and found that the Allee effect is a key mechanism for inducing Turing instability, and the small Allee effect usually plays a stabilization role (see Theorem 4.1). Furthermore, we have investigated its pattern dynamics and derived amplitude equation near the onset $d = d_c$. We have studied pattern transition among a hexagonal pattern, stripe pattern and a mixture of them and established an analytical formula to predict it, and we have preformed numerical simulations to illustrate the effectiveness of the theoretical results (see Figure 9). In addition, we found some interesting circular ring patterns in the HR and THR for a special choice of initial condition at the initial moment, which eventually transitioned to coldpattern or the mixture pattern consisting of coldspots and stripes (see Figures 11 and 12). This was different from the observed spatiotemporal chaotic pattern in a spatial slow-fast prey-predator system with the same initial condition ^[44].

It was imperative to emphasize that the outcomes of this theoretical investigation hold relevance as long as the fear factor $\alpha$ and the Allee parameter $\gamma$ remained within their practical range of values. Beyond this range, the system experienced permanent collapse due to negative growth in the prey species. Future enhancements of the model could involve the introduction of multiplicative or double Allee effects, offering potential for exciting developments in this challenging field of research. These future investigations are anticipated to introduce a fresh dimension to the domain of mathematical ecology.

Use of AI tools declaration

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

Acknowledgments

The authors are grateful to the anonymous reviewers for their valuable comments and suggestions. The first author gratefully acknowledges the financial support from Anhui Province University Excellent Talents Funding Project (Grant No. gxbjZD2021075) and Anhui Province University Research Project (Grant No. 2023AH051549). The last author is supported by Peak Discipline Construct Project of Zhejiang University of Science and Technology.

Conflict of interest

The authors declare there is no conflict of interest.

Appendix.   The proof of Theorem 3.12: a detailed derivation of generic B-T norm form for the temporal system (3.2)

Step-1: (Derivation of non-degeneracy condition)

Now, first of all it will be shown that the coexistence equilibrium point $E_{*} = (u_{*}, v_{*})$ is a cusp of codimension two. For this purpose, one may choose the following transformation $A_{1} = u-u_{*}, A_{2} = v-v_{*}$, which shifts the coexistence equilibrium point $E_{*}$ into the origin. By power series expansion of the modified system into second order terms, one obtains

where $A = (A_{1}, A_{2})$ and the coefficients $p_{ij} = c_{ij}, q_{ij} = d_{ij}$ for $i, j = 0, 1, 2, 3, \ldots$ are all provided in the Subsection $3.3.2$. Therefore, one must have

Now, let us consider the transformation $u_{1}^{'} = A_{1}, \; \; u_{2}^{'} = p_{10}A_{1}+p_{01}A_{2}.$ Using this transformation, (3.2) may be represented as

where $u^{'} = (u_{1}^{'}, u_{2}^{'})$ and the coefficients $r_{ij}, s_{ij}$ with $i, j = 0, 1, 2, 3, \ldots$ are given by

Step-2: (Introducing a $c^{\infty}$ transformation)

Now, there must be a $c^{\infty}$ invertible transformation given by

This transformation reduced (A.3) to the following one:

where $v^{'} = (v_{1}^{'}, v_{2}^{'})$ and the coefficients $\mu_{1}, \mu_{2}$ are given by

Step-3: (Introduction of a new transformation)

Now, let us introduce a new transformation in a very small neighborhood of $v^{'} = 0$ with the form

By use of this transformation, (A.5) may be represented as

where $\varkappa = (\varkappa_{1}, \varkappa_{2})$. Therefore, if $\mu_{1}\mu_{2} \neq 0$ (non-degeneracy condition) then a cusp bifurcation of codimension two occurs around the unique coexistence equilibrium point $E_{*} = (u_{*}, v_{*})$ ^[45].

Step-4: (Deduction of generic normal form)

Now, assume that the non-degeneracy condition has been well satisfied. So, one may try to calculate the generic normal form of the Bogdanov-Takens bifurcation according to the frame in ^[38]. First, we utilize the transformation $\bar{u}_1 = u-u_\ast$ and $\bar{v}_1 = v-v_\ast$ to shift the equilibrium point $E_\ast$ of (3.2) to the origin, and then the system (3.2) becomes

Let us perturb the bifurcating parameters $(\sigma, \gamma)$ in a small neighborhood of $\epsilon_{i}$ in such a way that $(\sigma, \gamma) = (\sigma^{[BT]}+\epsilon_{1}, \gamma^{[BT]}+\epsilon_{2})$. Then, we obtain

By expanding the Taylor's series of (A.8) at $(\bar u_1, \bar v_1) = (0, 0)$ up to the terms of second order, we get the following system

where the expression for $m_{ij}, n_{ij}$; $i, j = 0, 1, 2, 3, \ldots$ and $G_{10}, H_{10}$ are defined by

where the terms $\bar{H}_{1}(\bar{u}_{1}, \bar{v}_{1})$ and $\bar{H}_{2}(\bar{u}_{1}, \bar{v}_{1})$ are the power series of $\bar{u}_{1}$ and $\bar{v}_{1}$ and is the form for each of them $\bar{u}_{1}^{i}\bar{v}_{1}^{j}$; $i+j \geq 3$.

Step-5: (Introduction of a affine transformation)

Now, let us introduce a new affine transformation given by

By using this affine transformation, (A.9) reduced to the following one:

where, the coefficients $c_{ij}, d_{ij}$; $i, j = 0, 1, 2, 3, \ldots$, $a_{1}^{'}, b_{1}^{'}$ and $H_{10}^{'}$ are given by

while the remaining terms $H_{1}^{'}(p_{1}, p_{2})$ and $H_{2}^{'}(p_{1}, p_{2})$ are the power series expansion terms of $p_{1}$ and $p_{2}$ of the form $p_{1}^{i}p_{2}^{j}$; $i+j \geq 3$. Furthermore, for $\epsilon^{*} = (0, 0)$, we have $H_{10}^{'}(\epsilon^{*}) = 0$, $a_{1}^{'}(\epsilon^{*}) = 0$, $b_{1}^{'}(\epsilon^{*}) = 0$.

Step-6: (Setting of new variables)

Now, let us choose a transformation as described by

Employing this transformation into (A.10), one obtains

where $I(p_{1}^{'}, p_{2}^{'})$ is the power series of $p_{1}^{'}$ and $p_{2}^{'}$ of the form ${p_{1}^{'i}}{p_{2}^{'j}}$ with $i+j \geq 3$; whereas the remaining coefficients $L_{ij}$; $i, j = 0, 1, 2, \ldots$ are given by

Also, one may easily verify that

Therefore, the system of equations of (A.11) may be written as

where, the terms $\rho_{1}, \rho_{2}, \Phi_{1}^{'}$ are all smooth functions and $\rho_{1}(0, \epsilon^{*}) = L_{00}(\epsilon^{*}) = 0$, $\frac{\partial \rho_{1}}{\partial p_{1}^{'}}\bigg{\rvert}_{(0, \; \epsilon^{*})} = L_{10}(\epsilon^{*}) = 0$, $\frac{1}{2}\frac{\partial^{2} \rho_{1}}{\partial {p_{1}^{'}}^{2}}\bigg{\rvert}_{(0, \; \epsilon^{*})} = L_{20}(\epsilon^{*}) = d_{20}(\epsilon^{*})\neq 0$, $\rho_{2}(0, \epsilon^{*}) = L_{01}(\epsilon^{*}) = 0$. Also, one obtains

Since we have $\rho_{2}(0, \epsilon^{*}) = 0$ and $\frac{\partial \rho_{2}}{\partial p_{1}^{'}}\bigg{\rvert}_{(0, \; \epsilon^{*})} \neq 0$, by employing implicit function theorem, there exists a $c^{\infty}$ function $p_{1}^{'} = \phi_{1}^{'}(\epsilon)$ defined in a very small neighborhood $N(\epsilon^{*})$ of $\epsilon = \epsilon^{*}$ in such a way that $\phi_{1}^{'}(\epsilon^{*}) = 0$, $\rho_{2}(\phi_1^{'}(\epsilon), \epsilon) = 0$, for any choices of $\epsilon \in N(\epsilon^{*})$.

Step-7: (Parameter dependent shift)

Now, to annihilate the $p_{2}^{'}$-term from the second equation of (A.12), the following parameter dependent shift transformation has been taken into account

By using this shift transformation, (A.12) may be represented as

where $E = (E_{1}, E_{2})$, and

The coefficient of $E_{2}$ of the system of equations of (A.13) may be written as

From the first condition of $\mathrm{BT. 2}$, we conclude that $e_{01}(0, \epsilon^{*}) = L_{01}(\epsilon^{*}) = 0$, $\frac{\partial e_{01}}{\partial \phi_{1}^{'}}\bigg{|}_{(0, \epsilon^{*})} = L_{11}(\epsilon^{*}) = d_{11}(\epsilon^{*})+2c_{20}(\epsilon^{*}) = \theta \neq 0$, which indicates that there exists a smooth function $\phi_{1}^{'} = \phi_{1}^{'}(\epsilon)$ in such a way that $\phi_{1}^{'}(\epsilon^{*}) = 0$, $e_{01}(\phi_{1}^{'}(\epsilon), \epsilon) = \delta_{2}(0, \epsilon) = 0$ for any choices of $\epsilon \in N(\epsilon^{*})$. So, to annihilate the $E_{2}$-term from the second equation of (A.13), we consider $e_{01} = 0$ and one obtains

Therefore, the system of equations of (A.13) may be written as

which may be represented as

where $E = (E_{1}, E_{2})$.

Step-8: (Time reparametrization)

Furthermore, here $e_{00}(\epsilon^{*}) = e_{10}(\epsilon^{*}) = 0$. Introducing a new time scaling defined by $dt = (1+\omega E_{1})d\tau$, and by using it (A.14) becomes

One may consider a new transformation of the form

By use of this transformation, (A.15) may be written as

where, $H(P, Q)$ are the power series expansion of $P$ and $Q$ of the form $P^{i}Q^{j}$; $i+j \geq 3$. Also, the coefficients $f_{ij}$; $i, j = 0, 1, 2, 3, \ldots$ are given by

Now, to annihilate the $Q^{2}$-term, let us choose $\omega(\epsilon) = -e_{02}(\epsilon)$ and we obtain

where $\gamma_{1}(\epsilon) = e_{00}(\epsilon)$, $\gamma_{2}(\epsilon) = e_{10}(\epsilon)-2e_{00}(\epsilon)e_{02}(\epsilon)$, $\mathcal{D}(\epsilon) = e_{20}(\epsilon)-2e_{10}(\epsilon)e_{02}(\epsilon)+e_{00}(\epsilon)e_{02}^{2}(\epsilon)$, $\mathcal{E}(\epsilon) = e_{11}(\epsilon)$. Also, one may verify that, $\gamma_{1}(\epsilon^{*}) = 0$, $\gamma_{2}(\epsilon^{*}) = 0$, and $\mathcal{E}(\epsilon^{*}) = L_{11}(\epsilon^{*}) = \theta \neq 0$ because of ${{\bf{BT}}. \; {\bf{2}}}$. is satisfied.

Step-9: (Time rescaling)

Since the second condition of ${{\bf{BT}}. \; {\bf{2}}}$ is satisfied, $\mathcal{D}(\epsilon^\ast) = e_{20}(\epsilon^\ast) = L_{20}(\epsilon^\ast) = d_{20}(\epsilon^\ast)\neq 0$. Thus we can choose a $\delta^\ast-$neighborhood such that $\mathcal{D}(\epsilon)\neq 0$ when $\epsilon\in O(\epsilon_\ast, \delta^\ast)$. Let us introduce a new time scaling that is of the form $t = \bigg{|}\frac{\mathcal{E}(\epsilon)}{\mathcal{D}(\epsilon)}\bigg{|}\tau.$ Again, choose $\xi_{1} = \frac{\mathcal{D}(\epsilon)}{\mathcal{E}^{2}(\epsilon)}$, $\xi_{2} = s\bigg{(}\frac{\mathcal{D}^{2}(\epsilon)}{\mathcal{E}^{3}(\epsilon)}\bigg{)}Q$, where $s = sign\bigg{(}\frac{\mathcal{E}(\epsilon)}{\mathcal{D}(\epsilon)}\bigg{)} = sign \bigg{(}\frac{\mathcal{E}(\epsilon^{*})}{\mathcal{D}(\epsilon^{*})}\bigg{)} = sign\bigg{(}\frac{\theta}{d_{20}(\epsilon^{*})}\bigg{)} = \pm 1$.

Therefore, the reduced model system (A.16) takes the form

where $\xi = (\xi_{1}, \xi_{2})$,

In a very small neighborhood of the origin, (A.17) is topologically equivalent to the following one

which is the generic normal form of Bogdanov-Takens bifurcation. Combining BT. 3, we can define an invertible change of parameters near the origin, which is also equivalent to the regularity of the map $\epsilon\rightarrow (\vartheta_1(\epsilon), \vartheta_2(\epsilon))$ at $\epsilon = \epsilon^\ast = (0, 0)$. Therefore, one may conclude that (3.2) exhibits a Bogdanov-Takens bifurcation around the coexistence equilibrium point $E_{*} = (u_{*}, v_{*})$ whenever the system parameters $(\sigma, \gamma)$ attains the critical value $(\sigma^{[BT]}, \gamma^{[BT]})$, which implies the proof of the theorem as complete.