Bifurcation and Turing Pattern Formation in a Diffusion Modified Leslie-Gower Predator-Prey Model with Crowley-Martin Functional Response

Dong Wang, Yani Ma

College of Mathematics and Statistics, Northwest Normal University, Lanzhou, China.

DOI: 10.4236/jamp.2024.126133                

Abstract

In this paper, we study a modified Leslie-Gower predator-prey model with Smith growth subject to homogeneous Neumann boundary condition, in which the functional response is the Crowley-Martin functional response term. Firstly, for ODE model, the local stability of equilibrium point is given. And by using bifurcation theory and selecting suitable bifurcation parameters, we find many kinds of bifurcation phenomena, including Transcritical bifurcation and Hopf bifurcation. For the reaction-diffusion model, we find that Turing instability occurs. Besides, it is proved that Hopf bifurcation exists in the model. Finally, numerical simulations are presented to verify and illustrate the theoretical results.

Keywords

Modified Leslie-Gower Model, Crowley-Martin Function Response, Hopf Bifurcation, Transcritical Bifurcation, Turing Instability

1. Introduction

The interaction between predators and prey is an important process in maintaining ecological balance in ecosystems. To better understand the dynamic behavior of this interaction, researchers have proposed many mathematical models to describe the evolution of predator-prey systems. In recent years, the evolution of the population system has an important reference significance for the protection of the ecosystem and has obtained [1]-[17] many valuable results. The Leslie-Gower model is a classic predator-prey model that can describe how the populations of predators and prey change over time [3] [4]. However, this model assumes that the mortality rate of the prey population is independent of density, which may not be consistent with reality in some cases. To address this issue, researchers have introduced diffusion terms to modify the Leslie-Gower model, allowing it to better describe the death process of the prey population [18]-[23]. The Crowley-Martin function is a commonly used form to describe the response strength of predators to changes in prey density. This function has been widely applied in predator-prey models.

The growth function of the traditional logistic model is $h\left(u\right)=r\left(1-\frac{u}{k}\right)$ , where r represents the internal growth rate of the prey without the action of the predator. It is noted that the logistic equation $\frac{\text{d}u}{\text{d}t}=ru\left(1-\frac{u}{k}\right)$ and the average growth rate $\frac{u^{\prime }\left(t\right)}{u}$ is a linear function of the prey. However, under the action of

environmental poisons, this assumption is unrealistic for food limited populations. In the early 1960s, Smith [7] conducted experiments on large dapsophila populations in the laboratory, and the results showed that the hypothesis of linear growth rate was inconsistent with large dapsophila populations, and demonstrated that population growth needed food to sustain, but only needed food to maintain when the population reached saturation. Therefore, Smith modified the Logistic model and formed the “limited food” model, sometimes called the “limited resources” model [8] [9] [10] [11] [24]. The model is given as follows

$\frac{\text{d}u}{\text{d}t}=ru\frac{k-u}{k+au}\mathrm{,}$

where, a represents the resource limitation parameter of the population and $\frac{r}{a}$

represents the mass replacement rate of the population at k. This is a further development of the Logistic model. Adding the thinking of the above questions, Yue et al. [8] discussed the properties of the model solution by studying the diffusion Holling-Tanner predator-prey model with Smith growth, and obtained the existence of non-constant of the steady-state system. Jiang et al. [11] mainly studied the diffusion delay model with Smith growth and group behavior, analyzed the existence and stability of Hopf bifurcation in this model, and verified the theoretical results by numerical simulation. Han et al. [24] discussed the dynamic behavior of a space and time discrete predator-prey system with Smith growth function. By the stability analysis, the parametric conditions are gained to ensure the stability of the homogeneous steady state of the system. Xiaozhou Feng et al. [25] proposed a modified Leslie-Gower predator-prey ODE model with Smith growth rate and B-D functional response term, and the systematic study was conducted on the dynamic behavior of the model. Combining the advantages of the above models, the following models are proposed by us:

$\{\begin{array}{l}\frac{\text{d}u}{\text{d}t}\mathrm{<mo>=</mo>}{r}_{1}u\frac{1-u/k}{1+au/k}-\frac{Muv}{\left(1+Bu\right)\left(1+Cv\right)}\mathrm{,}\\ \frac{\text{d}v}{\text{d}t}\mathrm{<mo>=</mo>}{r}_{2}v\left(1-\frac{Dv}{u+E}\right)\mathrm{,}\end{array}$ (1)

where u and v are prey and predator densities and r₁ and r₂ denote their intrinsic growth rates, respectively; a represents the resource limitation parameter of the population, M is the rate at which the predator consumes the prey, B and C are normal numbers, D and E represent the conversion rate of prey into predators biomass and the carrying capacity of the population respectively, $r_1\mathrm{,<mi>\; </mi>}{r}_2\mathrm{,<mi>\; </mi>}a\mathrm{,<mi>\; </mi>}M\mathrm{,<mi>\; </mi>}B\mathrm{,<mi>\; </mi>}C\mathrm{,<mi>\; </mi>}D\mathrm{,<mi>\; </mi>}E$ are positive constants, and a is a non-negative constant, if $a=0$ , we retain a classical Holling type II Lotka-Volterra model with logistic growth of the prey.

For simplicity, we introduce the dimensionless variables as in [2],

$u\mapsto ku\mathrm{,}v\mapsto v\mathrm{,}t\mapsto \frac{t}{r_1}\mathrm{,}$

the dimensionless parameters

$b=BK,c=C,s=\frac{r_2}{r_1},d=\frac{D}{k},e=\frac{E}{k}$

system (1) can be simplified as follows:

$\{\begin{array}{l}\frac{\text{d}u}{\text{d}t}=\frac{u\left(1-u\right)}{1+au}-\frac{muv}{\left(1+bu\right)\left(1+cv\right)},\\ \frac{\text{d}v}{\text{d}t}=sv\left(1-\frac{dv}{u+e}\right),\\ u\left(0\right)=u_0>0,v\left(0\right)=v_0>0.\end{array}$ (2)

In real world, the spatial distributions of the predator and prey are inhomogeneous within a fixed bounded domain, and each species has a nature tendency to diffuse to areas of smaller population concentration [26]. Hence, we should use reaction-diffusion equations to describe spatial dispersal of each species, so the system (2) is further modified into

$\{\begin{array}{l}\frac{\partial u}{\partial t}=d_1\Delta u+\frac{u\left(1-u\right)}{1+au}-\frac{muv}{\left(1+bu\right)\left(1+cv\right)},\hfill \\ \frac{\partial v}{\partial t}=d_2\Delta v+sv\left(1-\frac{dv}{u+e}\right),\hfill \\ {\partial }_{\nu }u={\partial }_{\nu }v=0,\hfill \\ u\left(0\right)=u_0>0,v\left(0\right)=v_0>0.\hfill \end{array}$ (3)

where $\nu$ is the outer unit normal vectors of the boundary. The homogeneous Neumann boundary condition shows that the predator-prey system is self-contained and the population flow on the boundary is zero. The positive constants d₁ and d₂ are diffusion coefficients, and the initial values are nonnegative continuous functions.

The rest of this paper is organized as follows. In Section 2, we analyze the local asymptotic stability, and discuss various common bifurcation problems. In Section 3, we firstly consider the Turing instability of the coexistence equilibrium for reaction-diffusion system (3). Then we discuss the global asymptotic stability of the coexistence equilibrium for reaction-diffusion system (3).

2. Stability and Bifurcations of the ODE Model

2.1. Equilibria and Their Stability

Let

$\{\begin{array}{l}f\left(u,v\right)=\frac{u\left(1-u\right)}{1+au}-\frac{muv}{\left(1+bu\right)\left(1+cv\right)}=0,\\ g\left(u,v\right)=sv\left(1-\frac{dv}{u+e}\right)=0.\end{array}$ (4)

All the equilibrium solutions of the system are as follows

(1) $E_0=\left(\mathrm{0,0}\right)$ , which indicates both prey and predator are extinct.

(2) $E_1=\left(\mathrm{1,0}\right)$ , which indicates the predator is extinct.

(3) $E_2=\left(\mathrm{0,}\frac{e}{d}\right)$ , which indicates the prey is extinct.

(4) $E^{\mathrm{<mo>*</mo>}}=\left(u^{\mathrm{<mo>*</mo>}}\mathrm{,}{v}^{\mathrm{<mo>*</mo>}}\right)$ , which indicates the coexistence of prey and predator populations.

Where $\left(u^{\mathrm{<mo>*</mo>}}\mathrm{,}{v}^{\mathrm{<mo>*</mo>}}\right)$ is positive equilibrium point of (2), then $v_{\mathrm{<mo>*</mo>}}=\frac{1}{d}\left(e+u^{\mathrm{<mo>*</mo>}}\right)$ , $u^{\mathrm{<mo>*</mo>}}$ here is satisfied the following cubic equations

$f\left(u\right):={\rho }_3{u}^3+{\rho }_2{u}^2+{\rho }_{1}u+{\rho }_0=0,$ (5)

where

$\begin{array}{l}{\rho }_3=bc>0,\\ {\rho }_2=ma+c+bd+bce-bc,\\ {\rho }_1=m+mae+d+ce-bce-c-bd,\\ {\rho }_0=me-d-ce.\end{array}$

Similar to the references [27] [28] we obtain:

$\begin{array}{l}{\Delta }^{*}={\left(\frac{p}{3}\right)}^3+{\left(\frac{q}{2}\right)}^2,p=\frac{3{\rho }_3{\rho }_1-{\rho }_2^2}{3{\rho }_3^2},\\ q=\frac{27{\rho }_3^2{\rho }_0-9{\rho }_1{\rho }_2{\rho }_3+2{\rho }_2^3}{27{\rho }^3},\end{array}$

Case 1: When ${\Delta }^{*}>0$ , system (2) has a positive equilibrium $E^{\mathrm{<mo>*</mo>}}=\left(u^{\mathrm{<mo>*</mo>}}\mathrm{,}{v}^{\mathrm{<mo>*</mo>}}\right)$ .

Case 2: When ${\Delta }^{\mathrm{<mo>*</mo>}}=0$ ,

(i) $p=0$ ($\Rightarrow$ $q=0$ and $bc>ma+c+bd+bce$ ), system (2) has a triple real root:

$\overline{E}=\left(\overline{u}\mathrm{,}\overline{v}\right)=\left(\frac{bc-ma-c-bd-bce}{3bc}\mathrm{,}\frac{bc+2bce-ma-c-bd}{3bcd}\right)\mathrm{.}$

(ii) $p<0$ , system (2) has a non-degenerate single real root $\overline{E_1}$ and a degenerate double real root $\overline{E_2}$ .

Case 3: When ${\Delta }^{*}<0$ , system (2) has three unequal positive real roots $E_1^{\mathrm{<mo>*</mo>}}$ , $E_2^{\mathrm{<mo>*</mo>}}$ , $E_3^{\mathrm{<mo>*</mo>}}$ .

Next, the stability of the system is studied on this equilibrium solution.

By calculating the eigenvalues of Jacobin matrix J on (2) as follows

$J=\left(\begin{array}{ll}\frac{1-2u-a{u}^2}{{\left(1+au\right)}^2}-\frac{mv}{{\left(1+bu\right)}^2\left(1+cv\right)}\hfill & \mathrm{<mi>\; </mi>}-\frac{mu}{\left(1+bu\right){\left(1+cv\right)}^2}\hfill \\ \mathrm{<mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi>}\frac{ds{v}^2}{{\left(e+u\right)}^2}\hfill & \mathrm{<mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi>}s-\frac{2sdv}{u+e}\hfill \end{array}\right)\mathrm{,}$

We can establish the following conclusions.

(1) The eigenvalues of Jacobin matrix J at $E_0=\left(\mathrm{0,0}\right)$ are 1 and $-\delta$ , so equilibrium solution $E_0=\left(\mathrm{0,0}\right)$ is a saddle point, and it is unstable.

(2) The eigenvalues of Jacobin matrix J at $E_1=\left(\mathrm{1,0}\right)$ points the eigenmultinomial matrix is

$J\left(E_1\right)=\left(\begin{array}{ll}-\frac{1}{1+a}\hfill & \mathrm{<mi>\; </mi>0}\hfill \\ \mathrm{<mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi>0}\hfill & \mathrm{<mi>\; </mi>}s\hfill \end{array}\right)\mathrm{.}$

The characteristic polynomial corresponding to this matrix is

${\lambda }^2+\left(\frac{1}{1+a}-s\right)\lambda -\frac{s}{1+a}=0.$

Since $-\frac{s}{1+a}<0$ , so that $E_1$ is a saddle point, and the equilibrium point $E_1=\left(\mathrm{1,0}\right)$ is unstable.

(3) The eigenvalues of Jacobin matrix at point $E_2=\left(\mathrm{0,}\frac{e}{d}\right)$ points the eigenmultinomial matrix is

$J\left(E_2\right)=\left(\begin{array}{ll}1-\frac{e}{d+ce}\hfill & \mathrm{<mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi>0}\hfill \\ \mathrm{<mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi>}\frac{s}{d}\hfill & \mathrm{<mi>\; </mi>}-s\hfill \end{array}\right)\mathrm{.}$

The characteristic polynomial corresponding to this matrix is

${\lambda }^2+\left(s+\frac{e}{d+ce}-1\right)\lambda +\frac{se}{d+ce}-s=0.$

If $\frac{se}{d+ce}-s>0$ , there $e>d+ce$ , $s+\frac{e}{d+ce}-1>0$ , therefore, the equilibrium point $E_2=\left(\mathrm{0,}\frac{e}{d}\right)$ is asymptotically stable.

If $\frac{se}{d+ce}-s<0$ , there $e<d+ce$ , it follows that $E_2=\left(\mathrm{0,}\frac{e}{d}\right)$ is a saddle point, so the equilibrium point $E_2=\left(\mathrm{0,}\frac{e}{d}\right)$ is unstable.

(4) We mainly discuss the stability of the positive equilibrium point $\left(u^{\mathrm{<mo>*</mo>}}\mathrm{,}{v}^{\mathrm{<mo>*</mo>}}\right)$ . By calculation and bifurcation theory, we establish the following theorem.

In the following, we analyze the stability of the positive equilibria of model (2). Let $E^{\mathrm{<mo>*</mo>}}=\left(u^{\mathrm{<mo>*</mo>}}\mathrm{,}{v}^{\mathrm{<mo>*</mo>}}\right)$ be a positive equilibrium of model (2). Then the Jacobin matrix of model (2) at $E^{\mathrm{<mo>*</mo>}}$ is

$J\left(E^{\mathrm{<mo>*</mo>}}\right)=\left(\begin{array}{ll}{s}_0\hfill & \mathrm{<mi>\; </mi><mi>\; </mi>}\sigma \hfill \\ \frac{s}{d}\hfill & -s\hfill \end{array}\right)\mathrm{,}$ (6)

where

$\begin{array}{l}{s}_0=\frac{1-2u^{\mathrm{<mo>*</mo>}}-a{u}^{{\mathrm{<mo>*</mo>}}^2}}{{\left(1+a{u}^{\mathrm{<mo>*</mo>}}\right)}^2}-\frac{m{v}^{\mathrm{<mo>*</mo>}}}{{\left(1+b{u}^{\mathrm{<mo>*</mo>}}\right)}^2\left(1+c{v}^{\mathrm{<mo>*</mo>}}\right)}\mathrm{,}\\ \sigma =-\frac{m{u}^{\mathrm{<mo>*</mo>}}}{\left(1+b{u}^{\mathrm{<mo>*</mo>}}\right){\left(1+c{v}^{\mathrm{<mo>*</mo>}}\right)}^2}\mathrm{.}\end{array}$ (7)

The characteristic polynomical is

$P\left(\lambda \right)={\lambda }^2-\Theta \lambda +\Delta \mathrm{,}$ (8)

where $\Theta \mathrm{<mo>:</mo>}=s_0-s$ and

$\Delta :=-s\left(s_0+\frac{\sigma }{d}\right)=-\frac{s-sa{u}^{*}{}^2-2s{u}^{*}}{{\left(1+a{u}^{*}\right)}^2}+\frac{mds{v}^{*}\left(1+c{v}^{*}\right)+ms{u}^{*}\left(1+b{u}^{*}\right)}{d{\left(1+b{u}^{*}\right)}^2{\left(1+c{v}^{*}\right)}^2}.$

Thus, we have the following conclusions [29].

Theorem 2.1. Assume that ${\Delta }^{*}>0$ .

Define

$\begin{array}{l}T=mds{v}^{\mathrm{<mo>*</mo>}}\left(1+c{v}^{\mathrm{<mo>*</mo>}}\right){\left(1+a{u}^{\mathrm{<mo>*</mo>}}\right)}^2\mathrm{<mi>\; </mi><mi>\; </mi>}+ms{u}^{\mathrm{<mo>*</mo>}}\left(1+b{u}^{\mathrm{<mo>*</mo>}}\right){\left(1+a{u}^{\mathrm{<mo>*</mo>}}\right)}^2\mathrm{,}\\ S=d\left(s-sa{u}^{\mathrm{<mo>*</mo>}}{}^2-2s{u}^{\mathrm{<mo>*</mo>}}\right){\left(1+b{u}^{\mathrm{<mo>*</mo>}}\right)}^2{\left(1+c{v}^{\mathrm{<mo>*</mo>}}\right)}^2\end{array}$

(i) If $s>s_0$ and

(H₁) $T>S$ ,

then the positive equilibrium $\left(u^{\mathrm{<mo>*</mo>}}\mathrm{,}{v}^{\mathrm{<mo>*</mo>}}\right)$ is locally asymptotically stable.

(ii) If (H₁) and

(H₂) $\left(1-a{u}^{\mathrm{<mo>*</mo>}}{}^2-2u^{\mathrm{<mo>*</mo>}}\right){\left(1+b{u}^{\mathrm{<mo>*</mo>}}\right)}^2\left(1+c{v}^{\mathrm{<mo>*</mo>}}\right)>m{v}^{\mathrm{<mo>*</mo>}}{\left(1+a{u}^{\mathrm{<mo>*</mo>}}\right)}^2$ ,

are satisfied, then the positive equilibrium $\left(u^{\mathrm{<mo>*</mo>}}\mathrm{,}{v}^{\mathrm{<mo>*</mo>}}\right)$ is unstable when $s<s_0$ .

(iii) If

(H₃) $T<S$ ,

then the positive equilibrium $\left(u^{\mathrm{<mo>*</mo>}}\mathrm{,}{v}^{\mathrm{<mo>*</mo>}}\right)$ is a saddle point.

Remark 2.2. If the case (i) of Theorem 2.1 holds, $s_0$ is permitted to be positive or negative. If $s_0<0$ , that is,

$\left(1-a{u}^{\mathrm{<mo>*</mo>}}{}^2-2u^{\mathrm{<mo>*</mo>}}\right){\left(1+b{u}^{\mathrm{<mo>*</mo>}}\right)}^2\left(1+c{v}^{\mathrm{<mo>*</mo>}}\right)<m{v}^{\mathrm{<mo>*</mo>}}{\left(1+a{u}^{\mathrm{<mo>*</mo>}}\right)}^2$ (9)

then we have

$S<mds{v}^{*}{\left(1+a{u}^{*}\right)}^2\left(1+c{v}^{*}\right)<T,$

which indicates that (H₁) holds. $s_0$ must be positive when case (ii) in Theorem 2.1 holds. Moreover, we can see that conditions (H₁) and (H₂) can hold simultaneously. Furthermore, from (9) (H₁) and (H₃), we know that $s_0$ is positive in Theorem 2.1 case (iii).

2.2. Hopf Bifurcation

In the following, we analyze the existence of Hopf bifurcation at the interior equilibrium $E^{\mathrm{<mo>*</mo>}}$ (${\Delta }^{*}>0$ ) by choosing c as the bifurcation parameter [30] [31]. In fact, s can be regarded as the intrinsic growth rate of predators and plays an important role in determining the stability of the interior equilibrium and the existence of Hopf bifurcation. Denote

$c^{\mathrm{<mo>*</mo>}}=\frac{m{\left(1+a{u}^{\mathrm{<mo>*</mo>}}\right)}^2}{\left[1-2u^{\mathrm{<mo>*</mo>}}-a{u}^{\mathrm{<mo>*</mo>}}-s{\left(1+a{u}^{\mathrm{<mo>*</mo>}}\right)}^2\right]{\left(1+b{u}^{\mathrm{<mo>*</mo>}}\right)}^2}-\frac{1}{v^{\mathrm{<mo>*</mo>}}}\mathrm{.}$

From the above discussions, if $c=c^{\mathrm{<mo>*</mo>}}$ , then $\text{tr}\left[J\left(E^{\mathrm{<mo>*</mo>}}\right)\right]=0$ , which together with $\det \left[J\left(E^{*}\right)\right]>0$ yield that $J\left(E^{\mathrm{<mo>*</mo>}}\right)$ has a pair of purely imaginary eigenvalues ${\lambda }_{\mathrm{1,2}}=\pm i\sqrt{\text{det}\left[J\left(E^{\mathrm{<mo>*</mo>}}\right)\right]}$ .

We claim that the transeversal condition is satisfied. In fact, if ${\lambda }_{\mathrm{1,2}}=\kappa \left(c\right)+i\omega \left(c\right)$ are the roots of $P\left(\lambda \right)=0$ , then $\kappa \left(c\right)=\frac{1}{2}\text{tr}\left[J\left(E^{\mathrm{<mo>*</mo>}}\right)\right]$ , $\omega \left(c\right)=\frac{1}{2}\sqrt{4\det \left[J\left(E^{*}\right)\right]-\text{tr}{\left[J\left(E^{*}\right)\right]}^2}$ . Hence,

$\kappa \left(c^{*}\right)=0,\mathrm{<mi>\; </mi>}{\kappa }^{\prime }\left(c^{*}\right)=\frac{m{v}^{*}{}^2}{2{\left(1+b{u}^{*}\right)}^2{\left(1+c{v}^{*}\right)}^2}>0.$

This shows that the transversality condition holds. Thus (2) undergoes a Hopf bifurcation about $\left(u^{\mathrm{<mo>*</mo>}}\mathrm{,}{v}^{\mathrm{<mo>*</mo>}}\right)$ as c passes through the $c^{\mathrm{<mo>*</mo>}}$ .

To understand the detailed behaviour of model (2) around $c=c^{\mathrm{<mo>*</mo>}}$ , we need a further analysis of the normal form. We translate the interior equilibrium $\left(u^{\mathrm{<mo>*</mo>}}\mathrm{,}{v}^{\mathrm{<mo>*</mo>}}\right)$ to the origin by the transformation $\tilde{u}=u-u^{\mathrm{<mo>*</mo>}}$ , $\tilde{v}=v-v^{\mathrm{<mo>*</mo>}}$ . For the sake of convenience, we still denote $\tilde{u}$ and $\tilde{v}$ by u and v, respectively. Thus, the local system (2) is transformed into

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

Expanding model (10) in power series around the origin produces the following model

$\left(\begin{array}{l}\frac{\text{d}u}{\text{d}t}\hfill \\ \frac{\text{d}v}{\text{d}t}\hfill \end{array}\right)=J\left(E\right)\left(\begin{array}{l}u\hfill \\ v\hfill \end{array}\right)+\left(\begin{array}{l}f\left(u\mathrm{,}v\mathrm{,}c\right)\hfill \\ g\left(u\mathrm{,}v\mathrm{,}c\right)\hfill \end{array}\right)\mathrm{,}$ (11)

where

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

and

$\begin{array}{l}{a}_1=\frac{mb{v}^{*}}{{\left(1+b{u}^{*}\right)}^3\left(1+c{v}^{*}\right)}-\frac{1+a}{{\left(1+a{u}^{*}\right)}^3},a_2=-\frac{m}{{\left(1+b{u}^{*}\right)}^2{\left(1+c{v}^{*}\right)}^2},\\ a_3=\frac{a+a^2}{{\left(1+a{u}^{*}\right)}^4}+\frac{m{b}^2{v}^{*}}{{\left(1+b{u}^{*}\right)}^4\left(1+c{v}^{*}\right)},a_4=\frac{mb}{{\left(1+b{u}^{*}\right)}^3{\left(1+c{v}^{*}\right)}^2},\\ b_1=-\frac{s}{d^2{v}^{*}},b_2=\frac{2s}{d{v}^{*}},b_3=-\frac{s}{v^{*}},b_4=\frac{s}{d^3{v}^{*}},b_5=-\frac{2s}{d^2{v}^{{*}^2}}.\end{array}$

Set the matrix

$P\mathrm{<mo>:</mo>}=\left(\begin{array}{ll}N\hfill & 1\hfill \\ M\hfill & 0\hfill \end{array}\right)\mathrm{,}$

where

$M=-\frac{s}{\omega },N=-\frac{s_0+s}{2\omega }.$

It is easy to obtain that

$P^{-1}=P^{-1}JP=\Phi \left(c\right):=\left(\begin{array}{ll}\kappa \left(c\right)\hfill & -\omega \left(c\right)\hfill \\ \omega \left(c\right)\hfill & \mathrm{<mi>\; </mi><mi>\; </mi>}\kappa \left(c\right)\hfill \end{array}\right).$

When $c=c^{\mathrm{<mo>*</mo>}}$ , we have

$M_0:={M|}_{c=c^{*}},N_0:={N|}_{c=c^{*}},{\omega }_0:=\omega \left(c^{*}\right).$ (13)

By the transformation ${\left(u\mathrm{,}v\right)}^{{\rm T}}=P{\left(x\mathrm{,}y\right)}^{{\rm T}}$ , model (11) can be rewritten as

$\begin{array}{c}\left(\begin{array}{l}\frac{\text{d}x}{\text{d}t}\hfill \\ \frac{\text{d}y}{\text{d}t}\hfill \end{array}\right)=P^{-1}J\left(E^{\mathrm{<mo>*</mo>}}\right)P\left(\begin{array}{l}x\hfill \\ y\hfill \end{array}\right)+P^{-1}\left(\begin{array}{l}fP\left(x\mathrm{,}y\mathrm{,}c\right)\hfill \\ gP\left(x\mathrm{,}y\mathrm{,}c\right)\hfill \end{array}\right)\\ =\left(\begin{array}{ll}\kappa \left(c\right)\hfill & -\omega \left(c\right)\hfill \\ \omega \left(c\right)\hfill & \mathrm{<mi>\; </mi><mi>\; </mi>}\kappa \left(c\right)\hfill \end{array}\right)\left(\begin{array}{l}x\hfill \\ y\hfill \end{array}\right)+\left(\begin{array}{l}{f}^1\left(x\mathrm{,}y\mathrm{,}c\right)\hfill \\ g^1\left(x\mathrm{,}y\mathrm{,}c\right)\hfill \end{array}\right)\mathrm{,}\end{array}$ (14)

where

$\begin{array}{c}{f}^1\left(x,y,c\right)=\frac{1}{M}g\left(Nx+y,Mx,c\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\\ \mathrm{<mi>\; </mi><mi>\; </mi>}+\left(\frac{N^3}{M}{b}_4+N^2{b}_5\right)x^3+\left(\frac{3N^2}{M}{b}_4+2N{b}_5\right)x^{2}y\\ \mathrm{<mi>\; </mi><mi>\; </mi>}+\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}^1\left(x,y,c\right)=f\left(Nx+y,Mx,c\right)-\frac{N}{M}g\left(Nx+y,Mx,c\right)\\ =\left(N^2{a}_1+NM{a}_2+\frac{N^3}{M}{b}_1-N^2{b}_2-NM{b}_3\right)x^2\\ \mathrm{<mi>\; </mi><mi>\; </mi>}+\left(2N{a}_1+M{a}_2-\frac{2N^2}{M}{b}_1-N{b}_2\right)xy+\left(a_1-\frac{N}{M}{b}_1\right)y^2\\ \mathrm{<mi>\; </mi><mi>\; </mi>}+\left(N^3{a}_3+N^{2}M{a}_4-\frac{N^4}{M}{b}_4-N^3{b}_5\right)x^3\\ \mathrm{<mi>\; </mi><mi>\; </mi>}+\left(3N^2{a}_3+2NM{a}_4-\frac{3N^3}{M}{b}_4-2N^2{b}_5\right)x^{2}y\\ \mathrm{<mi>\; </mi><mi>\; </mi>}+\left(3N{a}_3+M{a}_4-\frac{3N^2}{M}{b}_4-N{b}_5\right)x{y}^2+\left(a_3-\frac{N}{M}{b}_4\right)y^3+\cdots ,\end{array}$

Rewrite (14) in the following polar coordinates form

$\{\begin{array}{l}\dot{r}=\kappa \left(c\right)r+a\left(c\right)r^3+\cdots \mathrm{,}\\ \dot{\theta }=\omega \left(c\right)+c\left(c\right)r^2+\cdots \mathrm{,}\end{array}$ (15)

then the Taylor expansion of (15) at $\delta =s_0$ yields

$\{\begin{array}{l}\dot{r}={\kappa }^{\prime }\left(c^{*}\right)\left(\delta -c^{*}\right)r+a\left(c^{*}\right)r^3+o\left({\left(c-c^{*}\right)}^{2}r,\left(c-c^{*}\right)r^3,r^5\right),\\ \dot{\theta }=\omega \left(c^{*}\right)+{\omega }^{\prime }\left(c^{*}\right)\left(c-c^{*}\right)+c\left(c^{*}\right)r^2+o\left({\left(c-c^{*}\right)}^2,\left(c-c^{*}\right)r^2,r^4\right).\end{array}$ (16)

In order to determine the stability of the Hopf bifurcation periodic solution, we need to calculate the sign of the coefficient $a\left(s_0\right)$ , which is given by

$\begin{array}{l}a\left(c^{*}\right):=\frac{1}{16}\left(f_{uuu}^1+f_{uvv}^1+g_{uuv}^1+g_{vvv}^1\right)\\ \mathrm{<mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi>}+\frac{1}{16{\omega }_0}\left[f_{uv}^1\left(f_{uu}^1+f_{vv}^1\right)-g_{uv}^1\left(g_{uu}^1+g_{vv}^1\right)-f_{uu}^1{g}_{uu}^1+f_{vv}^1{g}_{vv}^1\right],\end{array}$ (17)

where all the partial derivatives are evaluated at the bifurcation point $\left(u\mathrm{,}v\mathrm{,}s\right)=\left(\mathrm{0,0,}{c}^{\mathrm{<mo>*</mo>}}\right)$ , and

$\begin{array}{l}{f}_{uuu}^1\left(0,0,c^{*}\right)=6\left(\frac{N_0^3}{M_0}{b}_4+N_0^2{b}_5\right),f_{uvv}^1\left(0,0,c^{*}\right)=2\left(\frac{3N_0}{M_0}{b}_4+b_5\right),\\ g_{uuv}^1\left(0,0,c^{*}\right)=2\left(3N_0^2{a}_4+2N_0{M}_0{a}_5-\frac{3N_0^2}{M_0}{b}_4-2N_0^2{b}_5\right),\\ g_{vvv}^1\left(0,0,c^{*}\right)=6\left(a_3-\frac{N_0}{M_0}{b}_4\right),f_{uu}^1\left(0,0,c^{*}\right)=2\left(\frac{N_0^2}{M_0}{b}_1+N_0{b}_2+M_0{b}_3\right),\\ f_{uv}^1\left(0,0,c^{*}\right)=\frac{2N_0}{M_0}{b}_1+b_2,f_{vv}^1\left(0,0,c^{*}\right)=\frac{2}{M_0}{b}_1,g_{vv}^1\left(0,0,c^{*}\right)=2\left(a_1-\frac{N_0}{M_0}{b}_1\right),\\ g_{uu}^1\left(0,0,c^{*}\right)=2\left(N_0^2{a}_1+N_0{M}_0{a}_2+M_0^2{a}_3-\frac{N_0^2}{M_0}{b}_1-N_0^2{b}_2-N_0{M}_0{b}_3\right),\\ g_{uv}^1\left(0,0,c^{*}\right)=2N_0{a}_1+M_0{a}_2-\frac{2N_0^2}{M_0}{b}_1-N_0{b}_2.\end{array}$

Thus, we can determine the value and sign of $a\left(c^{*}\right)$ in (17).

Recall that

${\mu }_{*}=-\frac{a\left(c^{*}\right)}{{\kappa }^{\prime }\left(c^{*}\right)},\mathrm{<mi>\; </mi><mi>\; </mi>}\text{and}\mathrm{<mi>\; </mi><mi>\; </mi>}{\kappa }^{\prime }\left(c^{*}\right)>0,$

from the Poincaré-Andronov-Hopf Bifurcation theorem, we can summarize our results as following.

Theorem 2.3. Assume ${\Delta }^{\mathrm{<mo>*</mo>}}>0$ hold. Then system (2) undergoes a Hopf bifurcation at the interior equilibrium $\left(u^{\mathrm{<mo>*</mo>}}\mathrm{,}{v}^{\mathrm{<mo>*</mo>}}\right)$ when $c=c^{\mathrm{<mo>*</mo>}}$ . Furthermore,

(i) $a\left(c^{*}\right)$ determines the stabilities of the bifurcated periodic solutions: if $a\left(c^{\mathrm{<mo>*</mo>}}\right)<0$ (>0); then the bifurcating periodic solutions are stable(unstable);

(ii) ${\mu }_{\mathrm{<mo>*</mo>}}$ determines the directions of Hopf bifurcation: if ${\mu }_{\mathrm{<mo>*</mo>}}>0$ (<0) then the Hopf bifurcation is supercritical (subcritical).

2.3. Transcritical Bifurcation Analysis

Theorem 2.4. The prey-free equilibrium $E_2$ will experience a transcritical bifurcation around $c=c_t=1-\frac{e}{d}$ .

Proof. If $c=c_t$ is chosen as the bifurcation parameter, then system (2) has a zero eigenvalue and a negative eigenvalue $E_2$ , the Jacobian matrix becomes

$J\left(E_2\right)=\left(\begin{array}{ll}0\hfill & \mathrm{<mi>\; </mi><mi>\; </mi>0}\hfill \\ \frac{s}{d}\hfill & -s\hfill \end{array}\right)\mathrm{,}$

In this case, the eigenvectors of matrices $J\left(E_2\right)$ and $J^{{\rm T}}\left(E_2\right)$ are respectively:

$V=\left(\begin{array}{l}{v}_1\hfill \\ v_2\hfill \end{array}\right)=\left(\begin{array}{l}1\hfill \\ \frac{1}{d}\hfill \end{array}\right),\mathrm{<mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi>}W=\left(\begin{array}{l}{w}_1\hfill \\ w_2\hfill \end{array}\right)=\left(\begin{array}{l}1\hfill \\ 0\hfill \end{array}\right).$

Let $F\left(u\mathrm{,}v\right)={\left(f\left(u\mathrm{,}v\right)\mathrm{,}g\left(u\mathrm{,}v\right)\right)}^{{\rm T}}$ , then through calculation we get

$F_c\left(E_2\mathrm{<mo>;</mo>}{c}_t\right)={\left(\begin{array}{l}\frac{\partial f}{\partial c}\hfill \\ \frac{\partial g}{\partial c}\hfill \end{array}\right)}_{\left(E_2\mathrm{<mo>;</mo>}{c}_t\right)}=\left(\begin{array}{l}0\hfill \\ 0\hfill \end{array}\right)\mathrm{,}$

$D{F}_c\left(E_2\mathrm{<mo>;</mo>}{c}_t\right)V=\left(\begin{array}{ll}\frac{\partial f_c}{\partial u}\hfill & \mathrm{<mi>\; </mi><mi>\; </mi>}\frac{\partial f_c}{\partial v}\hfill \\ \frac{\partial g_c}{\partial u}\hfill & \mathrm{<mi>\; </mi><mi>\; </mi>}\frac{\partial g_c}{\partial v}\hfill \end{array}\right){\left(\begin{array}{l}{v}_1\hfill \\ v_2\hfill \end{array}\right)}_{\left(E_2\mathrm{<mo>;</mo>}{c}_t\right)}=\left(\begin{array}{l}\frac{ed-d^2}{e^2}\hfill \\ \mathrm{<mi>\; </mi><mi>\; </mi><mi>\; </mi>0}\hfill \end{array}\right)\mathrm{,}$

$\begin{array}{c}{D}^{2}F\left(E_2\mathrm{<mo>;</mo>}{c}_t\right)\left(V\mathrm{,}V\right)={\left(\begin{array}{l}\frac{{\partial }^{2}f}{\partial u^2}{v}_1{v}_1+2\frac{{\partial }^{2}f}{\partial u\partial v}{v}_1{v}_2+\frac{{\partial }^{2}f}{\partial v^2}{v}_2{v}_2\hfill \\ \frac{{\partial }^{2}g}{\partial u^2}{v}_1{v}_1+2\frac{{\partial }^{2}g}{\partial u\partial v}{v}_1{v}_2+\frac{{\partial }^{2}g}{\partial v^2}{v}_2{v}_2\hfill \end{array}\right)}_{\left(E_2\mathrm{<mo>;</mo>}{c}_t\right)}\\ \mathrm{<mo>=</mo>}\left(\begin{array}{l}-\left(1+a\right)+\frac{db{e}^2-2b^{2}d{e}^3+de-d}{e^4}\hfill \\ \mathrm{<mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi>}\frac{-2sd+2e-se}{d{e}^2}\hfill \end{array}\right)\mathrm{,}\end{array}$

and further there is $W^{{\rm T}}{F}_c\left(E_2\mathrm{<mo>;</mo>}{c}_t\right)=0$ , it means that system (2) doesn’t have a saddle-node bifurcation near $E_2$ . Also

$W^{{\rm T}}\left[D{F}_c\left(E_2\mathrm{<mo>;</mo>}{c}_t\right)V\right]=\frac{ed-d^2}{e^2}\ne \mathrm{0,}$

$W^{{\rm T}}\left[D^{2}F\left(E_2\mathrm{<mo>;</mo>}{c}_t\right)\left(V\mathrm{,}V\right)\right]=-\left(1+a\right)+\frac{db{e}^2-2b^{2}d{e}^3+de-d}{e^4}\ne 0.$

Hence, system (2) will produce a transcritical bifurcation at $c=c_t$ through the Sotomayor’s Theorem.

2.4. Numerical Simulations

In this subsection, we provide numerical simulations to support the analytical results obtained above. The ODE model (2) have seven parameters: a, b, c, d, e, s, m. We illustrate these results by fixing parameters a, b, d, e, s and m and taking the magnitude of interference among predators c as the control parameter.

Example 2.5. We choose parameters as follows

$a=0.2,b=3,d=1.2,e=0.2,m=5,s=\mathrm{0.0555.}$ (18)

Figure 1 shows the phase portraits of model (2) with parameters in (18). In this case, the model (2) has four equilibria: Two saddle points $E_1=\left(\mathrm{1,0}\right)$ , $E_2=\left(\mathrm{0,0.2}\right)$ , a nodal source point $E_0=\left(\mathrm{0,0}\right)$ , a unique coexistence point

$E^{\mathrm{<mo>*</mo>}}=\left(u^{\mathrm{<mo>*</mo>}}\mathrm{,}{v}^{\mathrm{<mo>*</mo>}}\right)$ . By simple computation, we can obtain that $c^{\mathrm{<mo>*</mo>}}\approx 1.169$ . In Figure 1(a), $c=1<c^{*}$ , $E^{\mathrm{<mo>*</mo>}}=\left(\mathrm{0.18929,0.32441}\right)$ is an unstable spiral source and the model exhibits a limit cycle. In Figure 1(b), $c=1.5>c^{*}$ , $E^{\mathrm{<mo>*</mo>}}=\left(\mathrm{0.2829,0.40242}\right)$ is a local asymptotically stable spiral sink. Moreover, when c passes through $c^{\mathrm{<mo>*</mo>}}$ from the left-hand side of $c^{\mathrm{<mo>*</mo>}}$ , $E^{\mathrm{<mo>*</mo>}}$ will lose its stability and Hopf bifurcation occurs, that is, a family of periodic solutions bifurcate from the positive equilibrium. Since $a\left(c^{*}\right)\approx 0.28698>0$ , it follows from Theorem 2.3 that the Hopf

Figure 1. Phase portraits of model (2) with parameters in (18). (a) $c=0.2$ ; (b) $c=1$ .

bifurcation is supercritical and the bifurcation periodic solutions are orbitally asymptotically stable.

3. Turing Instability and Bifurcation of the Reaction-Diffusion Model (19)

In this section, we investigate the reaction-diffusion model (3) to derive the Turing unstable parameter region of the positive equilibria, and the existences, direction and stability of the Hopf bifurcation periodic solutions which describes the spatiotemporal pattern formation.

3.1. Turing Instability of the Positive Equilibria

It is well-known that the equilibrium $\left(u^{\mathrm{<mo>*</mo>}}\mathrm{,}{v}^{\mathrm{<mo>*</mo>}}\right)$ is Turing unstable if it is stable equilibrium of the ODE model (1) but is unstable for the PDE model (3) [32] [33] [34] [35].

For simplicity, we take the spatial domain Ω as the one-dimensional interval $\Omega =\left(\mathrm{0,}\pi \right)$ , and study the following model

$\{\begin{array}{ll}\frac{\partial u}{\partial t}=d_1\Delta u+\frac{u\left(1-u\right)}{1+au}-\frac{muv}{\left(1+bu\right)\left(1+cv\right)},\hfill & x\in \Omega ,t>0,\hfill \\ \frac{\partial v}{\partial t}=d_2\Delta v+sv\left(1-\frac{dv}{u+e}\right),\hfill & x\in \Omega ,t>0,\hfill \\ {\partial }_{v}u={\partial }_{v}v=0,\hfill & x=0,\pi ,t>0,\hfill \\ u\left(x,0\right)=u_0\left(x\right)>0,v\left(x,0\right)=v_0\left(x\right)>0,\hfill & x\in \left(0,\pi \right).\hfill \end{array}$ (19)

Notice that the operator $\phi \mapsto -{\phi }^{″}$ with no-flux boundary conditions has the following eigenvalues and corresponding eigenfunctions

${\mu }_0=\mathrm{0,}{\phi }_0=\sqrt{\frac{1}{\pi }}\mathrm{,}{\mu }_k=k^2\mathrm{,}{\phi }_k\left(x\right)=\sqrt{\frac{2}{\pi }}\cos \left(kx\right)\mathrm{,}\text{for}k=\mathrm{1,<mi>\; </mi>2,<mi>\; </mi>3,}\cdots \mathrm{.}$

The linearized system of (19) at $E^{\mathrm{<mo>*</mo>}}=\left(n^{\mathrm{<mo>*</mo>}}\mathrm{,}{p}^{\mathrm{<mo>*</mo>}}\right)$ is

$\left(\begin{array}{l}\frac{\text{d}u}{\text{d}t}\hfill \\ \frac{\text{d}v}{\text{d}t}\hfill \end{array}\right)=M\left(\begin{array}{l}n\hfill \\ p\hfill \end{array}\right):=N\left(\begin{array}{l}\Delta u\hfill \\ \Delta v\hfill \end{array}\right)+J\left(E^{*}\right)\left(\begin{array}{l}u\hfill \\ v\hfill \end{array}\right),$ (20)

where $J\left(E^{\mathrm{<mo>*</mo>}}\right)$ is the Jacobian matrix defined in (6) and $N=\text{diag}\left(d_1\mathrm{,}{d}_2\right)$ . M is a linear operator with domain

$D_M=U_{ℂ}\mathrm{<mo>:</mo>}=U\oplus iU=\left\{u_1+i{u}_2\mathrm{<mo>|</mo>}{u}_1\mathrm{,}{u}_2\in U\right\}\mathrm{,}$

where

$X:=\left\{{\left(u,v\right)}^{{\rm T}}\in H^2\left[\left(0,\pi \right)\right]\times H^2\left[\left(0,\pi \right)\right]|u_x\left(0,t\right)=u_x\left(\pi ,t\right)=v_x\left(0,t\right)=v_x\left(\pi ,t\right)=0\right\}$

is a real-valued Sobolev space.

According to the standard linear operator theory, if all eigenvalues of the operator have negative real parts, then $\left(u^{\mathrm{<mo>*</mo>}}\mathrm{,}{v}^{\mathrm{<mo>*</mo>}}\right)$ is asymptotically stable. If some eigenvalues have positive real parts, then $\left(u^{\mathrm{<mo>*</mo>}}\mathrm{,}{v}^{\mathrm{<mo>*</mo>}}\right)$ is unstable. Define

$J_k\mathrm{<mo>:</mo>}=J\left(E^{\mathrm{<mo>*</mo>}}\right)-k^{2}N=\left(\begin{array}{ll}{s}_0-k^2{d}_1\hfill & \mathrm{<mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi>}\sigma \hfill \\ \mathrm{<mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi>}\frac{s}{d}\hfill & -s-k^2{d}_2\hfill \end{array}\right)\mathrm{.}$

It is clear that the eigenvalues of the operator M are given by the eigenvalues of the matrix $J_k$ . The characteristic equation of $J_k$ is

${\lambda }^2-T_k\lambda +D_k=0,\mathrm{<mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi>}k=0,1,2,\cdots ,$ (21)

where

$\begin{array}{l}{T}_k\mathrm{<mo>:</mo>}=\text{tr}\left(J_k\right)=\text{tr}\left[J\left(E_1^{\mathrm{<mo>*</mo>}}\right)\right]-k^2\left(d_1+d_2\right)\mathrm{,}\\ D_k\mathrm{<mo>:</mo>}=\text{det}\left(J_k\right)=d_1{d}_2{k}^4-\left(s_0{d}_2-s{d}_1\right)k^2+\text{det}\left[J\left(E_1^{\mathrm{<mo>*</mo>}}\right)\right]\mathrm{,}\end{array}$

By analyzing the root distribution of the characteristic Equation (21), we can draw the following theorem.

Theorem 3.1. Suppose that ${\Delta }^{*}<0$ holds, the equilibrium point $E_1^{\mathrm{<mo>*</mo>}}$ of the system (2) is asymptotically stable when $s>s_0$ , and the equilibrium point $E_1^{\mathrm{<mo>*</mo>}}$ is locally asymptotically stable in the system (24) if and only if the following assumptions are satisfied

(H₁) $d_1\ge \frac{s_0{d}_2}{s}$ ,

(H₂) $0<d_1<\frac{s_0{d}_2}{s}$ and $0<\frac{d_1}{d_2}<\frac{s{s}_0+2\det \left[J\left(E_1^{*}\right)\right]+2\sqrt{\det \left[J\left(E_1^{*}\right)\right]\left[\det \left[J\left(E_1^{*}\right)\right]+s{s}_0\right]}}{s^2}$ ,

whereas $E_1^{\mathrm{<mo>*</mo>}}$ is unstable with respect to the PDE model (19), that is, Turing instability occurs if

(H₃) $\frac{d_1}{d_2}>\frac{s{s}_0+2\det \left[J\left(E_1^{*}\right)\right]+2\sqrt{\det \left[J\left(E_1^{*}\right)\right]\left[\det \left[J\left(E_1^{*}\right)\right]+s{s}_0\right]}}{s^2}$

Proof. First, it is clear that, $T_{k+1}<T_k$ for $k\ge 0$ from the definition of $T_k$ , and $T_0<0$ . So $T_k<0$ for all $k\ge 0$ . Therefore, the signs of real parts of roots of (21) are determined by the signs of $D_k$ , respectively. For the sake of convenience, define

$D\left(k^2\right)\mathrm{<mo>:</mo>}=D_k=d_1{d}_2{k}^4-\left(s_0{d}_2-s{d}_1\right)k^2+\text{det}\left[J\left(E_1^{\mathrm{<mo>*</mo>}}\right)\right]\mathrm{,}$

which is a quadratic polynomial with respect to $k^2$ . The symmetric axis of graph $\left(k^2\mathrm{,}D\left(k^2\right)\right)$ is $k_{\min }^2=\frac{s_0{d}_2-s{d}_1}{2d_1{d}_2}$ . When $k_{\min }^2\le 0$ , $D\left(k^2\right)>0$ for all $k\ge 0$ since $D_0>0$ . When $k_{\min }^2>0$ , $D\left(k^2\right)$ will take the minimum value at $k^2=k_{\min }^2$ , and

$\underset{k}{\min }D\left(k^2\right)=D\left(k_{\min }^2\right)=\det \left[J\left(E_1^{*}\right)\right]-\frac{{\left(s_0{d}_2-s{d}_1\right)}^2}{4d_1{d}_2}.$

The hypothesis (H₁) implies that $k_{\min }^2\le 0$ , (H₂) implies that $k_{\min }^2>0$ and $\underset{k}{\min }D\left(k^2\right)>0$ . In both cases, $D\left(k^2\right)>0$ for all $k\ge 0$ . So all the roots of (21) will have negative real parts, $E_1^{\mathrm{<mo>*</mo>}}$ is asymptotically stable.

When (H₃) holds, $k_{\min }^2>0$ and $\underset{k}{\min }D\left(k^2\right)<0$ , (21) has at least one root with positive real part, therefore $E_1^{\mathrm{<mo>*</mo>}}$ is unstable. This indicates that the Turing instability occurs.

Theorem 3.2. Assume that ${\Delta }^{\mathrm{<mo>*</mo>}}>0$ and (9) hold. Then the unique positive $E^{\mathrm{<mo>*</mo>}}$ in (19) is locally asymptotically stable.

Proof. Let $N=\text{diag}\left(d_1\mathrm{,}{d}_2\right)$ , $E=\left(u\mathrm{,}v\right)$ , $M=N\Delta +J_E\left(E^{\mathrm{<mo>*</mo>}}\right)$ . Then the linearized system of (19) at $E^{\mathrm{<mo>*</mo>}}$ is

$E_t=ME\mathrm{,}$ (22)

and the eigenvalues of the operator M are the eigenvalues of the matrix $-{\mu }_{k}N+J_E\left(E^{\mathrm{<mo>*</mo>}}\right)$ , $\forall k\ge 1$ .

The characteristic equation of $-{\mu }_{k}N+J_E\left(E^{\mathrm{<mo>*</mo>}}\right)$ is

${\Phi }_k\left(\lambda \right)\triangleq \left|\lambda I-{\mu }_{k}N+J_E\left(E^{*}\right)\right|={\lambda }^2+A_k\lambda +B_k=0,$

where

$A_k={\mu }_k\left(d_1+d_2\right)-\Theta \mathrm{<mi>\; </mi>,<mi>\; </mi>}{B}_k=d_1{d}_2{\mu }_k^2-s_0{d}_2{\mu }_k+s{d}_1{\mu }_k+\Delta \mathrm{,}$

and $\Theta \mathrm{,<mi>\; </mi>}\Delta$ are defined as in Section 2.

If (9) holds, then $s_0<0,\mathrm{<mi>\; </mi>}{A}_k>0$ and $B_k>0$ . The roots ${\lambda }_{k\mathrm{,1}}$ , ${\lambda }_{k\mathrm{,2}}$ of ${\Phi }_k\left(\lambda \right)=0$ all have negative real parts.

We claim that there exists a positive constant $\stackrel{⌣}{\delta }$ such that

$Re\left\{{\lambda }_{k\mathrm{,1}}\right\}\mathrm{,}Re\left\{{\lambda }_{k\mathrm{,2}}\right\}\le -\stackrel{⌣}{\delta }\mathrm{,<mi>\; </mi>}\forall k\ge 1.$ (23)

Let $\lambda ={\mu }_k\varsigma$ , then

${\Phi }_k\left(\lambda \right)={\mu }_k^2{\varsigma }^2+A_k{\mu }_k\varsigma +B_k\triangleq {\overline{\Phi }}_k\left(\varsigma \right).$

and

$\underset{k\to \infty }{lim}\frac{\overline{\Phi }\left(\varsigma \right)}{{\mu }_k^2}={\varsigma }^2+\left(d_1+d_2\right)\varsigma +d_1{d}_2\mathrm{.}$

So, the two roots ${\varsigma }_1\mathrm{,}{\varsigma }_2$ of $\underset{k\to \infty }{lim}\frac{\overline{\Phi }\left(\varsigma \right)}{{\mu }_k^2}=0$ always have negative real parts. Let $\overline{d}=\min \left\{d_1,d_2\right\}$ then ${\varsigma }_1\mathrm{,}{\varsigma }_2\le -{\overline{d}}_1$ . By continuity, there exists $k_0$ such that the two roots ${\varsigma }_{k1}\mathrm{,}{\varsigma }_{k2}$ of ${\Phi }_k\left(\lambda \right)=0$ satisfy $Re\left\{{\lambda }_{k\mathrm{,1}}\right\}\mathrm{,}Re\left\{{\lambda }_{k\mathrm{,2}}\right\}\le -\frac{\overline{d}}{2}$ , $\forall k\ge k_0$ . As a result, $Re\left\{{\lambda }_{k\mathrm{,1}}\right\}\mathrm{,}Re\left\{{\lambda }_{k\mathrm{,2}}\right\}\le -\frac{{\mu }_k\overline{d}}{2}\le -\frac{\overline{d}}{2}$ , for $k\ge k_0$ . Let

$-\eta =\underset{1\le k\le k_0}{\max }\left\{Re\left\{{\lambda }_{k\mathrm{,1}}\right\}\mathrm{,}Re\left\{{\lambda }_{k\mathrm{,2}}\right\}\right\}\mathrm{.}$

and $\stackrel{⌣}{\delta }=\min \left\{\eta ,\frac{\overline{d}}{2}\right\}$ . Thus, we obtain the inequation (23), which completes the proof.

3.2. Hopf Bifurcation Analysis

In this subsection, we seek for the possible Hopf bifurcation points and explore the direction and stability of the bifurcation periodic solutions of model (19) with spatial domain $\left(\mathrm{0,}\pi \right)$ . We only discuss the Hopf bifurcation around $E_1^{\mathrm{<mo>*</mo>}}$ . The Hopf bifurcation around $E_2^{\mathrm{<mo>*</mo>}}$ and $E_3^{\mathrm{<mo>*</mo>}}$ can be similarly considered [36].

We take the transformation $\widehat{n}=n-n_1^{\mathrm{<mo>*</mo>}}$ , $\widehat{p}=p-p_1^{\mathrm{<mo>*</mo>}}$ , and still denote $\left(\widehat{n}\mathrm{,}\widehat{p}\right)$ by $\left(n\mathrm{,}p\right)$ . Then model (19) can be rewritten as

$\{\begin{array}{l}\frac{\text{d}u}{\text{d}t}-d_1\Delta u=\frac{u+u^{*}-{\left(u+u^{*}\right)}^2}{1+a\left(u+u^{*}\right)}-\frac{m\left(u+u^{*}\right)\left(v+v^{*}\right)}{\left(1+b\left(u+u^{*}\right)\right)\left(1+c\left(v+v^{*}\right)\right)},\\ \frac{\text{d}v}{\text{d}t}-d_2\Delta v=s\left(v+v^{*}\right)\left(1-\frac{d\left(v+v^{*}\right)}{\left(u+u^{*}\right)+e}\right),\end{array}$ (24)

Let $M^{\mathrm{<mo>*</mo>}}$ be the conjugate operator of L defined by (20). Then

$M^{\mathrm{<mo>*</mo>}}\left(s\right)\left(\begin{array}{c}u\\ v\end{array}\right)\mathrm{<mo>:</mo>}=N\left(\begin{array}{l}\Delta u\hfill \\ \Delta v\hfill \end{array}\right)+J^{\mathrm{<mo>*</mo>}}\left(s\right)\left(\begin{array}{c}u\\ v\end{array}\right)\mathrm{,}$ (25)

where $J^{\mathrm{<mo>*</mo>}}\left(s\right)=J^{{\rm T}}\left(s\right)$ with the domain $D_{M^{\mathrm{<mo>*</mo>}}}=U_C$ . Let

$q\mathrm{<mo>:</mo>}=\left(\begin{array}{l}{m}_0\hfill \\ n_0\hfill \end{array}\right)\mathrm{<mo>=</mo>}\left(\begin{array}{l}\mathrm{<mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi>1}\hfill \\ -\frac{s_0}{\sigma }+\frac{{\omega }_0}{\sigma }\text{i}\hfill \end{array}\right)\mathrm{,<mi>\; </mi><mi>\; </mi><mi>\; </mi>}{q}^{\mathrm{<mo>*</mo>}}\mathrm{<mo>:</mo>}=\left(\begin{array}{l}{m}_0^{\mathrm{<mo>*</mo>}}\hfill \\ n_0^{\mathrm{<mo>*</mo>}}\hfill \end{array}\right)=\frac{1}{2\pi {\omega }_0}\left(\begin{array}{l}{\omega }_0+s_0\text{i}\hfill \\ \mathrm{<mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi>}\sigma \text{i}\hfill \end{array}\right)\mathrm{,}$

It is easy to verify that $\langle M^{\mathrm{<mo>*</mo>}}m\mathrm{,}n\rangle =\langle m\mathrm{,}Mn\rangle$ for any $m\in D_{M^{\mathrm{<mo>*</mo>}}}$ , $n\in D_M$ , and $M\left(s_0\right)q=\text{i}{\omega }_{0}q$ , $M^{\mathrm{<mo>*</mo>}}\left(s_0\right)q^{\mathrm{<mo>*</mo>}}=-\text{i}{\omega }_0{q}^{\mathrm{<mo>*</mo>}}$ , $\langle q^{\mathrm{<mo>*</mo>}}\mathrm{,}\overline{q}\rangle =0$ , $\langle q^{\mathrm{<mo>*</mo>}}\mathrm{,}q\rangle =1$ , where $\langle m\mathrm{,}n\rangle ={\displaystyle {\int }_0^{\pi }}{\overline{m}}^{{\rm T}}n\text{d}x$ denotes the inner product in $L^2\left[\left(\mathrm{0,}\pi \right)\right]\times L^2\left[\left(\mathrm{0,}\pi \right)\right]$ . Similar to [29], let us decompose

$U=U^C\oplus U^S$

with $U^C=\left\{eq+\overline{e}\overline{q}\mathrm{<mo>:</mo>}e\in ℂ\right\}$ and $U^S=\left\{\omega \in U\mathrm{<mo>:</mo>}\langle q^{\mathrm{<mo>*</mo>}}\mathrm{,}\omega \rangle =0\right\}$ .

For any $\left(u\mathrm{,}v\right)\in U$ , there exist $e\in ℂ$ and $\omega =\left({\omega }_1\mathrm{,}{\omega }_2\right)\in U^S$ such that

${\left(u\mathrm{,}v\right)}^{{\rm T}}=eq+\overline{e}\overline{q}+{\left({\omega }_1\mathrm{,}{\omega }_2\right)}^{{\rm T}}\mathrm{,}e=\langle q^{\mathrm{<mo>*</mo>}}\mathrm{,}{\left(u\mathrm{,}v\right)}^{{\rm T}}\rangle \mathrm{.}$

Thus,

$\{\begin{array}{l}u=e+\overline{e}+{\omega }_1\mathrm{,}\\ v=e\left(-\frac{s_0}{\sigma }+\frac{{\omega }_0}{\sigma }\text{i}\right)+\overline{e}\left(-\frac{s_0}{\sigma }+\frac{{\omega }_0}{\sigma }\text{i}\right)+{\omega }_2\mathrm{.}\end{array}$

The model (19) in $\left(e\mathrm{,}\omega \right)$ coordinates becomes

$\{\begin{array}{l}\frac{\text{d}e}{\text{d}t}=\text{i}{\omega }_{0}e+\langle q^{\mathrm{<mo>*</mo>}}\mathrm{,}h\rangle \mathrm{,}\\ \frac{\text{d}\omega }{\text{d}t}=M\omega +h-\langle q^{\mathrm{<mo>*</mo>}}\mathrm{,}h\rangle q-\langle {\overline{q}}^{\mathrm{<mo>*</mo>}}\mathrm{,}h\rangle \overline{q}\mathrm{,}\end{array}$(26)

with $\tilde{h}={\left(f\mathrm{,}g\right)}^{{\rm T}}$ , where f and g are defined as (11). Some direct calculations give

$\begin{array}{l}\langle q^{\mathrm{<mo>*</mo>}}\mathrm{,}h\rangle =\frac{1}{2{\omega }^{\mathrm{<mo>*</mo>}}}\left[{\omega }_{0}f-\text{i}\left(s_{0}f+\sigma g\right)\right]\mathrm{,}\\ \langle {\overline{q}}^{\mathrm{<mo>*</mo>}}\mathrm{,}h\rangle =\frac{1}{2{\omega }_0}\left[{\omega }_{0}f+\text{i}\left(s_{0}f+\sigma g\right)\right]\mathrm{,}\\ \langle q^{\mathrm{<mo>*</mo>}}\mathrm{,}h\rangle q=\frac{1}{2{\omega }_0}\left(\begin{array}{l}\mathrm{<mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi>}{\omega }_{0}f-\text{i}\left(s_{0}f+\sigma g\right)\hfill \\ {\omega }_{0}g+\text{i}\left(\frac{{\omega }_0^2}{\sigma }f+\frac{s_0^2}{\sigma }f+s_{0}g\right)\hfill \end{array}\right)\mathrm{,}\\ \langle {\overline{q}}^{\mathrm{<mo>*</mo>}}\mathrm{,}h\rangle \overline{q}=\frac{1}{2{\omega }_0}\left(\begin{array}{l}\mathrm{<mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi><mi>\; </mi>}{\omega }_{0}f+\text{i}\left(s_{0}f+\sigma g\right)\hfill \\ {\omega }_{0}g-\text{i}\left(\frac{{\omega }_0^2}{\sigma }f+\frac{s_0^2}{\sigma }f+s_{0}g\right)\hfill \end{array}\right)\mathrm{,}\\ H\left(e\mathrm{,}\overline{e}\mathrm{,}\omega \right)\mathrm{<mo>:</mo>}=h-\langle q^{\mathrm{<mo>*</mo>}}\mathrm{,}h\rangle q-\langle {\overline{q}}^{\mathrm{<mo>*</mo>}}\mathrm{,}h\rangle \overline{q}=\left(\begin{array}{l}0\hfill \\ 0\hfill \end{array}\right)\mathrm{.}\end{array}$

By Appendix A of [35], model (26) possesses a center manifold, one can write $\omega$ in the form

$\omega =\frac{{\omega }_{20}}{2}{e}^2+{\omega }_{11}e\overline{e}+\frac{{\omega }_{02}}{2}{\overline{e}}^2+o\left({\left|e\right|}^3\right)\mathrm{.}$

Thus,

$\{\begin{array}{l}\left(2\text{i}{\omega }_{0}I-M\right){\omega }_{20}=\mathrm{0,}\\ \left(-M\right){\omega }_{11}=\mathrm{0,}\\ {\omega }_{02}={\overline{\omega }}_{20}\mathrm{.}\end{array}$

This implies that ${\omega }_{20}={\omega }_{02}={\omega }_{11}=0$ , so that $\left(e\mathrm{,}\overline{e}\right)$ , the equation becomes,

$\frac{\text{d}e}{\text{d}t}=\text{i}{\omega }_{0}e+\frac{1}{2}{g}_{20}{e}^2+g_{11}e\overline{e}+\frac{1}{2}{g}_{02}{e}^2\overline{e}+o{\left(\left|e\right|\right)}^4\mathrm{,}$

where

$\begin{array}{l}{g}_{20}=\frac{1}{2}\left[E_{20}+2E_{11}M\right],\mathrm{<mi>\; </mi><mi>\; </mi><mi>\; </mi>}{g}_{11}=\frac{1}{2}\left[E_{20}+E_{11}\overline{M}+E_{11}M\right],\\ g_{02}=\frac{1}{2}\left[E_{20}+2E_{11}\overline{M}\right],\mathrm{<mi>\; </mi><mi>\; </mi><mi>\; </mi>}{g}_{21}=\frac{1}{2}\left[E_{30}+E_{21}\overline{M}+2E_{21}M\right],\end{array}$

where

$E_{20}=\frac{{\partial }^{2}f}{\partial u^2}\left(u_1^{*},v_1^{*}\right),\mathrm{<mi>\; </mi><mi>\; </mi>}{E}_{11}=\frac{{\partial }^{2}f}{\partial u\partial v}\left(u_1^{*},v_1^{*}\right),\mathrm{<mi>\; </mi><mi>\; </mi>}{E}_{30}=\frac{{\partial }^{3}f}{\partial u^3}\left(u_1^{*},v_1^{*}\right),\mathrm{<mi>\; </mi><mi>\; </mi>}{E}_{21}=\frac{{\partial }^{3}f}{\partial u^2\partial v}\left(u_1^{*},v_1^{*}\right),$

and

$\begin{array}{l}{g}_{20}=\left(a_1-b_1\right)-a_{2}M,\mathrm{<mi>\; </mi><mi>\; </mi><mi>\; </mi>}{g}_{11}=\left(a_1-b_1\right)-\frac{a_2}{2}\left(M+\overline{M}\right),\\ g_{02}=\left(a_1-b_1\right)-a_2\overline{M},\mathrm{<mi>\; </mi><mi>\; </mi><mi>\; </mi>}{g}_{21}=-3a_4+a_3\overline{M}+a_{3}M.\end{array}$

with

$\begin{array}{l}{a}_1=\frac{mb{v}^{*}}{{\left(1+b{u}^{*}\right)}^3\left(1+c{v}^{*}\right)}-\frac{1+a}{{\left(1+a{u}^{*}\right)}^3},a_2=-\frac{m}{{\left(1+b{u}^{*}\right)}^2{\left(1+c{v}^{*}\right)}^2},\\ a_3=\frac{a+a^2}{{\left(1+a{u}^{*}\right)}^4}+\frac{m{b}^2{v}^{*}}{{\left(1+b{u}^{*}\right)}^4\left(1+c{v}^{*}\right)},a_4=\frac{mb}{{\left(1+b{u}^{*}\right)}^3{\left(1+c{v}^{*}\right)}^2},\end{array}$

$\begin{array}{l}{b}_1=-\frac{s}{d^2{v}^{*}},b_2=\frac{2s}{d{v}^{*}}.\\ M=\frac{{\left(1+b{u}^{*}\right)}^2{\left(1+c{v}^{*}\right)}^2-v^{*}{\left(1+a{u}^{*}\right)}^2\left(1+c{v}^{*}\right)}{m{u}^{*}{\left(1+a{u}^{*}\right)}^2\left(1+b{u}^{*}\right)}-\frac{{\omega }_0\left(1+b{u}^{*}\right)\left(1+c{v}^{*}\right)}{m{u}^{*}}\text{i},\end{array}$

It can conclude from the above calculation

${\eta }_1\left(s_0\right)=\frac{\text{i}}{2{\omega }_0}\left(g_{20}{g}_{11}-2{\left|g_{11}\right|}^2-\frac{1}{3}{g}_{02}^2+\frac{1}{2}{g}_{21}\right)\mathrm{,}$

then

$\begin{array}{c}\text{Re}\left({\eta }_1\left(s_0\right)\right)\mathrm{<mo>:</mo>}=\text{Re}\left\{\frac{i}{2{\omega }_0}\left(g_{20}{g}_{11}-2{\left|g_{11}\right|}^2-\frac{1}{3}{\left|g_{02}\right|}^2\right)+\frac{1}{2}{g}_{21}\right\}\\ =-\frac{1}{2{\omega }_0}\left[\text{Re}\left(g_{20}\right)\text{Im}\left(g_{11}\right)+\text{Im}\left(g_{20}\right)\text{Re}\left(g_{11}\right)\right]+\frac{1}{2}\text{Re}\left(g_{21}\right)\end{array}$ (27)

which together with ${\mu }_2=\frac{\text{Re}\left({\eta }_1\left(s_0\right)\right)}{{\kappa }^{\prime }\left(s_0\right)}$ determine the stability and direction of the bifurcated periodic solutions.

Theorem 3.3. Assume ${\Delta }^{*}<0$ hold. Then model (19) undergoes a Hopf bifurcation at $s=s_0$ .

(i) If $\text{Re}\left(a\left(c^{\mathrm{<mo>*</mo>}}\right)\right)<0$ , then the Hopf bifurcation is supercritical and the bifurcating periodic solutions are asymptotically stable on the centre manifold. Furthermore, they are orbitally asymptotically stable for model (19) if (H₁) or (H₂) holds, and unstable if (H₃) holds.

(ii) If $\text{Re}\left(a\left(c^{\mathrm{<mo>*</mo>}}\right)\right)>0$ , then the Hopf bifurcation is subcritical and the bifurcating periodic solutions are unstable.

3.3. Numerical Simulations

In this subsection, we give some numerical simulations to illustrate our theoretical analysis, we consider model (19) in one-dimensional space.

Example 3.4. We choose parameters as follows

$a=0.2,\mathrm{<mi>\; </mi>}b=3,\mathrm{<mi>\; </mi>}c=1.5,\mathrm{<mi>\; </mi>}d=1.2,\mathrm{<mi>\; </mi>}e=0.2,\mathrm{<mi>\; </mi>}m=5,\mathrm{<mi>\; </mi>}s=0.0555,\mathrm{<mi>\; </mi>}{d}_1=\mathrm{0.01.}$ (28)

Then $c<c^{*}\approx 1.169$ , and the unique positive equilibrium $E^{\mathrm{<mo>*</mo>}}=\left(\mathrm{0.2829,0.40242}\right)$ in model (2) is locally asymptotically stable. If we choose $d_2=0.1$ , then (H₂) holds, by Theorem 3.1, the equilibrium $E^{\mathrm{<mo>*</mo>}}$ in model (19) is still locally asymptotically stable (see Figure 2). If we choose $d_2=0.5$ , then (H₃) holds, by Theorem 3.1, the equilibrium $E^{\mathrm{<mo>*</mo>}}$ in model (19) becomes unstable, this means that the Turing instability of the equilibrium solution happens (see Figure 3).

Example 3.5. We choose parameters as follows

$a=0.2,\mathrm{<mi>\; </mi>}b=3,\mathrm{<mi>\; </mi>}c=1,\mathrm{<mi>\; </mi>}d=1.2,\mathrm{<mi>\; </mi>}e=0.2,\mathrm{<mi>\; </mi>}m=5,\mathrm{<mi>\; </mi>}s=0.0555,\mathrm{<mi>\; </mi>}{d}_1=\mathrm{0.01.}$ (29)

In this case, we have $c^{\mathrm{<mo>*</mo>}}\approx 1.169$ , $\mathrm{Re}\left(a\left(c^{*}\right)\right)\approx 0.28298>0$ . By Theorem 3.3, the supercritical Hopf bifurcation occurs at $c=c^{\mathrm{<mo>*</mo>}}$ . Notice that $c=1<c^{*}$ . The unique positive equilibrium $E^{\mathrm{<mo>*</mo>}}=\left(\mathrm{0.18929,0.32441}\right)$ in model (19) is unstable and bifurcating periodic orbits exist. If we choose $d_2=0.5$ , then (H₂) holds, the bifurcating periodic orbits are orbitally asymptotically stable (see Figure 4). If we choose $d_2=0.06$ , then (H₃) holds, so the bifurcating periodic orbits are

Figure 2. Numerical simulations of the stable equilibrium of model (19) with parameters in (28) and $d_2=0.1$ .

Figure 3. Numerical simulations of the Turing instability of the equilibrium of model (19) with parameters in (28) and $d_2=0.5$ .

Figure 4. Numerical simulations of the stable bifurcating periodic solution of model (19) with parameters in (29) and $d_2=0.1$ .

Figure 5. Numerical simulations of the Turing instability of bifurcating periodic solution of model (19) with parameters in (29) and $d_2=0.5$ .

Turing unstable (see Figure 5).

4. Conclusion and Discussion

The pattern formation of ecosystem has always been an important and fundamental topic in ecology. In this paper, we consider a diffused modified Leslie-Gower predator-prey system with a C-M functional response under homogeneous Neumann boundary conditions. Firstly, the local asymptotic stability and bifurcation in corresponding ODE systems are studied. (i) At the equilibrium point $E^{\mathrm{<mo>*</mo>}}$ , Hopf bifurcation occurs when $c=c^{\mathrm{<mo>*</mo>}}$ . (ii) At the equilibrium

point $E_2$ , transcritical bifurcation occurs when $c=1-\frac{e}{d}$ . Secondly, we consider

the Turing (diffusion-driven) instability of reaction-diffusion systems in co-equilibrium when the spatial domain is a bounded interval, which produces a spatially inhomogeneous pattern. Besides, we investigate the existence and direction of Hopf bifurcations and the stability of periodic solutions of bifurcations in a reaction-diffusion system, which exhibits a time-periodic pattern. Finally, we discuss the interaction of Turing instability and Hopf bifurcation in a reaction-diffusion system that exhibits a spatio-temporal pattern. Our theoretical results further suggest that predator-prey interference and predator feeding strategies are determinants of spatial and spatio-temporal patterns generated through predator-prey interactions in a uniform environment.

Since the C-M functional response function is more general than the B-D functional response function and contains the Holling II functional response function, we would like to point out that our results are still valid for the diffuse Leslie-Gower predator-prey system and the diffuse Holling-tanner predator-prey system with the B-D functional response.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	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
[2]	Yin, H., Xiao, X., Wen, X. and Liu, K. (2014) Pattern Analysis of a Modified Leslie-Gower Predator-Prey Model with Crowley-Martin Functional Response and Diffusion. Computers & Mathematics with Applications, 67, 1607-1621. https://doi.org/10.1016/j.camwa.2014.02.016
[3]	Leslie, P.H. and Gower, J.C. (1960) The Properties of a Stochastic Model for the Predator-Prey Type of Interaction between Two Species. Biometrika, 47, 219-234. https://doi.org/10.1093/biomet/47.3-4.219
[4]	Gierer, A. and Meinhardt, H. (1972) A Theory of Biological Pattern Formation. Kybernetik, 12, 30-39.
[5]	He, M. and Li, Z. (2023) Global Dynamics of a Leslie-Gower Predator-Prey Model with Square Root Response Function. Applied Mathematics Letters, 140, Article 108561. https://doi.org/10.1016/j.aml.2022.108561
[6]	Li, Y., He, M. and Li, Z. (2022) Dynamics of a Ratio-Dependent Leslie-Gower Predator-Prey Model with Allee Effect and Fear Effect. Mathematics and Computers in Simulation, 201, 417-439. https://doi.org/10.1016/j.matcom.2022.05.017
[7]	Smith, F.E. (1963) Population Dynamics in Daphnia Magna and a New Model for Population Growth. Ecology, 44, 651-663. https://doi.org/10.2307/1933011
[8]	Yue, Z. and Wang, W. (2013) Qualitative Analysis of a Diffusive Ratio-Dependent Holling-Tanner Predator-Prey Model with Smith Growth. Discrete Dynamics in Nature and Society, 2013, Article 267173. https://doi.org/10.1155/2013/267173
[9]	Misra, O.P. and Babu, A.R. (2014) A Model for the Effect of Toxicant on a Three Species Food-Chain System with “Food-Limited” Growth of Prey Population. Global Journal of Mathematical Analysis, 2, 120-145. https://doi.org/10.14419/gjma.v2i3.2990
[10]	Shi, Z., Cheng, H., Liu, Y. and Li, Y. (2019) A Cydia Pomonella Integrated Management Predator-Prey Model with Smith Growth and Linear Feedback Control. IEEE Access, 7, 126066-126076. https://doi.org/10.1109/access.2019.2938772
[11]	Jiang, H., Fang, H. and Wu, Y. (2020) Hopf Bifurcation in a Diffusive Predator-Prey Model with Smith Growth Rate and Herd Behavior. Advances in Difference Equations, 2020, Article No. 518. https://doi.org/10.1186/s13662-020-02879-4
[12]	Feng, X., Sun, H., Xiao, Y. and Xiao, F. (2020) Stability and Coexistence of a Diffusive Predator-Prey System with Nonmonotonic Functional Response and Fear Effect. Complexity, 2020, Article 8899114. https://doi.org/10.1155/2020/8899114
[13]	Dubey, V.P., Kumar, R. and Kumar, D. (2020) Numerical Solution of Time-Fractional Three-Species Food Chain Model Arising in the Realm of Mathematical Ecology. International Journal of Biomathematics, 13, Article 2050011. https://doi.org/10.1142/s1793524520500114
[14]	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
[15]	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
[16]	Luo, D. and Wang, Q. (2022) Global Bifurcation and Pattern Formation for a Reaction-Diffusion Predator-Prey Model with Prey-Taxis and Double Beddington-Deangelis Functional Responses. Nonlinear Analysis: Real World Applications, 67, Article 103638. https://doi.org/10.1016/j.nonrwa.2022.103638
[17]	Jiang, H., Wang, L. and Yao, R. (2015) Numerical Simulation and Qualitative Analysis for a Predator-Prey Model with B-D Functional Response. Mathematics and Computers in Simulation, 117, 39-53. https://doi.org/10.1016/j.matcom.2015.05.006
[18]	Min, N. and Wang, M. (2018) Dynamics of a Diffusive Prey-Predator System with Strong Allee Effect Growth Rate and a Protection Zone for the Prey. Discrete & Continuous Dynamical Systems-B, 23, 1721-1737. https://doi.org/10.3934/dcdsb.2018073
[19]	Ali, N. and Jazar, M. (2013) Global Dynamics of a Modified Leslie-Gower Predator-Prey Model with Crowley-Martin Functional Responses. Journal of Applied Mathematics and Computing, 43, 271-293. https://doi.org/10.1007/s12190-013-0663-3
[20]	Chen, M., Takeuchi, Y. and Zhang, J. (2023) Dynamic Complexity of a Modified Leslie-Gower Predator-Prey System with Fear Effect. Communications in Nonlinear Science and Numerical Simulation, 119, Article 107109. https://doi.org/10.1016/j.cnsns.2023.107109
[21]	Melese, D. and Feyissa, S. (2021) Stability and Bifurcation Analysis of a Diffusive Modified Leslie-Gower Prey-Predator Model with Prey Infection and Beddington Deangelis Functional Response. Heliyon, 7, E06193. https://doi.org/10.1016/j.heliyon.2021.e06193
[22]	Cosner, C., DeAngelis, D.L., Ault, J.S. and Olson, D.B. (1999) Effects of Spatial Grouping on the Functional Response of Predators. Theoretical Population Biology, 56, 65-75. https://doi.org/10.1006/tpbi.1999.1414
[23]	Wu, S., Shi, J. and Wu, B. (2016) Global Existence of Solutions and Uniform Persistence of a Diffusive Predator-Prey Model with Prey-Taxis. Journal of Differential Equations, 260, 5847-5874. https://doi.org/10.1016/j.jde.2015.12.024
[24]	Han, X. and Lei, C. (2023) Bifurcation and Turing Instability Analysis for a Space-and Time-Discrete Predator-Prey System with Smith Growth Function. Chaos, Solitons & Fractals, 166, Article 112910. https://doi.org/10.1016/j.chaos.2022.112910
[25]	Feng, X., Liu, X., Sun, C. and Jiang, Y. (2023) Stability and Hopf Bifurcation of a Modified Leslie-Gower Predator-Prey Model with Smith Growth Rate and B-D Functional Response. Chaos, Solitons & Fractals, 174, Article 113794. https://doi.org/10.1016/j.chaos.2023.113794
[26]	Okubo, A. and Levin, S.A. (2002) Diffusion and Ecological Problems: Modern Perspectives. Springer. https://doi.org/10.1007/978-1-4757-4978-6
[27]	Huang, J., Ruan, S. and Song, J. (2014) Bifurcations in a Predator-Prey System of Leslie Type with Generalized Holling Type III Functional Response. Journal of Differential Equations, 257, 1721-1752. https://doi.org/10.1016/j.jde.2014.04.024
[28]	Li, Y. and Xiao, D. (2007) Bifurcations of a Predator-Prey System of Holling and Leslie Types. Chaos, Solitons & Fractals, 34, 606-620. https://doi.org/10.1016/j.chaos.2006.03.068
[29]	Shi, H.-B. and Ruan, S. (2015) Spatial, Temporal and Spatiotemporal Patterns of Diffusive Predator-Prey Models with Mutual Interference. IMA Journal of Applied Mathematics, 80, 1534-1568. https://doi.org/10.1093/imamat/hxv006
[30]	Hale, J.K. and Kocak, H. (1991) Dynamics and Bifurcations. Springer. https://doi.org/10.1007/978-1-4612-4426-4
[31]	Takens, F. (1974) Forced Oscillations and Bifurcations. Communications of the Mathematical Institute, Rijksuniversiteit Utrecht, 3, 1-61.
[32]	Li, X., Jiang, W. and Shi, J. (2011) Hopf Bifurcation and Turing Instability in the Reaction-Diffusion Holling-Tanner Predator-Prey Model. IMA Journal of Applied Mathematics, 78, 287-306. https://doi.org/10.1093/imamat/hxr050
[33]	Capone, F., Carfora, M.F., De Luca, R. and Torcicollo, I. (2019) Turing Patterns in a Reaction-Diffusion System Modeling Hunting Cooperation. Mathematics and Computers in Simulation, 165, 172-180. https://doi.org/10.1016/j.matcom.2019.03.010
[34]	Gao, X., Ishag, S., Fu, S., Li, W. and Wang, W. (2020) Bifurcation and Turing Pattern Formation in a Diffusive Ratio-Dependent Predator-Prey Model with Predator Harvesting. Nonlinear Analysis: Real World Applications, 51, Article 102962. https://doi.org/10.1016/j.nonrwa.2019.102962
[35]	Alonso, D., Bartumeus, F. and Catalan, J. (2002) Mutual Interference between Predators Can Give Rise to Turing Spatial Patterns. Ecology, 83, 28-34. https://doi.org/10.1890/0012-9658(2002)083[0028:MIBPCG]2.0.CO;2
[36]	Hassard, B.D., Kazarinoff, N.D. and Wan, Y.H. (1981) Theory and Applications of Hopf Bifurcation. Cambridge University Press.