Stability Analysis and Hopf Bifurcation for ODE System of Predator-Prey Model with Mutual Interference

Khalid Ahmed Abbakar^1,2*, Yafei Yang³, Alhussein Mohamed¹, Songchen Xia¹, Mogahid Mamoon Abkar^2,4, Omer Bushra Elfadil Hassan²

¹College of Mathematics and Statistics, Northwest Normal University, Lanzhou, China.
²Department of Mathematics and Physics, Faculty of Education, University of Gadarif, Gadarif, Sudan.
³School of Mathematics and Statistics, Northeast Normal University, Changchun, China.
⁴Department of Mathematics, School of Mathematics, University of Witwatersrand, Johannesburg, South Africa.

DOI: 10.4236/am.2021.129053                

Abstract

In this paper, the temporal and spatial patterns of a diffusive predator-prey model with mutual interference under homogeneous Neumann boundary conditions were studied. Specifically, first, taking the intrinsic growth rate of the predator as the parameter, we give a computational and theoretical analysis of Hopf bifurcation on the positive equilibrium for the ODE system. As well, we have discussed the conditions for determining the bifurcation direction and the stability of the bifurcating periodic solutions.

Keywords

Predator-Prey Model, Mutual Interference, Hopf Bifurcation, Functional Response

1. Introduction

The interaction with predator-prey populations and their possible outcomes are probably the most studied topics in ecology, because of their existence and universal relevance, will continue to be one of the dominant topics in both ecology and mathematical ecology [1]. In spite of the predator-prey theory has undergone significant developments in the past several years, many long-standing mathematical and ecological problems still deserve researchers’ attention [1] [2] [3] [4]. The functional response of predators to their prey density refers to the change in the density of attached prey per unit time per predator with the change in prey density. It is the average number of prey killed per individual predator per unit of time. For the functional response functions, there are many types. The most commonly used functional response function, which was suggested by Holling (1959) and called Michaelis-Menten type or Holling type II functional re-

sponse [5]. And it is takes the form: $p\left(u\right)=\frac{mu}{a+bu}$. where u represents the density of the prey population, and the positive constants m (units: 1/times) and a (units: 1/prey) describe the effects of capture rate and handling time, respectively, on the feeding rate.

However, Prey-dependent functional responses cannot describe the interaction between predators, and have been challenged by the biology and physiology communities [6] [7]. Some biologists have argued that in many situations, especially when predators have to search for food (and therefore, have to share or compete for food), the functional response in a prey-predator model should be predator-dependent. There is significant evidence that predator dependence in functional response occurs frequently in laboratories and natural systems [8] [9]. Given that many experiments and observations show predators do indeed interfere with each other’s activities to trigger competition effects and that prey changes its behavior under increased predator threat, models with a functional predator-dependent response are reasonable alternatives to models with prey-dependent functional response.

Starting with this argument and the traditional prey-only model, to describe the mutual interference of a predator, Beddington (1975) and DeAngelis (1975) [10] [11] [12] [13] proposed that an individual from a population of over two predators not only allocate time in searching for and processing their prey but also takes time in encountering with other predator. This result in the so-called

Beddington-DeAngelis functional response $p\left(u\mathrm{,}v\right)=\frac{mu}{a+bu+cv}$. The Beddington-DeAngelis is like the well-known Holling type II functional response, but

it has an additional term cv in the denominator modelling mutual interference among predators and has some of the same qualitative features as the ratio-dependent form but avoids some of the singular behaviors of ratio-dependent models at low densities which have been the source of controversy.

The Beddington-DeAngelis functional response and the analogous response with diffusion in a constant environment have received much attention in the literatures [14] [15]. Such as Crowley and Martin [16] assumed that interference among predators occurs regardless of whether a particular predator is looking for prey or handling prey and proposed a functional response, which is called the Crowley-Martin-type functional response, and its takes the form:

$p\left(u,v\right)=\frac{mu}{1+au+bv+abuv}=\frac{mu}{\left(1+au\right)\left(1+bv\right)}$

It is assumed that the predator feeding rate decreases by higher predator density even when prey density is high, and therefore the effects of predator interference in the feeding rate remain important all the time whether an individual predator is handling or searching for a prey at a time. And in [17]. The authors consider a Leslie-Gower predator-prey model with a functional Crowley-Martin response describing predator mutual interference, which its takes the form:

$\{\begin{array}{l}\frac{\text{d}u}{\text{d}t}=ru\left(1-\frac{u}{K}\right)-\frac{muv}{\left(1+au\right)\left(1+bv\right)},\\ \frac{\text{d}v}{\text{d}t}=sv\left(1-\frac{v}{hu}\right),\end{array}$ (1.1)

Initially, their analysis focuses on local asymptotic stability and Hopf bifurcation of the spatially homogeneous model based on ODE system for the reaction-diffusion model.

Based on the above discussion, and therefore. In this paper, we consider the following predator-prey model with mutual interference with Beddington-De Angelis functional response:

$\{\begin{array}{l}\frac{\text{d}u}{\text{d}t}=ru\left(1-\frac{u}{K}\right)-\frac{muv}{a{u}^2+buv+c{v}^2},\\ \frac{\text{d}v}{\text{d}t}=sv\left(1-\frac{v}{hu}\right),\end{array}$ (1.2)

where $u\left(t\right)$ and $v\left(t\right)$ represent population densities of prey and predator at time t, respectively. r, K, s, m, a, b, and h are positive constants, K is the carrying capacity of the prey, and r and s denote to their intrinsic growth rate, respectively,

$\frac{u}{h}$ is a function on the prey population size (h is a measure of the food quality

of the prey for conversion into predator growth depending on the density of the population).

For the model (1.2), we mainly discuss about Hopf bifurcation on the positive equilibrium for the ODE system.

The aim of this article is to show that the diffusive predator-prey model (1.2) shows different spatial, temporal and spatiotemporal patterns across the mechanisms described above.

For simplicity, applying the following scaling:

$\begin{array}{l}\frac{u}{K}\to u,\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}rt\to t,\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}\frac{hKm}{r}\to m,\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}\frac{v}{hK}\to v,\\ h{K}^{2}b\to b,\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}\frac{s}{r}\to s,\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}{h}^2{K}^{2}c\to c,\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}{K}^{2}a\to a\end{array}$

System (1.2) can be simplified by:

$\{\begin{array}{ll}{u}_t-d_1\Delta u=u\left(1-u\right)-\frac{muv}{a{u}^2+buv+c{v}^2},\hfill & x\in \Omega ,\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}t>0,\hfill \\ v_t-d_2\Delta v=sv\left(1-\frac{v}{u}\right),\hfill & x\in \Omega ,\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}t>0,\hfill \\ \frac{\partial u}{\partial \nu }=\frac{\partial v}{\partial \nu }=0,\hfill & x\in \partial \Omega ,\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}t>0,\hfill \\ u\left(x,0\right)=u_0\left(x\right),\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}v\left(x,0\right)=v_0\left(x\right),\hfill & x\in \Omega ,\hfill \end{array}$ (1.3)

where $\Omega \subset {ℝ}^n\left(n\ge 1\right)$ is a bounded smooth domain and $\nu$ is the outward unit normal vector on $\partial \Omega$. The constants $d_1>0$ and $d_2>0$ are the diffusion coefficients of $u\left(x\mathrm{,}t\right)$ and $v\left(x\mathrm{,}t\right)$, which represent the natural dispersive force of movement of the prey and predator density, respectively. The condition of the homogeneous Neumann boundary means that the two species have zero-flux across the boundary $\partial \Omega$. The initial conditions $u_0\left(x\right)>0$, $v_0\left(x\right)>0$ are smooth functions satisfying $u_0\left(x\right)+v_0\left(x\right)>0$.

In section two of this paper, we investigate the asymptotical behaviour of the interior equilibrium and occurrence of Hopf bifurcation of the local ODE system of (1.3).

2. Hopf Bifurcation and Stability Analysis for ODE System

Firstly, we consider the stability of the following ODE system:

$u\left(1-u\right)-\frac{muv}{a{u}^2+buv+c{v}^2}=f\left(u,v\right)$, $sv\left(1-\frac{v}{u}\right)=g\left(u\mathrm{,}v\right)$,

When $f\left(u\mathrm{,}v\right)=g\left(u\mathrm{,}v\right)=0$, we can see that the system has a boundary equilibrium $E_0=\left(\mathrm{1,0}\right)$, two positive interior equilibrium

$\begin{array}{l}{E}_1=\left(u_1,v_1\right)=\left(\frac{\left(a+b+c\right)+\sqrt{\Delta }}{2\left(a+b+c\right)},v_1\right)\\ E_2=\left(u_2,v_2\right)=\left(\frac{\left(a+b+c\right)-\sqrt{\Delta }}{2\left(a+b+c\right)},v_2\right)\end{array}$

if and only if

$\left(H_1\right)a+b+c>4m,$

where $v_1=u_1,v_2=u_2$ and $\Delta ={\left(a+b+c\right)}^2-4m\left(a+b+c\right)$.

By the simple calculations, we note that the boundary equilibrium $E_0=\left(\mathrm{1,0}\right)$ is unstable, we will study the properties of the positive interior equilibrium. Thus, consider the stability of $E_i\left(i=1,2\right)$. Note that the Jacobian of (1.3) at $E_i\left(i=1,2\right)$ is

$J\left(E_i\right)=\left(\begin{array}{cc}1-2u_i+\frac{\left(a-c\right)\left(1-u_i\right)}{a+b+c}& -\frac{\left(a-c\right)\left(1-u_i\right)}{a+b+c}\\ s& -s\end{array}\right).$

For convenience, we denote that $s_0:=1-2u_1+\frac{\left(a-c\right)\left(1-u_1\right)}{a+b+c}$ and $\sigma :=-\frac{\left(a-c\right)\left(1-u_1\right)}{a+b+c}$, so $s_0>0$ if and only if

$\left(H_2\right)\left(2a+b\right)m>{\left(a+b+c\right)}^2{u}_1^2.$

Thus, we can obtain the following theorem.

Theorem 2.1. Suppose that (H₁).

1-1) if $s>s_0$, then $\left(u_1\mathrm{,}{v}_1\right)$ is locally asymptotically stable.

1-2) if (H₂) and $s<s_0$ hold, then $\left(u_1\mathrm{,}{v}_1\right)$ is unstable and there exists at least one periodic orbit for system (1.3).

2) $\left(u_2\mathrm{,}{v}_2\right)$ is unstable.

Proof

1) The characteristic equation corresponding to $\left(u_1\mathrm{,}{v}_1\right)$ is

${\lambda }^2-\text{tr}\left(J\left(E_1\right)\right)\lambda +\det \left(J\left(E_1\right)\right)=\mathrm{0,}$ (2.1)

where $\text{tr}\left(J\left(E_1\right)\right)=s-s_0$ and

$\det \left(J\left(E_1\right)\right)=\frac{\sqrt{{\left(a+b+c\right)}^2-4m\left(a+b+c\right)}s}{a+b+c}>0.$

Therefore, $s>s_0$ yield that $\text{tr}\left(J\left(E_1\right)\right)<0$ and combined with $\det \left(J\left(E_1\right)\right)>0$, we completed the (1-1). Otherwise, if $s<s_0$, then $\left(u_1\mathrm{,}{v}_1\right)$ is unstable.

2) Apply the same argument to the $\left(u_2\mathrm{,}{v}_2\right)$, we have

$\det \left(J\left(E_2\right)\right)=-\frac{\sqrt{{\left(a+b+c\right)}^2-4m\left(a+b+c\right)}s}{a+b+c}<0.$

thus, $\left(u_2\mathrm{,}{v}_2\right)$ is unstable.

In the following, we analyse the existence of Hopf bifurcation at the interior equilibrium $\left(u_1\mathrm{,}{v}_1\right)$, choose s as the bifurcation parameter.

Denote by $\lambda \left(s\right)=\alpha \left(s\right)\pm i\beta \left(s\right)$ the two roots of the characteristic Equation (2.1) with $\alpha \left(s\right)=\frac{1}{2}\left(s_0-s\right)$ and $\beta \left(s\right)=\frac{1}{2}\sqrt{-4s\sigma -{\left(s+s_0\right)}^2}$. Obviously, (2.1) has a pair of imaginary roots $\lambda \left(s\right)=\pm \sqrt{\det \left(J\left(E_1\right)\right)}i$ at $s=s_0$. Moreover, ${\alpha }^{\prime }\left(s_0\right)=-\frac{1}{2}<0$. This shows that the transversality condition holds. Thus undergoes Hopf bifurcation of $\left(u_1\mathrm{,}{v}_1\right)$ as s passes through the $s_0$.

Now, to make a better use of Hopf bifurcation theory, We translate the interior equilibrium $\left(u_1\mathrm{,}{v}_1\right)$ to the origin by the transformation $\tilde{u}=u-u_1$, $\tilde{v}=v-v_1$. For the sake of convenience, we still denote $\tilde{u}$ and $\tilde{v}$ by u and v, respectively. Thus, the system (1.3) is transformed into

$\{\begin{array}{l}\frac{\text{d}u}{\text{d}t}=\left(u+u_1\right)-{\left(u+u_1\right)}^2-\frac{m\left(u+u_1\right)\left(v+v_1\right)}{a{\left(u+u_1\right)}^2+b\left(u+u_1\right)\left(v+v_1\right)+c{\left(v+v_1\right)}^2},\\ \frac{\text{d}v}{\text{d}t}=s\left(v+v_1\right)\left(1-\frac{v+v_1}{u+u_1}\right).\end{array}$ (2.2)

Then, rewrite the system (2.2) as follows:

$\left(\begin{array}{c}\frac{\text{d}u}{\text{d}t}\\ \frac{\text{d}v}{\text{d}t}\end{array}\right)=J\left(E_1\right)\left(\begin{array}{c}u\\ v\end{array}\right)+\left(\begin{array}{c}f\left(u,v,s\right)\\ g\left(u,v,s\right)\end{array}\right),$ (2.3)

where

$\begin{array}{l}f\left(u,v,s\right)=a_1{u}^2+a_{2}uv+a_3{v}^2+a_4{u}^3+a_5{u}^{2}v+a_{6}u{v}^2+a_7{v}^3+\cdots ,\\ g\left(u,v,s\right)=b_1{u}^2+b_{2}uv+b_3{v}^2+b_4{u}^3+b_5{u}^{2}v+b_{6}u{v}^2+\cdots \end{array}$

and

$\begin{array}{l}{a}_1=-1+\frac{m\left(3ac+bc-a^2\right)}{{\left(a+b+c\right)}^3{u}_1^2},\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}{a}_2=\frac{m\left(a^2-ab-6ac-bc+c^2\right)}{{\left(a+b+c\right)}^3{u}_1^2},\\ a_3=\frac{m\left(ab+3ac-c^2\right)}{{\left(a+b+c\right)}^3{u}_1^2},\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}{a}_4=\frac{m\left(a^3-6a^{2}c-4abc+a{c}^2-b^{2}c\right)}{{\left(a+b+c\right)}^4{u}_1^3},\\ a_5=-\frac{m\left(a^3-2a^{2}b-14a^{2}c-4abc+9a{c}^2-b^{2}c+2b{c}^2\right)}{{\left(a+b+c\right)}^4{u}_1^3},\\ a_6=-\frac{m\left(2a^{2}b+9a^{2}c-a{b}^2-4abc-14a{c}^2-2b{c}^2+c^3\right)}{{\left(a+b+c\right)}^4{u}_1^3},\\ a_7=\frac{m\left(a^{2}c-a{b}^2-4abc-6a{c}^2+c^3\right)}{{\left(a+b+c\right)}^4{u}_1^3},\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}{b}_1=-\frac{s}{u_1},\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}{b}_2=\frac{2s}{u_1},\\ b_3=-\frac{s}{u_1},\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}{b}_4=\frac{s}{u_1^2},\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}{b}_5=-\frac{2s}{u_1^2},\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}{b}_6=\frac{s}{u_1^2}.\end{array}$

By setting the matrix

$P:=\left(\begin{array}{cc}N& 1\\ M& 0\end{array}\right),$

where $M=-\frac{s}{\beta }$ and $N=-\frac{s_0+s}{2\beta }$. It is easy to obtain that

$P^{-1}J\left(E_1\right)P=\Phi \left(s\right):=\left(\begin{array}{cc}\alpha \left(s\right)& -\beta \left(s\right)\\ \beta \left(s\right)& \alpha \left(s\right)\end{array}\right).$

When $s=s_0$, we have

$M_0:={M|}_{s=s_0}=-\frac{s_0}{{\beta }_0},\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}{N}_0:={N|}_{s=s_0}=-\frac{s_0}{{\beta }_0},\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}{\beta }_0:=\beta \left(s_0\right)=\sqrt{\det \left(J\left(E_1\right)\right)}.$ (2.4)

By the transformation ${\left(u\mathrm{,}v\right)}^{{\rm T}}=P{\left(x\mathrm{,}y\right)}^{{\rm T}}$, system (2.3) becomes

$\left(\begin{array}{c}\frac{\text{d}x}{\text{d}t}\\ \frac{\text{d}y}{\text{d}t}\end{array}\right)=\Phi \left(s\right)\left(\begin{array}{c}x\\ y\end{array}\right)+\left(\begin{array}{c}{f}^{\prime }\left(x,y,s\right)\\ g^{\prime }\left(x,y,s\right)\end{array}\right),$ (2.5)

where

$\begin{array}{c}{f}^{\prime }\left(x,y,s\right)=\frac{1}{M}g\left(Nx+y,Mx,s\right)\\ =\left(\frac{N^2}{M}{b}_1+N{b}_2+M{b}_3\right)x^2+\left(\frac{2N}{M}{b}_1+b_2\right)xy+\frac{b_1}{M}{y}^2\\ +\left(\frac{N^3}{M}{b}_4+N^2{b}_5+NM{b}_6\right)x^3+\left(\frac{3N^2}{M}{b}_4+2N{b}_5+M{b}_6\right)x^{2}y\\ +\left(\frac{3N}{M}{b}_4+b_5\right)x{y}^2+\frac{b_4}{M}{y}^3+\cdots ,\end{array}$

$\begin{array}{c}{g}^{\prime }\left(x,y,s\right)=f\left(Nx+y,Mx,s\right)-\frac{N}{M}g\left(Nx+y,Mx,s\right)\\ =\left(N^3{a}_4+N^{2}M{a}_5+N{M}^2{a}_6+M^3{a}_7-\frac{N^4}{M}{b}_4-N^3{b}_5-N^{2}M{b}_6\right)x^3\\ +\left(3N^2{a}_4+2NM{a}_5+M^2{a}_6-\frac{3N^3}{M}{b}_4-2N^2{b}_5-NM{b}_6\right)x^{2}y\\ +\left(N^2{a}_1+NM{a}_2+M^2{a}_3-\frac{N^3}{M}{b}_1-N^2{b}_2-NM{b}_3\right)x^2\\ +\left(3N{a}_4+M{a}_5-\frac{3N^2}{M}{b}_4-N{b}_5\right)x{y}^2+\left(a_1-\frac{N}{M}{b}_1\right)y^2\\ +\left(2N{a}_1+M{a}_2-\frac{2N^2}{M}{b}_1-N{b}_2\right)xy+\left(a_4-\frac{N}{M}{b}_4\right)y^3+\cdots \end{array}$

Rewrite the system (2.5) in the following polar coordinate form:

$\begin{array}{l}\stackrel{\dot{}}{\tau }=\alpha \left(s\right)\tau +a\left(s\right){\tau }^3+\cdots ,\\ \stackrel{\dot{}}{\theta }=\beta \left(s\right)\tau +c\left(s\right){\tau }^2+\cdots ,\end{array}$ (2.6)

And then, from the Taylor expansion of (2.6) at $s=s_0$ we have:

$\begin{array}{l}\stackrel{\dot{}}{\tau }={\alpha }^{\prime }\left(s_0\right)\left(s-s_0\right)\tau +a\left(s_0\right){\tau }^3+o\left({\left(s-s_0\right)}^2\tau ,\left(s-s_0\right){\tau }^3,{\tau }^5\right),\\ \stackrel{\dot{}}{\theta }=\beta \left(s_0\right)\tau +{\beta }^{\prime }\left(s_0\right)\left(s-s_0\right)+c\left(s_0\right){\tau }^2+o\left({\left(s-s_0\right)}^2,\left(s-s_0\right){\tau }^2,{\tau }^4\right).\end{array}$ (2.7)

Now, to determine the stability of the Hopf bifurcation periodic solution, So, we need to calculate the sign of the coefficient $a\left(s_0\right)$, which is given by the following formula:

$\begin{array}{l}a\left(s_0\right)=\frac{1}{16}\left(f_{xxx}^1+f_{xyy}^1+g_{xxy}^1+g_{yyy}^1\right)\\ +\frac{1}{16{\beta }_0}\left[f_{xy}^1\left(f_{xx}^1+f_{yy}^1\right)-g_{xy}^1\left(g_{xx}^1+g_{yy}^1\right)-f_{xx}^1{g}_{xx}^1+f_{yy}^1{g}_{yy}^1\right],\end{array}$ (2.8)

where all partial derivatives are evaluated at the bifurcation point $\left(x\mathrm{,}y\mathrm{,}s\right)=\left(\mathrm{0,0,}{s}_0\right)$, and

$\begin{array}{l}{f}_{xxx}^1\left(0,0,s_0\right)=6\left(\frac{N_0^3}{M_0}{b}_4+N_0^2{b}_5+N_0{M}_0{b}_6\right),\\ f_{xyy}^1\left(0,0,s_0\right)=2\left(\frac{3N_0}{M_0}{b}_4+b_5\right),\\ g_{xxy}^1\left(0,0,s_0\right)=2\left(3N^2{a}_4+2N_0{M}_0{a}_5+M_0^2{a}_6-\frac{3N_0^3}{M_0}{b}_4-2N_0^2{b}_5-N_0{M}_0{b}_6\right),\\ g_{yyy}^1\left(0,0,s_0\right)=6\left(a_4-\frac{N_0}{M_0}{b}_4\right),\\ f_{xx}^1\left(0,0,s_0\right)=f_{xy}^1\left(0,0,s_0\right)=0,\end{array}$

$\begin{array}{l}{g}_{xx}^1\left(0,0,s_0\right)\\ =2\left(N_0^2{a}_1+N_0{M}_0{a}_2+M_0^2{a}_3-\frac{N_0^3}{M_0}{b}_1-N_0^2{b}_2-N_0{M}_0{b}_3\right)f_{yy}^1\left(0,0,s_0\right)=\frac{2b_1}{M_0},\\ g_{yy}^1\left(0,0,s_0\right)=2\left(a_1-\frac{N_0}{M_0}{b}_1\right),\\ g_{xy}^1\left(0,0,s_0\right)=\left(2N_0{a}_1+M_0{a}_2-2\frac{N_0^2}{M_0}{b}_1-N_0{b}_2\right).\end{array}$

Noting that $N_0=M_0,b_1=b_3,b_5=-2b_6,b_2=-2b_1$, thus we can calculate that

$\begin{array}{c}a\left(s_0\right)=\frac{\left({\beta }_0^2+s_0^2\right)s_0}{4{\beta }_0^4}{a}_1^2+\frac{\left({\beta }_0^2+3s_0^2\right)s_0}{8{\beta }^4}{a}_1{a}_2+\frac{s_0^3}{4{\beta }_0^4}{a}_1{a}_3-\frac{{\beta }_0^2+s_0^2}{4{\beta }_0^2{s}_0}{a}_1{b}_3\\ +\frac{s_0^3}{8{\beta }_0^4}{a}_2^2+\frac{s_0^3}{8{\beta }_0^4}{a}_2{a}_3-\frac{s_0}{8{\beta }_0^2}{a}_2{b}_3+\frac{3\left({\beta }_0^2+s_0^2\right)}{8{\beta }_0^2}{a}_4-\frac{1}{4}{b}_6\\ +\frac{s_0^2}{4{\beta }_0^2}{a}_5+\frac{s_0^2}{8{\beta }_0^2}{a}_6+\frac{1}{4s_0}{b}_3^2.\end{array}$

Now by calculation, we can determine the value and sign of $a\left(s_0\right)<0$ in (2.8).

Recall that ${\mu }_0=\frac{a\left(s_0\right)}{{\alpha }^{\prime }\left(s_0\right)}$ and ${\alpha }^{\prime }\left(s_0\right)<0$, by Poincaré-Andronov-Hopf Bifurcation Theorem, we have the following result.

Theorem 2.2.

Assume that (H₁) and (H₂) hold. The system (1.3) undergoes a Hopf bifurcation at the coexistence equilibrium $\left(u_1\mathrm{,}{v}_1\right)$ when $s=s_0$. Furthermore:

1) $a\left(s_0\right)$ determines the stability of the bifurcated periodic solutions: if $a\left(s_0\right)<\mathrm{0\; <mtext>\; \hspace{0.17em}\; </mtext>}\left(>0\right)$ then the bifurcating periodic solutions are stable (unstable).

2) ${\mu }_0$ determines the directions of the Hopf bifurcation: if ${\mu }_0>\mathrm{0\; <mtext>\; \hspace{0.17em}\; </mtext>}\left(<0\right)$ then the Hopf bifurcation is supercritical (subcritical).

3. Conclusion

In this article, we showed that the predator-prey model with mutual interference exhibits a rich and interesting dynamic behavior. We first studied the local stability and Hopf bifurcation in the corresponding ODE system. Then we studied the existence and direction of Hopf bifurcation and the stability of the bifurcating periodic solution in the reaction-diffusion system. The system (1.3) undergoes a Hopf bifurcation at the interior equilibrium $\left(u_1\mathrm{,}{v}_1\right)$ when $s=s_0$.

Acknowledgements

The authors are greatly indebted to the anonymous referees for their careful reading and helpful comments, especially the referee who gave constructive suggestions for revision of the manuscript.

Conflicts of Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work re-ported in this paper.

References

[1]	Volkeltigoe, M. (1992) The Orgins and Evolution of Predator-Prey Theory. Ecology, 73, 1530-1535. https://doi.org/10.2307/1940005
[2]	Kuang, Y. (1990) Global Stability of Gause-Type Predator-Prey Systems. Journal of Mathematical Biology, 28, 463-474. https://doi.org/10.1007/BF00178329
[3]	Kuang, Y. (1995) Global Stability in Delay Differential Systems without Dominating Instantaneous Negative Feedbacks. Journal of Differential Equations, 119, 503-532. https://doi.org/10.1006/jdeq.1995.1100
[4]	Yang, K. and Freedman, H.I. (1988) Uniqueness of Limit Cycles in Gause-Type Models of Predator Prey Systems. Mathematical Biosciences, 88, 67-84. https://doi.org/10.1016/0025-5564(88)90049-1
[5]	Skalski, G.T. and Gilliam, J.F. (2001) Functional Responses with Predator Interference: Viable Alternatives to the Holling Type II Model. Ecology, 82, 3083-3092. https://doi.org/10.1890/0012-9658(2001)082[3083:FRWPIV]2.0.CO;2
[6]	Arditi, R. and Ginzburg, L.R. (1989) Coupling in Predator-Prey Dynamics: Ratio-Dependence. Journal of Theoretical Biology, 139, 311-326. https://doi.org/10.1016/S0022-5193(89)80211-5
[7]	Arditi, R. and Saiah, H. (1993) Empirical Evidence of the Role of Heterogeneity in Ratio-Dependent Consumption. Ecology, 74, 2020. https://doi.org/10.2307/1940847
[8]	Dolman, P.M. (1995) The Intensity of Interference Varies with Resource Density: Evidence from a Field Study with Snow Buntings, Plectrophenax nivalis. Oecologia, 102, 511-514. https://doi.org/10.1007/BF00341364
[9]	Jost, C. and Ellner, S.P. (2000) Testing for Predator Dependence in Predator-Prey Dynamics: A Non-Parametric Approach. Proceedings of the Royal Society of London. Series B: Biological Sciences, 267, 1611-1620. https://doi.org/10.1098/rspb.2000.1186
[10]	Beddington, J.R. (1975) Mutual Interference between Parasites or Predators and Its Effect on Searching Efficiency. The Journal of Animal Ecology, 44, 331-340. https://doi.org/10.2307/3866
[11]	Antwi-Fordjour, K., Parshad, R.D. and Beauregard, M.A. (2020) Dynamics of a Predator-Prey Model with Generalized Holling Type Functional Response and Mutual Interference. Mathematical Biosciences, 326, Article ID: 108407. https://doi.org/10.1016/j.mbs.2020.108407
[12]	DeAngelis, D.L., Goldstein, R.A. and O’Neill, R.V. (1975) A Model for Tropic Interaction. Ecology, 56, 881-892. https://doi.org/10.2307/1936298
[13]	Hwang, T.-W. (2004) Uniqueness of Limit Cycles of the Predator-Prey System with Beddington-DeAngelis Functional Response. Journal of Mathematical Analysis and Applications, 290, 113-122. https://doi.org/10.1016/j.jmaa.2003.09.073
[14]	Cantrell, R.S. and Cosner, C. (2001) Effects of Domain on the Persistence of Populations in a Diffusive Food-Chain Model with Beddington-Deangelis Functional Response. Natural Resource Modeling, 14, 335-367. https://doi.org/10.1111/j.1939-7445.2001.tb00062.x
[15]	Bai, D.Y., Yu, J.S., Fan, M. and Kang, Y. (2020) Dynamics for a Non-Autonomous Predator-Prey System with Generalist Predator. Journal of Mathematical Analysis and Applications, 485, Article ID: 123820. https://doi.org/10.1016/j.jmaa.2019.123820
[16]	Crowley, P.H. and Martin, E.K. (1989) Functional Responses and Interference within and between Year Classes of a Dragonfly Population. Journal of the North American Benthological Society, 8, 211-221. https://doi.org/10.2307/1467324
[17]	Rahmi, E., Darti, I. and Suryanto, A. (2021) A Modified Leslie-Gower Model Incorporating Beddington-DeAngelis Functional Response, Double Allee Effect and Memory Effect. Fractal and Fractional, 5, 84. https://doi.org/10.3390/fractalfract5030084