2 Stability analysis of system (1.1)

In this section, the future behaviour of prey and predator populations under the Allee effect will be investigated. For that, we first determine the equilibrium points of system (1.1). It is clear that $(\bar {N},\bar {P})=(0,0)$ is an equilibrium point of this system, but there is nothing biologically interesting about that equilibrium point. The model has been designed with the assumption that the prey population is the only source of food for the predator population, and the survival of the predator is dependent on the number of prey in the environment. Then, system (1.1) has an equilibrium point of the form

(2.1) $$ \begin{align} (\bar{N},\bar{P})=\bigg(\frac{\theta r_1}{\epsilon r_2}-\beta,\frac{\theta r_1}{\epsilon r_2}-\beta\bigg). \end{align} $$

Notice that this equilibrium point is positive so that the Allee coefficient stated in (2.1) can be taken as $ \beta =k{\theta r_1}/{\epsilon r_2} $ for $k \in \mathbb {R} $ and $k\in (0,1) $ . As a result,

$$ \begin{align*} (\bar{N},\bar{P})=\bigg((1-k)\frac{\theta r_1}{\epsilon r_2},(1-k)\frac{r_1}{\epsilon}\bigg) \end{align*} $$

is the unique positive equilibrium point of system (1.1) for a constant value of k. If we look closely, as the value of k increases, the Allee constant $\beta $ increases, so that the equilibrium point diminishes due to the Allee effect. This is in line with the biological interpretation of the decrease in the prey population that cannot find a mate owing to the Allee effect and, accordingly, the decrease in the predator population. The Jacobean matrix of system (1.1) at its steady state is

$$ \begin{align*} J(\bar{N},\bar{P})= \left(\begin{array}{@{}cc@{}} r_1 k (1-k) & -(1-k)\displaystyle\frac{\theta r_1}{r_2} \\[1mm] \displaystyle\frac{r^2_2}{\theta} & -r_2 \end{array}\right), \end{align*} $$

where $ tr(J)=r_1 k (1-k)-r_2 $ and $ \det (J)=r_1r_2(1-k)^2 $ . Let us define $ a_1:= -tr (J) $ and $ a_2:= \det (J) $ . So, we can write the characteristic polynomial associated with the Jacobean matrix as follows:

(2.2) $$ \begin{align} F(\lambda)=\lambda^2-tr (J) \lambda+\det (J) = \lambda^2+a_1\lambda+a_2. \end{align} $$

Then, the equilibrium point $ (\bar {N},\bar {P}) $ is locally asymptotically stable when $a_1>0$ (that is, $tr({J})<0$ ) and $a_2>0$ (that is, $\det ({J})>0$ ). Obviously, $a_2 =\det (J) =r_1r_2(1-k)^2$ is always positive. Therefore, let us determine the conditions required for the coefficient $a_1$ to be positive.

The equation $a_1=r_2-r_1k(1-k)=r_1k^2-r_1k+r_2$ is a second-order polynomial according to the parameter k. Let us express this equation by $ a_1(k) $ . When $r^2_1-4r_1r_2\geq 0$ , the equation $a_1(k)$ has two real roots of the form:

(2.3) $$ \begin{align} k_{1,2}=\frac{r_1 \mp \sqrt{r^2_1-4r_1r_2}}{2r_1}. \end{align} $$

We are dealing with two situations here. If $r_1>4r_2$ , then $r^2_1-4r_1r_2> 0$ and then $k_{1,2}>0$ is satisfied since $ \sqrt {r^2_1-4r_1r_2} < \lvert r_1 \rvert $ . Moreover, since

(2.4) $$ \begin{align} k_1=\frac{r_1 - \sqrt{r^2_1-4r_1r_2}}{2r_1} \quad \text{and} \quad k_2=\frac{r_1 + \sqrt{r^2_1-4r_1r_2}}{2r_1} \end{align} $$

are written as above, $k_1<k_2 $ is true and $a_1(k)>0 $ under the condition that $k \in (0,k_1) \cup (k_2, 1) $ . However, if $r_1=4r_2 $ , then $k_{1,2}={1}/{2}$ and $a_1(k)=0$ . Much as the characteristic polynomial $ F(\lambda ) $ has a pair of purely imaginary roots when $r_1=4r_2$ , that is, $a_1(k)=0 $ , $a_2>0 $ and $\lambda _{1,2}=\mp i \sqrt {a_2}$ , the stability condition of the equilibrium point does not change as k changes. Consequently, with the next theorem, we can express the required conditions for the stable point of system (1.1) to be a sink.

3 Bifurcation analysis

We have shown in Section 2 that the stability structure of the equilibrium point $ (\bar {N},\bar {P}) $ of system (1.1) changes as parameter k and, accordingly, the Allee constant, $\beta $ , changes. It may therefore be said that there are qualitative changes, namely bifurcations, in the dynamics of system (1.1) under the influence of Allee. Here, a change in the stability condition of the equilibrium point means that the negative real part of the eigenvalues is zero at the same value, which is called the bifurcation parameter, and the real part of the same eigenvalues is positive immediately following the bifurcation parameter. This occurs if the system has zero or a pair of purely imaginary eigenvalues, and Hopf bifurcation occurs in the system when there is a pair of purely imaginary eigenvalues. So, let us investigate under which conditions system (1.1) has Hopf bifurcation under the Allee effect.

Proof. Let us consider that $ \lambda = i w $ is a root of $ F(\lambda ) $ for $ w \in \mathbb {R} $ and $ w>0 $ . Since $ \lambda = i w $ satisfies the equation of $F(\lambda )=0$ , we obtain the following:

$$ \begin{align*} \begin{aligned} \lambda^2+a_1\lambda+a_2=0 &\Rightarrow (iw)^2+a_1(iw)+a_2=0 \Rightarrow -w^2+a_1iw+a_2=0 \\ &\Rightarrow (a_2-w^2)+i(a_1w)=0 \Leftrightarrow a_1=0 \hspace{0.2cm} ve \hspace{0.2cm} a_2=w^2>0. \end{aligned} \end{align*} $$

Now, suppose that $a_1=0 $ and $a_2>0 $ . Then, from $ F(\lambda ) = \lambda ^2+a_2 $ , we get the roots as $ \lambda _{1,2}= \mp i \sqrt {a_2} $ . Thus, we have concluded the proof.

Proof. First, if $k=k_1$ or $k=k_2$ , then $a_1=0$ . However, from Lemma 3.1, polynomial $F(\lambda )$ has a pair of purely imaginary eigenvalues since $a_2>0$ .

Now, suppose that the characteristic polynomial $ F(\lambda ) $ has a pair of complex conjugate eigenvalues in the form of $\lambda _{1,2}(k)=m(k)\mp in(k)$ , where $m(k), n(k) \in \mathbb {R}$ and $n(k) \neq 0$ . Here, parameter k is taken as a bifurcation parameter and to prove the existence of Hopf bifurcation arising from the stable point $(\bar {N},\bar {P})$ of system (1.1), we need to show that $m(k_c)=0$ and $n(k_c)=n_0>0$ are satisfied by choosing $k_c $ as a critical bifurcation parameter.

We next assume that $ \lambda _1= m(k) + in(k) $ , and take $ m(k)=m $ and $ n(k)=n $ to make representation easier. Since $ \lambda _1 $ is a root of $ F(\lambda ) $ ,

$$ \begin{align*} \begin{aligned} (m+in)^2+a_1(m+in)+a_2=0 &\Rightarrow m^2+2mni-n^2+a_1m+ia_1n+a_2=0 \\ &\Rightarrow (m^2-n^2+a_1m+a_2)+i(2mn+a_1n)=0. \end{aligned} \end{align*} $$

Hence, the equalities $ m^2-n^2+a_1m+a_2=0 $ and $ 2mn+a_1n =0 $ must be satisfied simultaneously. In a similar way, for $ \lambda _2= m- in $ ,

$$ \begin{align*} \begin{aligned} (m-in)^2+a_1(m-in)+a_2=0 &\Rightarrow m^2-2mni-n^2+a_1m-ia_1n+a_2=0 \\ &\Rightarrow (m^2-n^2+a_1m+a_2)+i(-2mn-a_1n)=0. \end{aligned} \end{align*} $$

Again, we have the following equations to provide the last equality:

(3.1) $$ \begin{align} n(2m+a_1)=0 & \Rightarrow a_1=-2m , \end{align} $$

(3.2) $$ \begin{align} m^2-n^2+a_1m+a_2=0 & \Rightarrow m^2+a_1m+a_2-n^2=0. \end{align} $$

If we substitute (3.1) into (3.2), we obtain the above equation,

$$ \begin{align*} m^2+m(-2m)+a_2-n^2=0 \Rightarrow m^2=a_2-n^2. \end{align*} $$

Obviously, we are dealing with three situations here.

• • If $a_2<n^2 $ , then $ m^2<0 $ and there is a contradiction to the fact that $ m \in \mathbb {R} $ .

• • If $a_2=n^2 $ , then $ m=0 $ and then $ \lambda _{1,2}=\mp i \sqrt {a_2} $ .

• • If $a_2>n^2 $ , then $ m=\mp \sqrt {a_2-n^2} \neq 0$ . However, in this situation, one cannot find any critical value $k_c $ that solves the equation $ m(k_c)=0 $ .

Let us now bear in mind (3.2). This equation is a second-order polynomial with real coefficients depending on the parameter $ m $ . Moreover, when $a_1^2-4(a_2-n^2)\geq 0 $ , (3.2) has two real roots as stated below:

$$ \begin{align*} \begin{aligned} m_1=\frac{-a1+\sqrt{a_1^2-4(a_2-n^2)}}{2} \quad \text{and}\quad m_2=\frac{-a1-\sqrt{a_1^2-4(a_2-n^2)}}{2}. \end{aligned} \end{align*} $$

If we pay attention, $a_1,a_2 $ and $ n $ are functions of parameter k, and $a_2 $ is always positive. So, in the next cases, we will look at the scenarios where the coefficient $a_1 $ is zero, negative or positive.

Case 1: If $a_1=0 $ , then from (3.1), we have $ m=0 $ and putting this value into (3.2), $ n^2=a_2 \Rightarrow n= \mp \sqrt {a_2}$ is found. Thus, the polynomial $ F(\lambda ) $ has a pair of purely imaginary eigenvalues in the form of $ \lambda _{1,2}=\mp i \sqrt {a_2} $ .

Case 2: Again consider (3.1); if $a_1>0 $ , then $ m=- {a_1}/{2} <0 $ must be satisfied. Moreover, let us consider the roots $ m_1 $ and $ m_2 $ separately for (3.2), and to do this, first assume that $a_1^2-4(a_2-n^2)>0 $ holds.

• • If $a_2=n^2 $ , then $ m_1=0 $ and this is inconsistent with $ m<0 $ .

• • If $a_2<n^2 $ , then $ m_1>0 $ and again it goes counter to $ m<0 $ .

• • If $a_2>n^2 $ , then $ m_1<0 $ is found.

As a result, if $ m$ is taken as $ m=m_1 $ when $a_2>n^2 $ under the condition that $a_1>0 $ and $a_1^2-4(a_2-n^2)>0 $ , then the polynomial $ F(\lambda ) $ has a complex conjugate root whose real part is negative. Therefore, the equilibrium point is locally asymptotically stable. Similarly:

• • $ m_2=-a_1<0 $ holds when $a_2=n^2 $ ;

• • $ m_2<0 $ whether $a_2<n^2 $ or $a_2>n^2 $ .

Then, for the values $a_1>0 $ and $a_1^2-4(a_2-n^2)>0 $ , if $ m=m_2 $ is chosen, the polynomial $ F(\lambda ) $ has a complex conjugate root with negative real part and this means that the equilibrium point of system (1.1) is a sink.

Case 3: If $a_1<0 $ , then $ m=-{a_1}/{2}>0 $ from (3.1). Also, $a_1^2-4(a_2-n^2)>0 $ must be held for the real value of $ m $ . Under these assumptions:

• • if $a_2=n^2 $ , then $ m_1=-a_1>0 $ ;

• • if $a_2<n^2 $ or $a_2>n^2 $ , then $ m_1>0 $ .

Therefore, if $ m=m_1 $ is taken when $a_1<0 $ or $a_1^2-4(a_2-n^2)>0 $ , the equilibrium point becomes unstable since the real part of the polynomial $ F(\lambda ) $ has a positive complex conjugate root. Now, let us investigate the sign of $ m_2 $ under the same conditions:

• • $ m_2=0 $ if $a_2=n^2 $ holds, but this contradicts with $ m>0 $ ;

• • $ m_2<0 $ is found for $a_2<n^2 $ and again there is a contradiction to the hypothesis $ m>0 $ ;

• • $ m_2>0 $ if $a_2>n^2 $ is taken.

As a result of the above findings, $ F(\lambda ) $ has a complex conjugate root with the positive real part if we choose $ m=m_2 $ when $a_2>n^2 $ under the assumption that $a_1<0 $ and $a_1^2-4(a_2-n^2)>0 $ . This indicates that the equilibrium point is unstable.

Notice that $a_1 $ , that is $a_1(k) $ , determines the critical bifurcation value. Therefore, the solution of $a_1(k)=0$ will give the critical bifurcation value. Let us denote this critical bifurcation value as $k_c $ . When $r_1>4r_2 $ , if $k_c=k_1 $ or $k_c=k_2 $ are taken (here, $k_1$ and $k_2 $ are the roots of the equation $a_1(k)=0$ are given as (2.3)), then $ m(k_c)=0 $ , $ n(k_c)=\sqrt {a_2} $ , and consequently $ \lambda _{1,2}= \mp i \sqrt {a_2} $ is found. This completes the proof.

The characteristic polynomial $ F(\lambda ) $ has a pair of purely imaginary eigenvalues under the condition $ (C2) $ given in Theorem 3.2. Also, when $ 0<k<k_1 $ or $k_2<k<1 $ (that is, the condition $ (C1) $ of Theorem 2.1 holds), the equilibrium point of system (1.1) is locally asymptotically stable (see Case 2). However, if $k_1<k<k_2 $ , then the equilibrium point of system (1.1) is unstable (see Case 3).

Now, let us give the following results for the complex eigenvalues we achieved under the above assumptions.

Proof. The proof is completed by taking the derivative of both sides of the equation $a_1 = -2m$ . Indeed, for

$$ \begin{align*} \frac{d(a_1)}{dk}=(2r_1k-r_1), \quad \frac{d(a_1)}{dk}\bigg|_{k_c} \neq 0, \end{align*} $$

since $k_c\neq {1}/{2}.$ If we take a close look, $r_1=4r_2 $ must be held to satisfy the equation

$$ \begin{align*} k_c=k_{1,2}=\frac{r_1 \mp \sqrt{r^2_1-4r_1r_2}}{2r_1}= \frac{1}{2}.\end{align*} $$

However, $r_1>4r_2 $ from the definition of $k_c $ . So, $k_c $ cannot be equal to $ 1/2 $ .

3.1 Hopf bifurcation

In the former section, we have concluded that the positive equilibrium point $ (\bar {N},\bar {P})=((1-k){\theta r_1}/{\epsilon r_2},(1-k){r_1}/{\epsilon }) $ of system (1.1) is locally asymptotically stable by Theorem 2.1. We also know from Theorem 3.2 that the characteristic polynomial $ F(\lambda ) $ has a pair of purely imaginary eigenvalues if $k=k_1 $ or $k=k_2 $ . Hence, Hopf bifurcation occurs when $r_1>4r_2 $ in this case.

We now show the existence of the Hopf bifurcation by using the centre manifold theorem and bifurcation theory given in [Reference Kuznetsov23] (see the Appendix). First, let us shift the equilibrium point $ (\bar {N},\bar {P}) $ of system (1.1) to the origin for simplicity and, for this purpose, denote the right-hand sides of the equations in system (1.1) as follows:

(3.3) $$ \begin{align} \begin{aligned} f(N(t),P(t)):&= r_1N(t)\frac{N(t)}{N(t)+\beta}-\epsilon P(t)N(t);\\ g(N(t),P(t)):&=r_2P(t)-\theta\frac{P^2(t)}{N(t)}. \end{aligned} \end{align} $$

After that, using Taylor expansion of f and g defined by (3.3) at $(\bar {N},\bar {P})$ and following transformations

$$ \begin{align*} \begin{aligned} U(t)&=N(t)-\bar{N} ,\\[1mm] V(t)&=P(t)-\bar{P} ,\\[1mm] \tilde{k}&=k-k_c, \end{aligned} \end{align*} $$

we obtain the following system with equilibrium point $(0, 0)$ :

(3.4) $$ \begin{align} \left(\begin{array}{@{}c@{}} \dot{U}(t) \\[2mm] \dot{V}(t) \end{array}\right)= \left(\begin{array}{@{}cc@{}} r_1 k_c (1-k_c) & -(1-k_c)\displaystyle\frac{\theta r_1}{r_2} \\[1mm] \displaystyle\frac{r^2_2}{\theta} & -r_2 \end{array}\right)\left(\begin{array}{@{}c@{}} U(t) \\[1mm] V(t) \end{array}\right)+\left(\begin{array}{c} f_1(U(t),V(t),\tilde{k}, k_c)\\[1mm]f_2(U(t),V(t),\tilde{k}, k_c) \end{array}\right), \end{align} $$

where

$$ \begin{align*} f_1(U(t),V(t),\tilde{k}, k_c)&=r_1(\tilde{k}-\tilde{k}^2-2\tilde{k}k_c) U(t) +\tilde{k} \frac{\theta r_1}{r_2} V(t) \\ &\quad +\frac{1}{2}\bigg[2(\tilde{k}+2\tilde{k}k_c+k^2_c)\frac{\epsilon r_2}{\theta}U^2(t)-2\epsilon U(t)V(t)\bigg]\\ &\quad +\frac{1}{6}\bigg[-6(\tilde{k}+2\tilde{k}k_c+k^2_c)\frac{\epsilon^2 r^2_2}{\theta^2r_1}U^3(t)\bigg]+ \text{H.O.T.},\\[1mm] f_2(U(t),V(t),\tilde{k}, k_c)&=\frac{1}{2} \bigg[-\frac{2\epsilon r^3_2}{(1-(\tilde{k}+k_c))\theta^2 r_1}U^2(t)+ \frac{4 \epsilon r_2^2}{(1-(\tilde{k}+k_c))\theta r_1} U(t)V(t) \\ &\quad-\frac{ 2\epsilon r_2}{(1-(\tilde{k}+k_c))r_1}V^2(t)\bigg] +\frac{1}{6}\bigg[\frac{6\epsilon^2 r^4_2}{(1-(\tilde{k}+k_c))^2 \theta^3 r^2_1}U^3(t) \\ &\quad -\frac{12\epsilon^2 r^3_2}{(1-(\tilde{k}+k_c))^2 \theta^2 r^2_1}U^2(t)V(t) +\frac{ 2\epsilon^2 r_2^2}{(1-(\tilde{k}+k_c))^2 \theta r_1^2}U(t)V^2(t)\bigg] \\ &\quad +\text{H.O.T.}, \end{align*} $$

in which H.O.T. represents the higher order terms in $f_1$ and $f_2$ . Notice that not only the equilibrium point but also the critical bifurcation value $ \tilde {k} $ (when $k=k_c $ ) has shifted to the origin in system (3.4). Therefore, we can rewrite the functions $f_1$ and $f_2$ as follows:

$$ \begin{align*} \begin{aligned} f_1(U(t),V(t),k_c)&=s_1 k^2_c U^2(t)-s_2 U(t)V(t) -s_3k^2_cU^3(t); \\[1mm] f_2(U(t),V(t),k_c)&=-\frac{s_4}{(1-k_c)}U^2(t) +\frac{s_5}{(1-k_c)} U(t)V(t) - \frac{s_6}{(1-k_c)} V^2(t)\\[1mm] &\quad+\frac{s_7}{(1-k_c)^2}U^3(t) -\frac{s_8}{(1-k_c)^2} U^2(t)V(t)+\frac{s_9}{(1-k_c)^2}U(t)V^2(t). \end{aligned} \end{align*} $$

Here, the coefficients $ s_i $ for $ i=1, \ldots , 9 $ are stated:

$$ \begin{align*} \begin{aligned} s_1&=\frac{\epsilon r_2}{\theta} ,\quad s_2=\epsilon ,\qquad s_3=\frac{\epsilon^2 r^2_2}{\theta^2r_1}, \qquad s_4=\frac{\epsilon r^3_2}{\theta^2 r_1},\quad s_5=\frac{2 \epsilon r_2^2}{\theta r_1},\\[1mm] s_6&=\frac{\epsilon r_2}{r_1}, \quad s_7=\frac{\epsilon^2 r^4_2}{ \theta^3 r^2_1},\quad s_8= \frac{2\epsilon^2 r^3_2}{ \theta^2 r^2_1},\quad s_9=\frac{\epsilon^2 r_2^2}{\theta r_1^2}. \end{aligned} \end{align*} $$

Now, assume that $k_c=k_1 $ , then $ \lambda _{1,2}(k_c)=m(k_c)\mp in(k_c), $ where $ m(k_c)=0 $ and $ n(k_c)=\sqrt {a_2} $ based on Theorem 3.2. In addition, Lemma 3.4 is also satisfied and the conditions of (H.1) and (H.2) of the Hopf bifurcation theorem are held. So, let us reduce system (3.4) to normal form so that we can apply the Hopf bifurcation theorem and calculate the Lyapunov coefficient $ c_1(k_c) $ , which gives the directional analysis of this bifurcation.

The Jacobian matrix of system (3.4) at the critical bifurcation value is

$$ \begin{align*} J(k_c)=\left(\begin{array}{@{}cc@{}} r_1 k_c (1-k_c) & -(1-k_c)\displaystyle\frac{\theta r_1}{r_2} \\[2mm] \displaystyle\frac{r^2_2}{\theta} & -r_2 \end{array}\right):=\left(\begin{array}{@{}cc@{}} j_{11} & j_{12} \\[2mm] j_{21} & j_{22} \end{array}\right). \end{align*} $$

We now construct a matrix T whose columns consist of real and imaginary parts of an eigenvector corresponding to $ \lambda _1=-in(k_c) $ or $ \lambda _1=-in $ in short as follows:

$$ \begin{align*} T=\left( \begin{array}{@{}cc@{}} j_{12} & 0 \\[1mm] -j_{11} & -n \end{array}\right). \end{align*} $$

Using the following transformation:

$$ \begin{align*} \left( \begin{array}{@{}c@{}} U(t) \\[1mm] V(t) \end{array} \right)=T \left( \begin{array}{@{}c@{}} X(t) \\[1mm] Y(t) \end{array} \right), \end{align*} $$

system (3.4) is reduced to the following system of equations:

(3.5) $$ \begin{align} \left( \begin{array}{@{}c@{}} \dot{X}(t) \\[1mm] \dot{Y}(t) \end{array} \right)= \left( \begin{array}{@{}cc@{}} 0 & -n \\[1mm] n & 0 \end{array}\right) \left( \begin{array}{@{}c@{}} X (t)\\[1mm] Y(t) \end{array} \right) + \left( \begin{array}{@{}c@{}} \tilde{f_1}(X(t),Y(t),k_c) \\[2mm] \tilde{f_2}(X(t),Y(t),k_c) \end{array} \right), \end{align} $$

where

$$ \begin{align*} \begin{aligned} \tilde{f_1}(X(t),Y(t),k_c)&= t_1X^2(t)+t_2X(t)Y(t)+t_3Y^2(t) +\Theta(\lvert X(t)\rvert+\lvert Y(t) \rvert+ \lvert k_c\rvert)^3;\\[1mm] \tilde{f_2}(X(t),Y(t),k_c)&=l_1X^2(t)+l_2X(t)Y(t)+l_3Y^2(t) +l_4X^3(t)+l_5X^2(t)Y(t) \\[1mm] &\quad+l_6X(t)Y^2(t)+\Theta(\lvert X(t)\rvert+\lvert Y(t) \rvert+ \lvert k_c\rvert)^4 , \end{aligned} \end{align*} $$

in which the coefficients are given as follows:

$$ \begin{align*} \begin{aligned} t_1&=s_1j_{12}k^2_c +s_2j_{11};\\[1mm] t_2&=s_2 n;\\[1mm] t_3&=-s_3j^2_{12}k^2_c;\\[1mm] l_1&=\frac{s_4j^2_{12}}{(1-k_c)n}-\frac{s_1j_{11}j_{12}k^2_c, }{n}-\frac{s_2j^2_{11}}{n}+\frac{s_5j_{11}j_{12} }{(1-k_c)n}+\frac{s_6j^2_{11} }{(1-k_c)n};\\[1mm] l_2&=\frac{s_5j_{12}}{(1-k_c)}-s_2j_{11}+\frac{2s_6j_{11}}{(1-k_c)};\\[1mm] l_3&=\frac{s_6}{(1-k_c)n}; \\[1mm] l_4&=\frac{s_3j_{11}j^2_{12}k^2_c }{n}-\frac{s_7j^2_{12} }{(1-k_c)n}-\frac{s_8j_{11}j^2_{12} }{(1-k_c)^2n}-\frac{s_9j^2_{11}j_{12} }{(1-k_c)^2n}; \\[1mm] l_5&=-\frac{s_8j^2_{12} }{(1-k_c)^2}-\frac{2s_9j_{11}j_{12} }{(1-k_c)^2}; \\[1mm] l_6&=-\frac{s_9j_{12}n }{(1-k_c)^2}. \end{aligned} \end{align*} $$

Again, considering system (3.5) and applying the complex transformation $ Z(t)=X(t)+iY(t) $ , we obtain

(3.6) $$ \begin{align} \begin{aligned} \dot{Z}(t)&=\dot{X}(t)+i\dot{Y}(t)\\[1mm] &=-nY(t)+t_1X^2(t)+t_2X(t)Y(t)+t_3X^3(t)\\[1mm] &\quad +i(nX(t)+l_1X^2(t)+l_2X(t)Y(t)+l_3Y^2(t)+l_4X^3(t)+l_5X^2(t)Y(t)+l_6X(t)Y^2(t))\\[1mm] &=in(X(t)+iY(t)) +\tilde{f_1}(X(t),Y(t),k_c)+i\tilde{f_2}(X(t),Y(t),k_c)\\[1mm] &=\lambda (k_c)Z(t) + g(Z(t),\bar{Z}(t),k_c), \end{aligned}\nonumber\\[-16pt] \end{align} $$

where

$$ \begin{align*} \begin{aligned} g(Z(t),\bar{Z}(t),k_c)&=\frac{g_{20}}{2} Z^2(t)+g_{11}Z(t)\bar{Z}(t)+\frac{g_{02}}{2} \bar{Z}^2(t)+\frac{g_{12}}{2} Z(t)\bar{Z}^2(t)\\[1mm] &\quad +\frac{g_{21}}{2} Z^2(t)\bar{Z}(t) +\frac{g_{30}}{6} Z^3(t) +\frac{g_{03}}{6} \bar{Z}^3(t). \end{aligned} \end{align*} $$

Here, the coefficients are calculated as

$$ \begin{align*} \begin{aligned} \frac{g_{20}}{2}&=\frac{t_1}{4}-i\frac{t_2}{4}+i\frac{l_1}{4}+\frac{l_2}{4}-i\frac{l_3}{4}, \\[1mm] g_{11}&=\frac{t_1}{2}+i\frac{l_1}{2}+i\frac{l_3}{2},\\[1mm] \frac{g_{02}}{2}&=\frac{t_1}{4}+i\frac{t_2}{4}+i\frac{l_1}{4}-\frac{l_2}{4}-i\frac{l_3}{4}, \\[1mm] \frac{g_{12}}{2}&=\frac{3t_3}{8}+i\frac{3l_4}{8}-\frac{l_5}{8}+i\frac{l_6}{8}, \\[1mm] \frac{g_{21}}{2}&=\frac{3t_3}{8}+i\frac{3l_4}{8}+\frac{l_5}{8}+i\frac{l_6}{8}, \\[1mm] \frac{g_{30}}{6}&=\frac{t_3}{8}+i\frac{l_4}{8}+\frac{l_5}{8}-i\frac{l_6}{8}, \\[1mm] \frac{g_{03}}{6}&=\frac{t_3}{8}+i\frac{l_4}{8}-\frac{l_5}{8}-i\frac{l_6}{8}. \end{aligned} \end{align*} $$

Now, we reduce system (3.6) to its normal form in the following two steps.

Step 1. Since $ \lambda (k_c)=in(k_c) $ in (3.6), both $ \lambda =in $ and $ \bar {\lambda }=-in $ are not equal to zero. Hence, for a sufficiently small value of k, using an invertible parameter-dependent complex coordinate transformation

$$ \begin{align*} Z(t)=Q(t)+\frac{h_{20}}{2}Q^2(t)+h_{11}Q(t)\bar{Q}(t)+\frac{h_{02}}{2}\bar{Q}^2(t), \end{align*} $$

and taking

(3.7) $$ \begin{align} h_{20}=\frac{g_{20}}{\lambda},\quad h_{11}=\frac{g_{11}}{\bar{\lambda}} \quad \text{and} \quad h_{02}=\frac{g_{02}}{2\bar{\lambda}-\lambda}, \end{align} $$

(3.6) is obtained as the following equation without quadratic terms:

(3.8) $$ \begin{align} \dot{Q}(t)= \lambda Q(t) +\frac{\tau_{30}}{6} Q^3(t)+ \frac{\tau_{12}}{2} Q(t)\bar{Q}^2(t)+ \frac{\tau_{21}}{2} Q^2(t)\bar{Q}(t)+\frac{\tau_{03}}{6} \bar{Q}^3(t) , \end{align} $$

where the coefficients are given as

$$ \begin{align*} \frac{\tau_{30}}{6}&=\frac{g_{30}}{6}-\lambda h_{11}\frac{h_{02}}{2}-h_{20}\frac{g_{20}}{2}-h_{11}\frac{g_{02}}{2},\nonumber\\ \frac{\tau_{12}}{2}&=2(\bar{\lambda}-\lambda)h_{02}h_{11}-\lambda h_{11}\bigg(\frac{h_{20}}{2}+h_{11} \bigg)+\frac{g_{12}}{2}-h_{20}\frac{g_{02}}{2}\nonumber\\ &\quad -h_{11}\bigg(g_{11}+\frac{g_{20}}{2}\bigg)-h_{02}g_{11},\\ \frac{\tau_{21}}{2}&=(\bar{\lambda}-\lambda)h^2_{02}-\lambda h_{11}\bigg(h_{11}+\frac{h_{20}}{2}\bigg)+\frac{g_{21}}{2}-h_{20}g_{11}\nonumber\\ &\quad -h_{11}\bigg(g_{11}+\frac{g_{20}}{2}\bigg)-h_{02}\frac{g_{02}}{2},\nonumber \\ \frac{\tau_{03}}{6}&=\frac{g_{03}}{6}-\lambda h_{11}\frac{h_{02}}{2}+(\bar{\lambda}-\lambda)h^2_{02}-h_{02}\frac{g_{02}}{2}-h_{11}\frac{g_{02}}{2}.\nonumber \end{align*} $$

Step 2. Again, for sufficiently small values of k, using an invertible parameter- dependent change of complex coordinate transformation as follows:

$$ \begin{align*} Q(t)=W(t)+\frac{h_{30}}{6}W^3(t)+\frac{h_{12}}{2}W(t)\bar{W}^2(t)+ \frac{h_{21}}{2}W^2(t)\bar{W}(t)+\frac{h_{03}}{6}\bar{W}^3(t)+ \text{H.O.T.}, \end{align*} $$

and choosing

$$ \begin{align*} h_{30}=\frac{\tau_{30}}{2\lambda},\quad h_{12}=\frac{\tau_{12}}{2\bar{\lambda}} \quad \text{ and}\quad h_{03}=\frac{\tau_{03}}{3\bar{\lambda}-\lambda}, \end{align*} $$

cubic terms disappear except $ W^2(t)\bar {W}(t) $ and (3.8) can transform into the following equation:

(3.9) $$ \begin{align} \dot{W}(t)=\lambda W(t) +\tfrac{1}{2}(\tau_{21}-(\lambda+\bar{\lambda}) h_{21})W^2(t)\bar{W}(t) + \text{H.O.T.}. \end{align} $$

Notice that if $ h_{21}={\tau _{21}}/({\lambda +\bar {\lambda }}) $ is chosen, then the leading coefficient of the term $ W^2(t)\bar {W}(t) $ can be zero. However, this option cannot be made because $ \lambda (k_c)=\sqrt {a_2} $ and $\bar {\lambda }(k_c)=-\sqrt {a_2} $ would be $ \lambda +\bar {\lambda }=0 $ . Thus, since $ \lambda +\bar {\lambda }=0$ in (3.9), the coefficient of the term $ W^2(t)\bar {W}(t) $ is found as ${\tau _{21}}/{2}$ and taking this coefficient as $ c_1 $ , one obtains the next equation:

$$ \begin{align*} \dot{W}(t)=\lambda W(t) +c_1W^2(t)\bar{W}(t) + \text{H.O.T.}, \end{align*} $$

where the coefficient $ c_1=c_1(k_c) $ has the following form using (3.7) and (3.8) with $ \lambda (k_c)=in(k_c)=in $ and $ \bar {\lambda }(k_c)=-in(k_c)=-in $ :

$$ \begin{align*} c_1(k_c)&=\frac{g_{20}g_{11}(2\lambda+\bar{\lambda})}{2 \lvert \lambda \rvert^2}+\frac{\lvert g_{11} \rvert^2}{\lambda}+\frac{\lvert g_{20} \rvert^2}{2(2\bar{\lambda}-\lambda)}+\frac{g_{21}}{2}\notag \\[1mm] &=\frac{i}{2n}g_{20}g_{11}-\frac{i}{n}\lvert g_{11} \rvert^2-\frac{i}{6n}\lvert g_{20} \rvert^2+\frac{g_{21}}{2}. \end{align*} $$

As it is known, the mate-finding Allee effect is particularly important for species that are sensitive to low population densities. Heretofore, the stability structure of system (1.1) was examined depending on the k parameter and, consequently, the Allee constant $\beta $ value, while the other parameters remained constant. Does the Allee effect, that is, the parameter $\beta $ , play a more effective role in the dynamics of system (1.1) than the other parameters? To what extent is system (1.1) sensitive to the parameters? In the former section, we attempted to answer these questions.

4 Sensitivity analysis

The parameters of the equations that comprise the model determine how species interact with one another and vary over time in the prey–predator system. So, we concentrate on a sensitivity analysis to find out how resilient the model is to parameter values that are associated with the $k_1 $ and $k_2 $ , at which the dynamics of the system changes. Additionally, to do this, we calculate the sensitivity indices according to the sensitivity approach given in [Reference Martcheva27] and derive an analytical expression for the sensitivity index (elasticity) of $k_1 $ and $k_2 $ with respect to each parameter according to the explicit formula of these values as follows:

$$ \begin{align*} \begin{aligned} \bullet \hspace{0.1cm} E^{r_1}_{k_1}&=\frac{\partial k_1}{\partial r_1} \times \frac{r_1}{k_1}=-\frac{2r_1r_2}{r_1\sqrt{r^2_1-4r_1r_2}-(r^2_1-4r_1r_2)};\\ \bullet \hspace{0.1cm} E^{r_2}_{k_1}&=\frac{\partial k_1}{\partial r_2} \times \frac{r_2}{k_1}=\frac{2r_1r_2}{r_1\sqrt{r^2_1-4r_1r_2}-(r^2_1-4r_1r_2)};\\ \bullet \hspace{0.1cm} E^{r_1}_{k_2}&=\frac{\partial k_2}{\partial r_1} \times \frac{r_1}{k_2}=\frac{2r_1r_2}{r_1\sqrt{r^2_1-4r_1r_2}+(r^2_1-4r_1r_2)};\\ \bullet \hspace{0.1cm} E^{r_2}_{k_2}&=\frac{\partial k_2}{\partial r_2} \times \frac{r_2}{k_2}=-\frac{2r_1r_2}{r_1\sqrt{r^2_1-4r_1r_2}+(r^2_1-4r_1r_2)}. \end{aligned} \end{align*} $$

Since $k_1,k_2 \in (0,1) $ and $r_1,r_2>0 $ , we have the following corollary.

According to the above result, since $E^{r_1}_{k_1} < 0$ , the increase in $r_1$ , which expresses the growth rate of the prey population, will reduce $k_1$ , that is, the Allee effect. It is true that a higher $r_1$ would result in a higher prey population birth rate and, consequently, a higher rate of mate finding, which would mitigate the Allee effect. However, since $E^{r_2}_{k_1}> 0$ , an increase in $r_2$ , which represents the growth rate of the predator population, will result in a decrease in the number of prey and a lower chance of finding a mate. For other sensitivity formulae, similar inferences can be made.

5 Numerical examples and simulations

In this section, the theoretical results found in Sections 2, 3 and 4 will be supported by numerical examples, and the impact of the Allee effect on the dynamics of the given system will be demonstrated. To perform numerical calculations, we used the MATLAB ODE package.

To begin with, we consider the next example to exhibit the local stability of system (1.1). We use the parameter values inspired by [Reference Baydemir, Merdan, Karaoglu and Sucu7], which are $r_1=0.5$ , $r_2=0.08$ , $\epsilon =0.03$ and $\theta =0.05$ .

Notice that the equilibrium point of system (1.1) with constant $r_1,r_2,\epsilon $ and $ \theta $ parameter values is

$$ \begin{align*} (\bar{N},\bar{P})=\bigg((1-k)\frac{\theta r_1}{\epsilon r_2},(1-k)\frac{r_1}{\epsilon}\bigg), \end{align*} $$

and it varies depending on the k parameter. In addition, as the value of k increases (correspondingly, as the Allee constant increases), the equilibrium point decreases. In Table 1, the Allee coefficient $\beta $ and equilibrium point $ (\bar {N},\bar {P}) $ of system (1.1) are given, corresponding to some k values.

Table 1 The Allee constant and equilibrium point of system (1.1) according to the chosen k values.

Since $r_1>4r_2 $ , $k_1=0.2 $ and $k_2=0.8 $ , if we take $k=0.02 $ , then condition $ (C1) $ holds in Theorem 2.1. Therefore, based on the fixed parameter values above, the unique equilibrium point $ (\bar {N}, \bar {P})=(10.21,16.33) $ is locally asymptotically stable. We have chosen the initial conditions as $N(0)=10$ and $P(0)=15$ for the numerical simulations of system (1.1). Figure 1(a) shows that if k is taken as less than $k_c =k_1$ , then the equilibrium point is locally asymptotically stable. As seen in Figure 1(b), the equilibrium point loses its stability and periodic solutions occur when k is equal to the critical bifurcation value $k_1 $ . However, Figure 1(c) illustrates that the equilibrium point is unstable if we select $k=0.3 $ as larger than $k_1 $ .

Figure 1 Density-time graphs of prey (left), density-time graphs of predator (middle) and phase portrait of system (1.1) (right) for the different values of k: (a) $k=0.02 $ ; (b) $k=0.2 $ and (c) $k=0.3 $ . The initial conditions are $ N(0)=10 $ and $ P(0)=15 $ for each simulation.

In Table 1, one can easily see that the equilibrium point is calculated as $ (\bar {N},\bar {P})=(8.3333, 13.3333) $ and the Allee constant is calculated as $ \beta =2.0833 $ for the crucial bifurcation value of $k_1=0.2 $ . At this critical bifurcation value, system (1.1) has a pair of purely imaginary roots as $ \lambda _{1,2}=m(k_1)\mp i n(k_1)=\mp 0.16 i $ . Figure 2(a) and (b) depict the decline in the equilibrium number of prey and predators as the Allee constant rises, and “H” is the point that designates the Allee value at which Hopf bifurcation appears. In addition to these, in accordance with the formulae in the Hopf bifurcation theorem expressed in the Appendix, the coefficients giving the stability, direction and period of the Hopf bifurcation are calculated as follows:

$$ \begin{align*} &\star\quad c_1(k_1)=-0.0011- 0.0003i;\\ &\star \quad Re(c_1(k_1))=-0.0011<0;\\ &\star \quad m' (k_1)=0.15\neq 0;\\ &\star \quad \alpha(k_1)=-\frac{Re(c_1(k_1))}{m' (k_1)}=0.0075>0;\\ &\star \quad n' (k_1)=-0.2;\\ &\star \quad \tilde{T} = -\frac{Im (c_1(k_1))+\alpha(k_1) n'(k_1)}{n(k_1)}=0.0114>0. \end{align*} $$

Here, $ m'(k_1)>0$ and $\alpha (k_1)>0 $ . So, the direction of the Hopf bifurcation is supercritical, and bifurcating periodic solutions are stable since $ Re(c_1(k_1))<0 $ . Moreover, the period of the limit cycles increases since $ \tilde {T}>0 $ . The aforementioned Hopf bifurcation diagram is shown in Figure 2(c). This diagram was created by using the MATCONT package, which is a graphical MATLAB software, for which further details can be found in [Reference Dhooge, Govaerts and Kuznetsov15].

Figure 2 Panels (a) and (b) represent prey and predator equilibria as the Allee coefficient $\beta $ increases, and point H presents the Hopf bifurcation point. (c) Limit cycles that undergo Hopf bifurcation in $(\beta , N, P) $ space.

Now, let us examine Figure 3, which includes the solution graphs of system (1.1) and the phase portrait of this system for the values of k that are greater than the first critical bifurcation value $k_1 $ at the same parameter values and initial conditions. Figure 3 presents that the Allee effect changes the stability of the equilibrium point according to each specific k value from unstable to stable. One can also see from Figure 3(a) and (b) that equilibrium point $ (\bar {N},\bar {P})=(3.1250, 5.0000) $ is unstable for $k=0.7<k_2=0.8 $ , and periodic solutions arise in a neighbourhood of the equilibrium point $ (\bar {N},\bar {P})=(2.0833, 3.3333) $ when k passes through another critical bifurcation value $k_c=k_2=0.8 $ .

Figure 3 The trajectory of prey and predator density of system (1.1) versus time and the phase portrait of prey density versus predator density, respectively, with the initial conditions $ N(0)=10 $ and $ P(0)=15 $ when (a) $k=0.7 $ , (b) $k=0.8 $ and (c) $k=0.9 $ .

Finally, the equilibrium point $ (\bar {N},\bar {P})=(1.0417, 1.6667) $ is stable for the larger values of the bifurcation values and this case is simulated for $k=0.9>k_2 $ in Figure 3(c). Therefore, Hopf bifurcation occurs at the critical bifurcation value of $k_2 $ . Indeed, for $k=k_2=0.8 $ , the eigenvalues of system (1.1) have the form of $ \lambda _{1,2}=m(k_2)\mp i n(k_2)=\mp 0.04 i $ .

Again, similar calculations are done to determine the stability, direction and period of the Hopf bifurcation for the value of $k_2 $ :

$$ \begin{align*} &\star \quad c_1(k_2)=0.0139 - 1.4367i,\\ &\star \quad Re(c_1(k_2))=0.0139>0,\\ &\star \quad m' (k_2)=-0.15\neq 0,\\ &\star \quad \alpha(k_2)=-\frac{Re(c_1(k_2))}{m' (k_2)}=0.0925>0,\\ &\star \quad n' (k_2)=-0.2,\\ &\star \quad \tilde{T} = -\frac{Im (c_1(k_2))+\alpha(k_2) n'(k_2)}{n(k_2)}=-35.4544<0. \end{align*} $$

The fourth case in the Hopf bifurcation theorem exists because $ m'(k_2)<0 $ and $\alpha (k_2)>0$ . Since $Re(c_1(k_2))>0$ , unstable periodic solutions are seen. Also, $ \tilde {T} -35.4544<0 $ , and as can be seen from Figure 3, periodic solutions are unstable (that is, subcritical) with a decreasing period.

Note that the Allee effect has a significant impact on the dynamical behaviour of system (1.1). Since critical values $k_1$ and $k_2$ (accordingly Allee constant) are functions of the parameters $r_1$ and $r_2$ , it is expected that densities of populations are sensitive to the changing of these values. We provide numerical evidence for this notion in Figure 4, which also serves as a numerical example to bolster Corollary 4.1. It should be noted that in the provided sensitivity analysis example, the other parameter values are kept constant while the impact of a single parameter value change on the system solution is examined.

Figure 4 (a,c) Trajectories of prey and predator densities versus time when $r_1=0.5 $ (square), $r_1=0.55 $ (diamond) and $r_1=0.65 $ (star), while $r_2, \epsilon , \theta $ are constant values in the case of $k=k_1 $ . (b,d) Trajectories of prey and predator densities versus time where $r_2=0.08$ , $0.09$ , $0.12$ are marked with square, diamond and star, respectively, while other parameters are fixed in the same case. The initial conditions are $N(0)=10$ and $P(0)=15$ for these simulations.

The question is, which parameter, when all other parameter values fluctuate over time, has the greatest effect on the populations of prey and predators? The FAST (Fourier amplitude sensitivity test) approach of the MATLAB SAFE (Sensitivity analysis for everyone) program [Reference Pianosi, Sarrazin and Wagener35] was used to provide an answer to this query. Here, the sensitivity index of the model parameters was calculated based on the ranges given in Tables 2 and 3 for the variation of parameter values, and the results are displayed in Figures 5 and 6. As can be seen from Figure 5, when the parameter values change so that they remain in the stability region of system (1.1), the most effective parameter for the system is the $r_1$ parameter, which shows the growth rate of the prey. Changes in the $r_1$ and $r_2$ parameters also impact the $\beta $ parameter since $\beta =k{\theta r_1}/{\epsilon r_2}$ . Consequently, it seems inevitable that the dynamics of system (1.1) would be more susceptible to variations in these three parameters.

Table 2 Upper and lower limits for parameter values in system (1.1).

Table 3 Upper and lower limits for parameter values in system (1.1) where prey population is low.

Figure 5 (a) Time-dependent sensitivity index of each parameter of system (1.1) for the parameter ranges in Table 2. (b) Sensitivity index bar graph.

Figure 6 (a) Time-dependent sensitivity index of each parameter of system (1.1) for the parameter ranges in Table 3. (b) Sensitivity index bar graph.

Another issue that occurs in real life is the potential for the prey population to not increase enough once the system is stable. Table 3 shows example values where the prey growth rate $r_1$ is reduced and its range is restricted.

For these values, the sensitivity indices are shown in Figure 6. Here, the Allee effect $\beta $ has a greater sensitivity index. From a biological perspective, the Allee effect drives a species towards extinction if the prey population grows slowly or does not have a sufficient growth rate, which results in difficulties for the population in finding a mate and reproducing. In this case, the dynamics of the system is more sensitive to changes in the Allee parameter compared with other parameters. However, we would anticipate a decrease in the Allee effect if the prey population is sufficiently large, the reproductive rate is high and finding a partner is not difficult. This is in line with Figure 5.