Analysis of a prey-predator system incorporating the additive Allee effect and intraspecific cooperation

1.   Introduction

From single species to different species interacting in the same habitat, there is a growing interest in the study of population evolution. Population dynamics is one of the most discussed topics among scholars. The study of the predator-prey system focuses on qualitative aspects, such as equilibrium, local and global stability, and the presence of limit cycles. The predator-prey system, originally proposed by Alfred J. Lotka and Vito Volterra in 1925, holds a prominent position in the field of ecological dynamics. The underlying foundation of this model rests upon the premise of linear behaviors. Since its inception, a significant corpus of theoretical research has emerged. This system has been extensively employed and implemented in various fields, including chemical reactions ^[1], plasma physics ^[2], hydrodynamics ^[3], as well as social science and economics ^[4]. Therefore, population models, especially predator-prey models, are extremely important in both theoretical research and other fields of application. From a theoretical research perspective, studying them is also highly intriguing. In fact, the predator-prey model has consistently been a popular topic in the field of mathematical biology research, especially in two-dimensional space. This refers to a situation in which there is one predator population and one prey population.

Stability is a paramount concern in qualitative research on ecosystems, and scholars have conducted numerous investigations in this field, resulting in a wealth of comprehensive findings. In certain ecological systems, prey may exhibit fear responses toward predator, thereby altering their behavior and making it more challenging for predators to successfully capture their prey. For example, playing the sounds of known predators in the wild reduces the spawning rate of songbirds by forty percent. Some researchers have also investigated various models of fear effect in conjunction with different functional responses ^[5,6,7].

The Allee effect manifests when population densities are at a low level, thereby posing challenges for individuals in terms of reproduction and survival. The Allee effect is a well-documented biological phenomenon that provides insight into the correlation between the population's growth rate and its size or density. The proposal of this concept can be traced back to W.C. Allee, as documented in the reference ^[8]. According to its density dependence at low density, the Allee effect can be divided into two types: strong Allee effect and weak Allee effect. For example, the island fox (Urocyon littoralis) and the gypsy moth (Lymantria dispar), an invasive species, will suffer from a strong Allee effect ^[9,10], while another invasive species, the smooth cordgrass (Spartina alterniflora), will show a weak Allee effect ^[11,12]. In recent years, there has been a notable focus among researchers on the investigation of the Allee effect. The Allee effect significantly impacts the dynamic behavior of the system. For instance, the Allee effect has the potential to modify the local stability of the system singularity or induce system instability ^[13,14]. Additionally, it can prolong the time required for the system to reach a stable equilibrium point ^[15,16].

Mathematically, the Allee effect is commonly characterized by modifying the growth function, usually achieved by introducing a multiplication factor. In the study conducted by the author ^[17], the Lotka-Volterra predator-prey system was examined, with particular attention given to the incorporation of the Allee effect:

where $x$ and $y$ represent the densities of the prey population and the predator population, respectively. The term $\alpha(x) = x/(\beta+x)$ represents the Allee effect. The system always has an unstable trivial equilibrium point and an unstable semi-trivial equilibrium point, which is consistent with the results of the system without the Allee effect($\beta = 0$). In classical models, a coexisting equilibrium point is ultimately present and globally asymptotically stable. However, the Allee effect disrupts this by only allowing a positive equilibrium to exist and be locally asymptotically stable when $r-a\beta > 0$. The Allee effect, while not altering the positive equilibrium point of the system, exerts two distinct influences on the dynamics of the system. Firstly, the prolongation of the time required for the system to achieve equilibrium occurs. Secondly, it leads to a decrease in the population density of both populations. If we consider the fact that species higher up in food chains are more likely to face extinction, Guan et al. in ^[18] focused their study on the Allee effect on the predator population, rather than the prey population. The proposed system is analogous to the one described in ^[17], as it considers the identical function $\alpha$:

Their research shows that in cases where $a > b$, the system exhibits persistence, indicating that both the trivial and semi-trivial equilibrium points are characterized by instability. Regarding the local stability of the positive equilibrium point, the authors obtained identical findings to those reported in the study by ^[17]. The Bendixon-Dulac criterion was employed to confirm the absence of a limit cycle, thus proving that the only positive equilibrium point is globally asymptotically stable. At the end of the study, the authors' numerical simulations demonstrated that the presence of the Allee effect did not have any impact on the final density of both predator and prey populations. This outcome exhibits disparity when compared to the findings reported in ^[17]. Afterward, Chen et al. in ^[19] chose a more realistic growth model for prey birth rate, known as the Beverton-Holt function with saturation effect, instead of the logical model that was previously used. Similar to previous systems, the conclusions regarding the stability of equilibrium points are the same.

Wang et al. in ^[20] proposed a predator-prey model with the Allee effect that is more comprehensive in nature. They are considering the following system:

The function $f(x)$ satisfies the following conditions: $f(x)\in\mathcal{C}^{1}([0, \infty), f(b) = f(1) = 0$, with $0 < b < 1$, $f(x)$ is positive for $b < x < 1$ and negative in other case, and there exists $\bar{\lambda}\in(b, 1)$ such that $f'(x) > 0$ in $[b, \bar{\lambda})$ and $f'(x) < 0$ in $(\bar{\lambda}, 1]$. Their assumption of $f$ indicates that the prey population exhibits a strong Allee effect. Of course, there are some models that cannot be described by the generic model mentioned above. For instance, ^[21] studied a model of prey population affected by the Allee effect, which is based on the Leslie-Gower model:

In addition, another common way to introduce the Allee effect in the model is to include an additional term in the net population growth rate. In ^[22], Dennis first proposed the single species model with an additive Allee effect:

The term $m/(x+a)$ represents the additive Allee effect, where $m$ and $a$ reflect the degree of Allee effect. Later, many scholars also adopted this modeling idea in their studies. For example, the additive Allee effects were considered in the predator-prey model ^[23] and the spatiotemporal prey-predator model ^[24], two-patch model ^[25], and the Lotka-Volterra commensalism model ^[26]. The future of the Allee effect awaits us not only at the empirical level but also in theoretical research. There are still many questions that need to be addressed through models.

Intraspecific cooperation is widespread in nature ^[27]. For example, Emperor penguins engage in a behavior known as huddling, where they gather closely together and regularly rotate positions from the outer edges to the inner core. This collective action serves to regulate their body temperature and offers them protection against the extreme climatic conditions prevalent in Antarctica. Wolves exhibit cooperative behavior by coordinating their movements to acquire sustenance and various other essential resources. In ^[28], Jacobs proposed to illustrate the dynamics of intraspecific cooperation:

The term $\frac{\xi\nu x}{x+\nu}$ indicates that intraspecific cooperation exhibits a Holling-type saturation, where $\xi > 0$ quantifies the effectiveness of intraspecific cooperation. The authors state that the model exhibits rich dynamical behaviors and has numerous potential applications. There are various environmental disturbances in nature, and neglecting the unpredictability of the environment may be one of the main reasons why many species are extinct or on the brink of extinction ^[29]. In ^[30], Zhang et al. developed a model that incorporates intraspecific cooperation and stochastic perturbation:

They provided sufficient criteria for species extinction and permanence and analyzed the existence of UESD. They also derived an explicit form of the density function of UESD under certain conditions. It has also been shown that perturbations in the growth rate, such as white noise, can lead to species extinction. Additionally, the presence of telephone noise may cause species to disappear or persist.

In general, intraspecific cooperation can improve survival by enhancing food acquisition efficiency (evening bats), increasing prey capture effectiveness (Harris hawks, lions, spotted hyena), providing nourishment and protection to injured individuals (African wild dogs), and collectively defending attacked individuals (White-winged choughs, Florida scrub jays, banded mongoose, Geoffroys marmoset). It can also increase reproduction by extending the reproductive lifespan of the breeders (Chimpanzees, pied kingfisher, splendid fairy wren, dwarf mongoose) or by reducing the time interval between reproductive events (prairie voles, house mice) and increasing the number of litters produced (Black-backed jackals). It can further increase the survival rate of offspring by feeding (Arabian babblers, meerkats) or protecting the offspring (African wild dogs).

We believe that intraspecific cooperation plays an important role at higher levels of the food chain. This is based on the fact that species at higher levels are more vulnerable to extinction compared to species at lower levels, particularly when the higher-level species are specialists rather than omnivorous. Therefore, we will consider a predator-prey system in which the obligate predator exhibits intraspecific cooperation. Unlike ^[28] and ^[30], which directly incorporate a function related to population density into the model to describe intraspecific cooperation, we will represent intraspecific cooperation of the predator by modifying predator mortality as a decreasing function of population size. For example, emperor penguins huddle together in low temperatures and harsh winds to ensure they do not freeze to death. This form of intraspecific cooperation directly reduces the mortality rate among individuals in harsh environments. So, we investigated the following system:

where $x$ and $y$ represent the density of prey and predator populations, respectively, and all parameters are positive constants. The term $m/(x+\alpha)$ represents the additive Allee effect, while the term $e/(y+\beta)$ represents intraspecific cooperation. The meaning of the parameters is explained in Table 1. The rest of the paper is organized as follows: In Sec 2, the existence and stability of the equilibria are discussed. In Sec 3, the occurrence of Hopf bifurcation and the limit cycle is proved. In Section 4, the analysis focuses on the influence of the Allee effect and intraspecific competition on the system. In Sec 5, numerical simulations are provided to verify the validity of the conclusions.

Table 1.  Meaning of all parameters.

Parameter	Meaning
$a$	The intrinsic growth rate
$b$	The coefficient of intraspecific
$c$	The predation rate
$e$	The death rate
$\varepsilon \in (0, 1)$	The conversion factor
$\alpha$, $m$	The Allee constants
$\beta$	The intraspecific cooperation strength

---

 | Show Table

DownLoad: CSV

2.   The existence and stability of equilibria

Firstly, all possible equilibria of system (1.1) are discussed. The equilibria of system (1.1) are determined by the following equations:

Obviously, Equation (2.1) always has a trivial equilibrium $E_{0}(0, 0)$. From the second equation of (2.1), we have its nonnegative roots: $y = 0$ and $x = \frac{e}{c\varepsilon}$. (1) Substituting $y = 0$ into the equation $a-bx-cy-\frac{m}{x+\alpha} = 0$ and simplifying it, we obtain the following quadratic equation

Then we consider the positive root of Eq (2.2).

(i) If $m-a\alpha < 0$, then Eq (2.2) has a unique positive root $x_{1} = \frac{a-b\alpha+\sqrt{(a+b\alpha)^{2}-4bm}}{2b}$.

(ii) If $m-a\alpha = 0$ and $b\alpha-a < 0$, then Eq (2.2) has a unique positive root $x_{2} = \frac{a-b\alpha}{b}$.

(iii) If $m-a\alpha > 0$, $b\alpha-a < 0$ and $(a+b\alpha)^{2}-4bm > 0$, then Eq. (2.2) has two positive roots

(iv) If $m-a\alpha > 0$, $b\alpha-a < 0$ and $(a+b\alpha)^{2}-4bm = 0$, i.e. $a > b\alpha$ and $m = \frac{(a+b\alpha)^{2}}{4b}$, then Eq (2.2) has only one positive root $x_{5} = \frac{a-b\alpha}{2b}$.

(2) Substituting $x^{*} = \frac{e}{c\varepsilon}$ into the equation $a-bx-cy-\frac{m}{x+\alpha} = 0$, we have $y^{*} = \frac{1}{c}(a-\frac{be}{c\varepsilon}-\frac{mc\varepsilon}{e+c\alpha\varepsilon}) > 0$ if $a > \frac{be}{c\varepsilon}+\frac{mc\varepsilon}{e+c\alpha\varepsilon}$.

The existence of all equilibria of system (1.1) is shown in Table 2.

Table 2.  Equilibria existence and stability.

Equilibria	$a$	$m$	Type
$E_{0}(0, 0)$			$m > a\alpha$ or $m=a\alpha=b\alpha^{2}$, stable node;
			$m < a\alpha$, saddle point;
			$m=a\alpha\neq b\alpha^{2}$, saddle node point.
$E_{1}(x_{1}, 0)$		$m < a\alpha$	$m < min\{a\alpha, b\alpha^{2}\}$, saddle point.
$E_{2}(x_{2}, 0)$	$a > b\alpha$	$m=a\alpha$	$a-b\alpha < eb/c\varepsilon$, stable node;
			$a-b\alpha > eb/c\varepsilon$, saddle point.
$E_{3}(x_{3}, 0)$	$a > b\alpha$	$a\alpha < m < \frac{(a+b\alpha)^{2}}{4b}$	unstable node or saddle point.
$E_{4}(x_{4}, 0)$			$a-b\alpha < eb/c\varepsilon$ or $x_{4} < \frac{e}{c\varepsilon}$, stable node;
			$a-b\alpha > 2eb/c\varepsilon$ or $x_{4} > \frac{e}{c\varepsilon}$, saddle point.
$E_{5}(x_{5}, 0)$	$a > b\alpha$	$m=\frac{(a+b\alpha)^{2}}{4b}$	$J_{11}(E_{5})=0$, nonhyperbolic equilibrium point.
$E^{*}(x^{*}, y^{*})$	$a > \frac{be}{c\varepsilon}+\frac{mc\varepsilon}{e+c\alpha\varepsilon}$		$m < b(\frac{e}{c\varepsilon}+\alpha)^{2}$, stable node;
			$m > b(\frac{e}{c\varepsilon}+\alpha)^{2}$, unstable node.

---

 | Show Table

DownLoad: CSV

Secondly, the stability of all equilibria is taken into consideration. The Jacobian matrix of the system (1.1) is

where $J_{11}=a-2bx-cy-\frac{m\alpha}{(x+\alpha)^{2}}$, $J_{12}=-cx$, $J_{21}=\frac{c\varepsilon y}{y+\beta}$ and $J_{22}=\frac{\beta(-e+c\varepsilon x)}{(y+\beta)^{2}}$.

Theorem 1. $E_{0}(0, 0)$ is a saddle if $m < a\alpha$ holds, a stable node if $m > a\alpha$ or $m = a\alpha = b\alpha^{2}$ holds, and a saddle node if $m = a\alpha\neq b\alpha^{2}$ holds.

Proof. The Jacobian matrix at $E_{0}(0, 0)$ is

whose eigenvalues are

Consequently, $E_{0}(0, 0)$ is a saddle if $m < a\alpha$ holds, and it is a stable node if $m > a\alpha$ holds.

If $m = a\alpha$, the system (1.1) can be expanded into a power series around $E_{0}$ by using the formula $\frac{1}{1+x} = \sum\limits_{i = 0}^{\infty}(-1)^{i}x^{i}$ $(|x| < 1)$ and making a time transformation $d\tau = -\frac{e}{\beta}dx$.

Based on the implicit function theorem, we know from $y+Q(x, y) = 0$ that there exists only one function $y = \varphi(x)$ such that $\varphi(0) = \varphi'(0) = 0$. We can solve for $y$ to obtain $y = \varphi(x)\equiv0$. Substituting this into the first equation, we get

Consequently, $E_{0}(0, 0)$ is a stable node if $m = a\alpha = b\alpha^{2}$ holds, and it is a saddle node if $m = a\alpha\neq b\alpha^{2}$ holds. □

Let $h_{1}(x) = \frac{m}{(x+\alpha)^{2}}-b, \, \, h_{2}(x) = -e+\varepsilon cx.$

Theorem 2. $E_{i}(x_{i}, 0)$ is a saddle if $h_{1}(x_{i})h_{2}(x_{i}) < 0$ holds, and is a node if $h_{1}(x_{i})h_{2}(x_{i}) > 0$ holds.

Proof. The Jacobian matrix at $E_{i}$ is

where

and

whose eigenvalues are $\lambda_{1} = J_{11}(E_{i})$ and $\lambda_{2} = J_{22}(E_{i})$. Hence, $E_{i}(x_{i}, 0))$ is a saddle if $J_{11}(E_{i})J_{22}(E_{i}) < 0$ holds, and it is a node if $J_{11}(E_{i})J_{22}(E_{i}) > 0$ holds.

In particular, $E_{1}(x_{1}, 0)$ is a saddle if $m < min\{a\alpha, b\alpha^{2}\}$ holds. This is due to the fact that if $m < min\{a\alpha, b\alpha^{2}\}$, we have

Now, we consider the sign of $J_{22}(E_{1})$. Let $f(x, y) = a-bx-cy-\frac{m}{x+\alpha}$, then $f(x^{*}, y^{*}) = 0$ and for any $x > 0$, $\frac{\partial f(x, y)}{\partial y} = -c < 0$. Consequently,

then

If $m < min\{a\alpha, b\alpha^{2}\}$ holds, then $b-\frac{m}{(x_{1}+\alpha)(x^{*}+\alpha)} > b-\frac{m}{\alpha^{2}} > 0$, therefore, $x_{1}-x^{*} > 0$. So

For $E_{2}(x_{2}, 0)$,

So, $E_{2}$ is a stable node if $c\varepsilon(a-b\alpha) < eb$ holds, and it is a saddle if $c\varepsilon(a-b\alpha) > eb$ holds. Combining the conditions for the existence of $E_{2}$, we obtain that when $m = a\alpha$, $E_{2}$ is a stable node if $0 < a-b\alpha < \frac{eb}{c\varepsilon}$ holds and it is a saddle point if $a-b\alpha > \frac{eb}{c\varepsilon}$ holds.

For $E_{4}(x_{4}, 0)$,

So, $E_{4}$ is a stable node if $c\varepsilon x_{4} < e$ holds and it is a saddle if $c\varepsilon x_{4} > e$ holds. In particular, $E_{4}$ is a stable node when $a-b\alpha < \frac{eb}{c\varepsilon}$ holds. Notice that $m > a\alpha$, then

Similarly, one can determine that $E_{4}$ is a saddle point when $a-b\alpha > \frac{2eb}{c\varepsilon}$.

In order to discuss the stability of the $E_{3}$, we first need to analyze it in combination with the stability of the $E_{0}$ and $E_{4}$. Assume that $E_{3}$ is a hyperbolic equilibrium point, there are two additional equilibrium points in the system. From Theorem 1, we know that $E_{0}$ is always a stable node when $m > a\alpha$. Therefore, $E_{3}$ is an unstable node if $E_{4}$ is a stable node, and $E_{3}$ is an unstable node or saddle point if $E_{4}$ is a saddle point. In other words, $E_{3}$ is always unstable. □

Remark: In fact, the type of semi-trivial equilibrium point is closely related to the position of the horizontal coordinate $x_{i}$ and the position of the curves $h_{1}(x)$ and $h_{2}(x)$ in the plane Figure 1.

Figure 1.  Four possible positional relationships exist between the curve and the line. The blue curve represents $h_{1}(x) = 0$ ($\lambda_{1} = 0$), while the black line represents $h_{2}(x) = 0$ ($\lambda_{2} = 0$). (a) The semi-trivial equilibrium is either a saddle point if $x_{i}\in (0, \lambda_{2})$ or an unstable node if $x_{i}\in(\lambda_{2}, +\infty)$. (b) The semi-trivial equilibrium can be a stable node if $x_{i}\in (0, \lambda_{1})$, a saddle point if $x_{i}\in (\lambda_{1}, \lambda_{2})$, or an unstable node if $x_{i}\in(\lambda_{2}, +\infty)$. (c) The semi-trivial equilibrium is either a stable node if $x_{i}\in (0, \lambda_{1})$ or an unstable node if $x_{i}\in(\lambda_{1}, +\infty)$. (d) The semi-trivial equilibrium can be a stable node if $x_{i}\in (0, \lambda_{2})$, a saddle point if $x_{i}\in (\lambda_{2}, \lambda_{1})$, or an unstable node if $x_{i}\in(\lambda_{1}, +\infty)$. If the semi-trivial equilibrium point happens to fall at the intersection of the curve and line with the positive half-axis of $x$, then there is a singular trajectory passing through the equilibrium point. It may consist of a stable singular trajectory (left intersection point in (b) and (d)) and an unstable singular trajectory (intersection point in (a) and right intersection point in (b) and (d)), or parallel trajectories with different directions above and below the singular trajectory (intersection point in (c)).

DownLoad: Full-Size Img PowerPoint

Theorem 3. If $a > \frac{be}{c\varepsilon}+\frac{mc\varepsilon}{e+c\alpha\varepsilon}$, then there is a unique positive equilibrium. It is a stable node if $m < b(\frac{e}{c\varepsilon}+\alpha)^{2}$ holds and an unstable node if $m > b(\frac{e}{c\varepsilon}+\alpha)^{2}$ holds.

Proof. From the previous analysis, we know that the system has a unique positive equilibrium $E^{*}(x^{*}, y^{*})$ if $a > \frac{be}{c\varepsilon}+\frac{mc\varepsilon}{e+c\alpha\varepsilon}$ holds. The Jacobian matrix at $E^{*}$ is

where

and

The two eigenvalues of $J(E^{*})$ satisfy

Hence, it can be concluded that the two eigenvalues possess a negative real part when $J_{11}(E^{*}) < 0$, and possess a positive real part when $J_{11}(E^{*}) > 0$. The positive equilibrium can be categorized as a stable node if the condition $m < b(\frac{e}{c\varepsilon}+\alpha)^{2}$ is satisfied. Conversely, when the condition $m > b(\frac{e}{c\varepsilon}+\alpha)^{2}$ is met, the positive equilibrium can be categorized as an unstable node. □

3.   Bifurcation analysis

In this section, we will investigate the occurrence of transcritical bifurcation at $E_{0}$, saddle-node bifurcation at $E_{5}$, Hopf bifurcation at $E^{*}$ and the existence of a limit cycle emerging through Hopf bifurcation.

Theorem 4. If $a\neq b$ holds, then the system (1.1) experiences a transcritical bifurcation at the trivial equilibrium $E_{0}$ as the parameter $m$ passes through the bifurcation value $m = m^{SN} = a\alpha$.

Proof. The Jacobian matrix at $E_{0}$ with $m = m^{SN}$ is

Obviously, matrix $A$ has a simple eigenvalue $\lambda = 0$ with eigenvector $V$, and $A^{T}$ has an eigenvector $W$ corresponding to the eigenvalue $\lambda$. After a simple calculation, we obtain

Therefore,

Hence, according to Sotomayor's theorem in [Perko, 2001], the system (1.1) experiences a transcritical bifurcation at $E_{0}$ as the parameter $m$ passes through the bifurcation value $m^{SN}$. □

Theorem 5. If $a > b\alpha$ and $2be\neq c\varepsilon (a-b\alpha)$ hold, then the system (1.1) experiences a saddle-node bifurcation at the semi-trivial equilibrium $E_{5}$ as the parameter $m$ passes through the bifurcation value $m = m_{SN} = \frac{(a+b\alpha)^{2}}{4b}$.

Proof. The Jacobian matrix at $E_{5}$ with $m = m_{SN}$ is

Obviously, matrix $A$ has a simple eigenvalue $\lambda = 0$ under the given condition. Let $V$ and $W$ be two eigenvectors corresponding to the eigenvalue $\lambda = 0$ for the matrices $A$ and $A^{T}$. After a simple calculation, we obtain

Therefore,

Hence, according to Sotomayor's theorem in [Perko, 2001], system (1.1) experiences a saddle-node bifurcation at $E_{5}$ as the parameter $m$ passes through the bifurcation value $m_{SN}$. □

Theorem 6. If $a > \frac{be}{c\varepsilon}+\frac{mc\varepsilon}{e+c\alpha\varepsilon}$ holds, then the system (1.1) experiences a Hopf bifurcation around $E^{*}$ at $m = m^{*} = b(\frac{e}{c\varepsilon}+\alpha)^{2}$.

Proof. The secular equation of the Jacobian matrix at the equilibrium point $E^{*}$ can be expressed as follows

where

When $m = m^{*}$, we have

The characteristic equation can be reduced to $\lambda^{2}+det(J(E^{*})) = 0.$ The equation has a pair of purely imaginary roots. By a simple computation, we obtain

Hence, from (3.1) and (3.2), it can be concluded that the system (1.1) experiences a Hopf bifurcation around $E^{*}$ at $m = m^{*}$. When the trace of the Jacobian matrix of the system undergoes a sign change from negative to positive, the equilibrium point $E^{*}$ becomes unstable, resulting in the occurrence of Hopf bifurcation. □

Theorem 7. If $a > \frac{be}{c\varepsilon}+\frac{mc\varepsilon}{e+c\alpha\varepsilon}$ and $m^{*} < m < min\{a\alpha, b(x_{1}+\alpha)(x^{*}+\alpha), b(x_{1}+\alpha)^{2}\}$, then the system (1.1) exhibits the presence of a single limit cycle.

Proof. If $a > \frac{be}{c\varepsilon}+\frac{mc\varepsilon}{e+c\alpha\varepsilon}$ and $m^{*} < m < min\{a\alpha, b(x_{1}+\alpha)(x^{*}+\alpha), b(x_{1}+\alpha)^{2}\}$ are satisfied, it can be inferred that $E_{0}$ and $E_{1}$ are saddle points, while $E^{*}$ is an unstable node. Suppose that

where $\bar{x} = \frac{a}{b}$, $\bar{y} > y^{*}$. Obviously, $x^{*} < x_{1} < \bar{x}$. Therefore,

Due to the quadratic coefficient of $g(x)$ being $-\frac{cx^{*}}{2} < 0$, let us now choose a sufficiently large value of $\bar{y}$ (or $\bar{y} > \frac{2}{cx^{*}}(\frac{2(\bar{x}-x^{*})(-e+c\varepsilon \bar{x})}{\beta}+\frac{a^{2}}{4b}-\frac{mx^{*}}{x^{*}+\alpha})$). This ensures that $g(y) < 0$, meaning $\frac{dL_{3}}{dt}\left|\right._{(1.1)} < 0$ for $x^{*} < x < \bar{x}$, $\frac{\bar{y}}{2} < y < \bar{y}$. By the Poincar$\acute{e}$-Bendixson theorem, system (1.1) admits a single limit cycle in the above region. □

4.   The influence of intraspecific cooperation and the Allee effect

Two key questions that need to be answered are how intraspecific cooperation and the Allee effect affect the dynamic behavior of the system. The following discussion provides the answers.

We will first discuss the effect of intraspecific cooperation on the predator. We modeled intraspecific cooperation in the predator population using the term -$\frac{ey}{y+\beta}$. From the previous discussion in Section 2 and Section 3, we know that the parameter $\beta$ will not affect the final states of the local stability of all equilibria, the occurrence of bifurcation at the positive equilibrium, and the existence of a limit cycle around positive equilibria. However, we notice that the predator population net growth rate of the predator population is a monotonically decreasing function with respect to $\beta$,

That is, the growth rate of the predator population decreases gradually with an increase in the parameter $\beta$. Thus, the growth rate of the prey population accelerates while the growth rate of the predator population slows down. This means that the larger the value of $\beta$, the longer it takes for the predator population to reach a stable state. In other words, even when the system is stable, it will take much longer for it to reach a steady state.

We then discuss the effect of the Allee effect on the system. We introduced the Allee effect in the prey population by the term $\frac{mx}{x+\alpha}$. Firstly, due to the existence of the Allee effect in the prey population, the number of equilibria in the system will increase. Secondly, in general, the introduction of the Allee effect in prey population will lead to more complex dynamics of the system. From Theorem 1, we know that $E_{0}$ is a saddle if $\frac{\alpha}{m} > \frac{1}{a}$, and a stable node if $\frac{\alpha}{m} < \frac{1}{a}$. In other words, the local stability of the trivial equilibrium will change from unstable to stable as the Allee effect strength decreases ($E_{0}$ is just a saddle when $m = 0$). The Allee effect can increase the risk of extinction for both predator and prey populations.

From Theorem 3, we know that $E^{*}$ is an unstable node if $\frac{1}{a-bx^{*}}-\frac{x^{*}}{m} < \frac{\alpha}{m} < \frac{\alpha}{b(x^{*}+\alpha)^{2}}$. On the other hand, $E^{*}$ is a stable node if $\frac{\alpha}{m} > max\{\frac{\alpha}{b(x^{*}+\alpha)^{2}}, \frac{1}{a-bx^{*}}-\frac{x^{*}}{m}\}$. The critical value is given by $\frac{\alpha}{m} = \frac{\alpha}{b(x^{*}+\alpha)^{2}}$. The coexistence equilibrium becomes stable as the strength of the Allee effect increases. In other words, the presence of the Allee effect in the prey population can stabilize the previously unstable positive equilibrium. Therefore, the Allee effect plays a crucial role in maintaining the stability of the system.

On the other hand, since the prey population $x^{*} = \frac{e}{c\varepsilon}$ is not affected by the Allee effect strength $\frac{\alpha}{m}$, we will only focus on the impact on the predator population. $y^{*}$ is a continuous function of $\frac{\alpha}{m}$. Rewriting the equation, we have $y^{*} = \frac{1}{c}(a-bx^{*}-\frac{1}{x^{*}/m+\alpha/m})$, then $\frac{dy^{*}}{d(\alpha/m)} = \frac{1}{c}\cdot\frac{1}{(x^{*}/m+\alpha/m)^{2}} > 0$. So, the predator population $y^{*}$ increases with the increase of the Allee effect strength. The predator population becomes greater and greater, but it is still less than $\frac{1}{c}(a-bx^{*})$.

5.   Numerical simulations

Example 1. Let $b = 1$, $c = 2$, $e = 1$, $m = 1$, $\beta = 2$ and $\varepsilon = 0.5$. The unique positive equilibrium point $E^{*}(1, 0.8333)$ of the system is stable when $a = 3$ and $\alpha = 2$ (Figure 2). By performing simple calculations, we can determine that the system has two equilibrium points. The first is a trivial equilibrium point, denoted as $E_{0}(0, 0)$, and the second is a semi-trivial equilibrium point, denoted as $E_{1}(2.7913, 0)$. Both of these equilibrium points are saddle points.

Figure 2.  Coexistence in system (1.1).

DownLoad: Full-Size Img PowerPoint

To understand the influence of the additive Allee effect on the system, let $a = 1.4$ and let $\alpha$ take on various values (Figure 3). Firstly, compared to the absence of the Allee effect, the presence of the Allee effect will certainly cause a decrease in the density of predators, while changes in the density of prey will depend on the strength of the Allee effect. In addition, if the positive equilibrium point of the system in the absence of the Allee effect is stable, then the density of both species is positively correlated with the strength of the Allee effect.In other words, as the strength of the Allee effect increases, the population densities of both predator and prey will increase. Secondly, the Allee effect can act as a destabilizing force. When $\alpha = 2$ and $\alpha = 1.5$, the positive equilibrium point is stable, indicating that the predator and prey can coexist. However, when $\alpha = 1$, the positive equilibrium point loses its stability and the semi-trivial equilibrium point becomes stable, resulting in the extinction of the predator but the survival of the prey. when $\alpha = 0.5$, the semi-trivial equilibrium point loses its stability and the trivial equilibrium point becomes stable, leading to the extinction of both the predator and the prey. Finally, the Allee effect can result in prolonged time for populations to reach their ultimate state of survival, while also accelerating the extinction process for populations that ultimately become extinct.

Figure 3.  Influence of the additive Allee effect.

DownLoad: Full-Size Img PowerPoint

Example 2. The main objective of this example is to analyze the influence of the strength of intraspecific cooperation on the system from the perspective of view of numerical simulations. Let $b = 1$, $c = 2$, $e = 1$, $m = 1$, $\beta = 2$ and $\alpha = 2$. Figure 4 shows that if an equilibrium point (either positive or a boundary equilibrium point) is stable, variations in the strength of intraspecific cooperation do not affect the position and stability of that equilibrium point. However, the presence of varying strength of intraspecific cooperation can result in divergent time for species to achieve their final state. As the strength of intraspecific cooperation increases, the time it takes for predators to reach the final state of survival or extinction becomes longer (Figure 4(a)-4(c)). Similarly, the time it takes for prey to reach the final state of survival also increases (Figure 4(a), 4(b)). On the other hand, the time it takes to reach the final state of extinction becomes shorter (Figure 4(c), 4(d)).

Figure 4.  Influence of intraspecific cooperation.

DownLoad: Full-Size Img PowerPoint

Example 3. Let $a = 3.1$, $b = 1$, $c = 2$, $e = 2$, $m = 3.5$, $\alpha = 1$, $\beta = 5$ and $\varepsilon = 0.5$. In this particular scenario, the system exhibits three equilibrium points. These points include the trivial equilibrium point denoted as $E_{0}(0, 0)$, the semi-trivial equilibrium point $E_{3}(0.2118, 0)$, and the equilibrium point $E_{4}(1.882, 0)$. The simulation results are presented in Figure 5(a). It is evident from the figure that there exist two stable nodes and one unstable saddle point. This finding suggests that, despite the prey population not reaching a density of zero, it is still inadequate to offset the mortality rate of the predator caused by predation. Consequently, the predator is likely to face eventual extinction, despite the energy benefit it gains from consuming prey. Given the inevitable extinction of the predator population, the system described by the system (1.1) will experience degeneration, leading to the emergence of a single population system comprised exclusively of prey:

Figure 5.  The bi-stability behavior of the system.

DownLoad: Full-Size Img PowerPoint

In this case, the system (1.1) is equivalent to the system (5.1). The numerical simulations of system (5.1) are shown in Figure 5(b). In fact, the final state of the prey depends on its initial density. If the initial density falls in the interval $(0, x_{3})$, the prey will eventually become extinct. On the other hand, if the initial density falls in the interval $(x_{3}, +\infty)$, the prey will eventually survive. In summary, the final dynamic behavior of the system (1.1) depends on whether the initial density of the prey exceeds the critical value $x_{3}$. Furthermore, it can be concluded that the presence of the additive Allee effect has a destabilizing effect on the original system. The global asymptotic stability of the positive equilibrium point in the single population Logistic system is contrasted with the local stability of the corresponding positive equilibrium point when the additive Allee effect is present. Furthermore, we can conclude that the additive Allee effect destabilises the original system. This is because the positive equilibrium point of the single population logistic system is globally asymptotically stable, while the corresponding positive equilibrium point with the additive Allee effect is locally stable.

Example 4. Let $a = 3$, $b = 0.9$, $c = 2$, $e = 1$, $m = 2.3$, $\alpha = 1$, $\beta = 2$ and $\varepsilon = 0.9$. In this case, the semi-trivial equilibria $E_{0}(0, 0)$ and $E_{1}(2.6292, 0)$ are saddle points and the coexisting equilibrium $E^{*}(0.5556, 0.5107)$ is an unstable node. The system exhibits one limit cycle. We can clearly observe that the trajectories of an initial value inside and outside the limit cycle approach the limit cycle. The simulation results are shown in Figure 6.

Figure 6.  The limit cycle.

DownLoad: Full-Size Img PowerPoint

Use of AI tools declaration

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

Acknowledgments

This study was funded by the Natural Science Foundation of Fujian Province (No.2021J011155) and the Young and Middle-Aged Scientific research project of Ningde Normal University (No. 2021Q101).

Conflict of interest

The authors declare that they have no conflict of interest.