Discontinuity, Nonlinearity, and Complexity
Stationary Pattern in a Predator-Prey Model with Predator-Harvesting Policy
Discontinuity, Nonlinearity, and Complexity 12(2) (2023) 411--436 | DOI:10.5890/DNC.2023.06.013
Tian Wang$^{1}$, Ranchao Wu$^{1}$, Mengxin Chen$^{2}$
$^1$ School of Mathematical Sciences, Anhui University, Hefei 230601, China
$^2$ College of Mathematics and Information Science, Henan Normal University, Xinxiang 453007, China
Download Full Text PDF
Abstract
Stationary patterns in a predator-prey model with Holling-III functional response and harvesting policy are investigated in this work. For the proposed model, nonnegativity, uniformly boundedness, permanence, stability, the existence and direction of Hopf bifurcation are analyzed. For the spatial system, the existence conditions for the Turing instability are also established. Then using weakly nonlinear analysis, amplitude equations near critical values of the Turing instability are derived. Different kinds of solutions can be shown analytically by analyzing amplitude equations. Based on numerical analysis, the stationary patterns can be found, such as hexagonal patterns, stripe patterns and mixed states of hexagonal pattern and stripe pattern. Numerical simulations are well consistent with theoretical results. It is further noted that the behavior of harvesting is a factor for existence and stability of equilibria, the occurrence of the transcritical bifurcation, pattern formation and the permanence of the system.}[\hfill Predator-prey model\par
\hfill Stationary patterns\par
\hfill Turing instability\par
\hfill Harvesting][S.Y. Xing][16 April 2021][23 July 2021][1 April 2023][2023 L\&H Scientific Publishing, LLC. All rights reserved.]
\maketitle
%\thispagestyle{fancy}
\thispagestyle{firstpage}
\renewcommand{\baselinestretch}{1}
\normalsize
\section{Introduction}
\noindent Interactions between organisms and their surroundings play an important part during ecological and biological
processes. Predation is one of the typical interactions in which one species kills or eats another, since animals need to hunt in order to survive and make possible the continuation of race. A lot of predator-prey models describing this interaction have been proposed and studied. One of the classical predator-prey models is the Leslie-Gower predator-prey model, for details we can refer to various results \cite{1,2,3} and so on. In the species interaction model, the function response is also of great significance in describing population dynamics. Functional response is the intake rate of per predator as a function of food density. There are different kinds of functional responses, for instance, Holling type, Ivlev-type, Hassell-Varley type and so on, see related references \cite{4,5,6,7,8}. Different functional responses will lead to different dynamics and show the complex relationships between species.
There are many interesting phenomena and relationships in species, the inhibition of the prey is generally associated with the foraging behavior, the complex predation may exists, such as chasing behavior. So we choose the Holling type functional responses. The above-mentioned Holling type functional responses are only prey-dependent. However, predators have to share or compete for food when they search for food. We need to establish a functional response that the per capita predator growth rate should be a function of the ratio of prey to predator abundance, that is to say, ratio-dependent functional response \cite{9}. Considering the following Leslie-Gower predator-prey model with ratio-dependent Holling-III functional response
\begin{eqnarray}
\left\{
\begin{split}
&\frac{\partial u}{\partial t}-d_1\Delta u=u(1-u)-\frac{\beta u^2v}{u^2+mv^2},\\
&\frac{\partial v}{\partial t}-d_2\Delta v=r v(1-\frac{v}{u}),
\end{split}\right.
\end{eqnarray}
where $u$ and $v$ represent the densities of the prey and predator population respectively. $r$ denotes the intrinsic growth of the predator, $d_1$ and $d_2$ stand for the positive diffusion coefficients of the prey and predator, $\beta$, $m$ are positive constants. The function $\frac{u}{v}$ is known as the Leslie-Gower term \cite{10}, the term $\frac{\beta u^2v}{u^2+mv^2}$ is called ratio-dependent Holling-III response, it means that the per capita predator growth rate should be a function of the ratio of prey to predator abundance. $\Delta=\frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2}$ stands for the Laplacian operator. Notice that the spatiotemporal dynamics of Eq. (1) has been investigated in \cite{11}, including spatial patterns, temporal patterns and spatiotemporal patterns. Chang et al \cite{12} presented the study of predator-prey system with ratio-dependent functional response and showed the conditions for Hopf bifurcation and its properties. Many researchers have studied the dynamic behaviors of the predator-prey model with functional response, we can refer to \cite{13,14,15}.
As we know, the harvesting is also a factor of the dynamic evolution of population. The benefits of harvesting far outweigh the costs, it protects species population, brings economic gain and provides entertainment. There are two types of harvesting regimes, the constant-yield harvesting and proportional harvesting, however, the well-known nonlinear type harvesting is Michaelis-Menten type harvesting. Harvesting in a ratio-dependent predator-prey model has been considered in some studies. Kar et al \cite{16} investigated the impacts of maximum sustainable yield (MSY) policy on two different predator-prey systems under different harvesting schemes. The results suggest that the harvesting of predator specie at MSY level may be a sustainable policy for two systems. Huang et al \cite{17} studied the dynamical behavior of Leslie-Gower predator-prey model with constant-yield predator harvesting and discussed the effect of harvesting on the number and kind of equilibria and the type of bifurcations. For more related works, see \cite{18,19,20,21,31}.
With the harvesting of predator, the target species harvesting not only promotes economic development, but also protects species diversity. In this paper, we investigate the effect of the linear harvesting on predators population, adding the harvesting term to Eq. (1). The corresponding model is
\begin{eqnarray}
\left\{
\begin{split}
&\frac{\partial u}{\partial t}=d_1\Delta u+u(1-u)-\frac{\beta u^2v}{u^2+mv^2},~&(x,y)\in\Omega,~t>0,\\
&\frac{\partial v}{\partial t}=d_2\Delta v+rv(1-\frac{v}{u})-av,~&(x,y)\in\Omega,~t>0,\\
&\frac{\partial u}{\partial \mathbf{n}}=\frac{\partial v}{\partial \mathbf{n}}=0,~&(x,y)\in\partial\Omega,~t\geq0,\\
&u(x,y,0)=u_0(x,y)\geq0,~v(x,y,0)=v_0(x,y)\geq0,~&(x,y)\in\Omega,
\end{split} \right.
\end{eqnarray}
where $u(x,y,t)$ and $v(x,y,t)$ are the densities of the prey and predator in time $t$ and space $(x,y)\in\Omega$, respectively, $\Omega$ is a bounded domain in $R^2$ with smooth boundary $\partial\Omega$. $\mathbf{n}$ is the outward unit normal vector on $\partial\Omega$. $av$ represents linear harvesting term, where $a$ is harvesting rate, $r$ denotes the intrinsic growth of the predator, $d_1$ and $d_2$ stand for the positive diffusion coefficients of the prey and predator, and $\beta$, $m$ are all positive constants. Now we consider that no external input is imposed from the outside, so take the homogeneous Neumann boundary conditions, which indicates that the predator-prey system is self-contained with zero population flux across the boundary.
One of interesting research topics in population models is the Turing instability \cite{22}. Over the years, the problems of the Turing instability have been widely discussed and applied to different fields, such as biology, chemistry and physics. In \cite{23},
the authors investigated the persistence, Turing instability and bifurcation analysis for the diffusive Leslie-Gower predator-prey model, and patterns formation in one dimension space. Han et al \cite{24} considered the toxic-phytoplankton-zooplankton model with nonmonotonic function response, and showed the different kinds of spatiotemporal patterns, such as spot, stripe, spot-stripe pattern. In the work \cite{25}, the conditions of secondary instability of Gierer-Meinhardt Model were obtained. The patterns in one-dimensional or two-dimensional space were also studied, for more results about Turing instability and induced patterns, see references \cite{26,27,28,29}. The purpose of this paper is to explore the effect of harvesting on the model (2) and investigate the different kinds of Turing patterns in two-dimensional space.
The outline of the paper as follows. In Section 2, we analyze the uniform boundedness and permanence of the solutions, the existence and stability of the equilibria and the direction of Hopf bifurcation and then focus on the investigation of transcritical bifurcation. In Section 3, the existence conditions of the Turing instability will be studied. In Section 4, the amplitude equations at the Turing bifurcation point are established by weakly nonlinear analysis and analyzed. In Section 5, some numerical simulations for the patterns are carried out. Finally, some conclusions are obtained in Section 6.
\section{Stability analysis}
\subsection{Uniform boundedness and permanence}
\noindent
In this part, consider Eq. (2) without diffusion, then it becomes the following system
\begin{eqnarray}
\left\{
\begin{split}
&\frac{d u}{d t}=u(1-u)-\frac{\beta u^2v}{u^2+mv^2},\\
&\frac{d v}{d t}=rv(1-\frac{v}{u})-av.
\end{split} \right.
\end{eqnarray}
\begin{lemma}
Assume that the conditions: $u(0)\geq0$, $v(0)\geq0$ hold, then all the solutions of Eq. (3) are nonnegative and uniformly bounded, for all $t\geq0$.
\end{lemma}
\begin{proof}
From Eq. (3), when $u(0)\geq0$, $v(0)\geq0$, it holds that,
\begin{eqnarray*}
&&u(t)=u(0)e^{\int^{t}_{0}(1-u-\frac{\beta uv}{u^2+mv^2})ds}\geq0,\\
&&v(t)=v(0)e^{\int^{t}_{0}(r-a-r\frac{v}{u})ds}\geq0,
\end{eqnarray*}
so it is not difficult to infer that the solutions are nonnegative, for all $t\geq0$. Next prove that solutions are uniformly bounded. Define a function $\chi(t)=u(t)+v(t)$ and take the derivative of this function with respect to $t$, we can get
\begin{eqnarray*}
\frac{d \chi(t)}{dt}&&=\frac{d u(t)}{d t}+\frac{d v(t)}{d t}\\
&&=u(1-u)-\frac{\beta u^2v}{u^2+mv^2}+rv(1-\frac{v}{u})-av,
\end{eqnarray*}
let $\eta\leq a-r$, we have
\begin{eqnarray*}
\frac{d \chi(t)}{dt}+\eta \chi&&=u(1-u)-\frac{\beta u^2v}{u^2+mv^2}+rv(1-\frac{v}{u})-av+\eta u+\eta v\\
&&=-u^2+u+\eta u+rv-av+\eta v-\frac{\beta u^2v}{u^2+mv^2}-\frac{rv^2}{u}\\
&&\leq-u^2+u+\eta u+(r-a+\eta) v\\
&&\leq u(1+\eta-u)\\
&&\leq\frac{(1+\eta)^2}{4},
\end{eqnarray*}
it can be deduced that $\lim _{t \to \infty}\sup\chi(t)\leq\frac{(1+\eta)^2}{4\eta}$, so the solutions of the system are uniformly bounded.
\end{proof}
\begin{lemma}
(Permanence) Suppose $\beta\leq2\sqrt{m}$, $r-a\geq0$, then solutions of Eq. (3) satisfy the following conditions
\begin{eqnarray*}
&&l_1\leq \lim_{t \to \infty}\inf u(t)\leq \lim_{t \to \infty}\sup u(t)\leq 1,\\
&&l_2\leq \lim_{t \to \infty}\inf v(t)\leq \lim_{t \to \infty}\sup v(t)\leq 1,
\end{eqnarray*}
where $l_1=\frac{2\sqrt m-\beta}{2\sqrt m}$, $l_2=\frac{(r-a)(2\sqrt{m}-\beta)}{2r\sqrt{m}}$.
\end{lemma}
\begin{proof}
For $\beta\leq2\sqrt{m}$, one has
\begin{eqnarray*}
\frac{d u}{d t}&&=u(1-u)-\frac{\beta u^2v}{u^2+mv^2}\\
&&\geq u(1-u)-\frac{\beta u^2v}{2\sqrt{m}uv}\\
&&\geq u(1-\frac{\beta}{2\sqrt{m}}-u),
\end{eqnarray*}
we can get $\lim _{t \to \infty}\inf u(t)\geq\frac{2\sqrt m-\beta}{2\sqrt m}$, and then considering the following
\begin{equation*}
\begin{split}
\frac{d u}{d t}&=u(1-u)-\frac{\beta u^2v}{u^2+mv^2}\\
&\leq u(1-u),
\end{split}
\end{equation*}
we can deduce $\lim _{t \to \infty}\sup u(t)\leq 1$.
For $r-a\geq0$, it is true that
\begin{eqnarray*}
\frac{d v}{d t}&&=rv(1-\frac{v}{u})-av\\
&&=v(r-a-\frac{rv}{u})\\
&&\geq v(r-a-\frac{2r\sqrt{m}v}{2\sqrt{m}-\beta}),
\end{eqnarray*}
we can know $\lim _{t \to \infty}\inf v(t)\geq \frac{(r-a)(2\sqrt{m}-\beta)}{2r\sqrt{m}}$, and
then considering the following form
\begin{equation*}
\begin{split}
\frac{d v}{d t}&=rv(1-\frac{v}{u})-av\\
&=v(r-a-\frac{rv}{u})\\
&\leq v(r-\frac{rv}{u}),
\end{split}
\end{equation*}
we can infer $\lim _{t \to \infty}\sup v(t)\leq\lim _{t \to \infty}\sup u(t)\leq 1$.
Obviously, based on the above results, the system is permanent.
\end{proof}
\newpage
\subsection{Equilibria and their stability}
\noindent
In this subsection, we will explore the existence and the stability of equilibria of Eq. (3). Obviously, Eq. (3) has a boundary equilibrium $E_{1}=(u_1,v_1)=(1,0)$. In order to obtain the positive equilibria $E_*=(u_2,v_2)$, let
\begin{eqnarray}
\left\{
\begin{split}
&f(u,v)\triangleq u(1-u)-\frac{\beta u^2v}{u^2+mv^2}=0,\\
&g(u,v)\triangleq rv(1-\frac{v}{u})-av=0.
\end{split} \right.
\end{eqnarray}
Suppose $b=\frac{r-a}{r}$, where $0\frac{(1+br)(1+mb^2)^2}{2b}$, $E_*$ is unstable.
(iv)~If $\beta=\beta_{H}$, where $\beta_{H}=\frac{(1+br)(1+mb^2)^2}{2b}$, system (3) undergoes the Hopf bifurcation.
\end{theorem}
\begin{proof}
For the equilibrium $E_{1}=(1,0)$, the Jacobian matrix of Eq. (3) at $E_{1}$ is written as
\begin{eqnarray*}
J_{E_{1}}=&\left(
\begin{array}{cc}
-1 & -\beta \\
0 & r-a
\end{array}
\right),
\end{eqnarray*}
then the eigenvalues $\lambda_{1}=-1<0$, $\lambda_{2}=r-a>0$, hence $E_{1}=(1,0)$ is a hyperbolic saddle.
The Jacobian matrix at $E_*=(1-\frac{\beta b}{1+mb^2},b-\frac{\beta b^2}{1+mb^2})$ is
\begin{eqnarray*}
J_{E_*}=&&\left(
\begin{array}{cc}
\frac{2\beta b}{(1+mb^2)^2}-1 & \frac{\beta (mb^2-1)}{(1+mb^2)^2} \\
rb^2 & -rb
\end{array}
\right)\\
\triangleq&&\left(
\begin{array}{cc}
r_0 & r_1\\
r_2 & -r_3
\end{array}
\right),
\end{eqnarray*}
where $r_0=\frac{2\beta b}{(1+mb^2)^2}-1$, $r_1=\frac{\beta (mb^2-1)}{(1+mb^2)^2}$, $r_3=rb$, $r_2=br_3$, then
\begin{eqnarray*}
\left\{ {\begin{array}{*{20}{l}}
T_0(\beta)=r_0-r_3=\frac{2\beta b-(1+mb^2)^2}{(1+mb^2)^2}-rb ,\\
D_0(\beta)=-r_3(r_0+br_1)=-rb(\frac{\beta b}{1+mb^2}-1),
\end{array}} \right.
\end{eqnarray*}
when $\beta <\frac{1+mb^2}{b}$, then $\frac{\beta b}{1+mb^2}-1<0$, it is easy to know that $D_0(\beta)>0$;
if $\beta <\frac{(1+mb^2)^2}{2b}$, which is $r_0<0$, we know $r_3>0$. It follows that $r_0-r_3<0$, so $T_0(\beta)<0$. It can conclude that $\beta<\min\{\frac{1+mb^2}{b}, \frac{(1+mb^2)^2}{2b}\}$, we can get that $T_0(\beta)<0,D_0(\beta)>0$, therefore $E_*$ is locally asymptotically stable.
If $00$ , and when $\beta <\beta_H$, it is obtained that $\frac{2\beta b-(1+mb^2)^2}{(1+mb^2)^2}-br<0$, so $T_0(\beta)<0,D_0(\beta)>0$, it follows that the eigenvalues are negative or have negative real parts. If $\beta >\beta_H$, then $T_0(\beta)>0$, thus Eq. (3) has positive eigenvalues or eigenvalues with positive real parts.
If $\beta =\beta_H$, the transverse condition of the Hopf bifurcation can be satisfied
\begin{equation*}
\frac{dRe\lambda(\beta)}{d\beta}\mid _{\beta=\beta_{H}}=\frac{b}{(1+mb^2)^2}>0.
\end{equation*}
Therefore, when $\beta =\beta_H$, the Hopf bifurcation may occur.
\end{proof}
When $r=a$, the corresponding eigenvalues are $\lambda_{1}=-1$ and $\lambda_{2}=0$, so the equilibrium $E_1$ is nonhyperbolic and by using center manifold reduction theory, we investigate the stability of the equilibrium $E_1$ in \cite{30}.
Let $u\rightarrow u+1$, $v\rightarrow v$, then
\begin{eqnarray}
\left(
\begin{array}{*{20}{l}}
\frac{du}{dt} \\[5pt]
\frac{dv}{dt}
\end{array}
\right)
=
\left(
\begin{array}{cc}
f(u+1,v) \\
g(u+1,v)
\end{array}
\right)
=&&J_{E_1}\left(
\begin{array}{cc}
u \\
v
\end{array}
\right)
+\left(
\begin{array}{cc}
f_2(u,v)\\
g_2(u,v)
\end{array}
\right),
\end{eqnarray}
where
\begin{equation*}
f_2(u,v)=-u^2+\beta mv^3, \qquad g_2(u,v)=-rv^2+ruv^2,
\end{equation*}
and $J_{E_1}$ can be found as before. Let
\begin{eqnarray*}
T=\left(
\begin{array}{cc}
\beta & 1 \\
-1 & 0
\end{array}
\right),
\qquad
T^{-1}=\left(
\begin{array}{cc}
0 & -1 \\
1 & \beta
\end{array}
\right),
\end{eqnarray*}
then,
\begin{eqnarray*}
T^{-1}J_{E_1}T
=\left(
\begin{array}{cc}
0 & 0 \\
0 & -1
\end{array}
\right).
\end{eqnarray*}
The transformation $\left(\begin{array}{cc}u \\ v\end{array}\right)=T\left(\begin{array}{cc}z_1 \\ z_2\end{array}\right)$, then the system
can be represented as
\begin{eqnarray}
\left(
\begin{array}{*{20}{l}}
\frac{dz_1}{dt} \\
\frac{dz_2}{dt}
\end{array}
\right)
=&&\left(
\begin{array}{cc}
0 & -1 \\
1 & \beta
\end{array}
\right)
\left(
\begin{array}{cc}
z_1 \\
z_2
\end{array}
\right)\nonumber\\
&&+\left(
\begin{array}{cc}
rz_1^2-r\beta z_1^3-rz_1^2z_2 \\
-(\beta^2+r\beta)z_1^2-2\beta z_1z_2-z_2^2-(\beta m-r\beta^2)z_1^3+\beta z_1^2z_2
\end{array}
\right).
\end{eqnarray}
There exists a local center manifold which can be stated as
\begin{equation*}
W^c(0)=\{(z_1,z_2)\in R^2|z_2=\gamma(z_1),|z_1|<\delta,\gamma(0)=D\gamma(0)=0\},
\end{equation*}
for $\delta>0$ sufficiently small. Assume that $\gamma(z_1)$ has the following form
\begin{equation}
\gamma(z_1)=ez_1^2+fz_1^3+\mathcal{O}(z_1^4).
\end{equation}
Substituting Eq. (9) into Eq. (8), equating coefficients on each power of $s$ to zero, then we have
\begin{equation*}
e=\frac{(r+\beta)\beta}{2},\quad f=\frac{(r\beta-m)\beta-(r+\beta)^2\beta}{2},
\end{equation*}
therefore,
\begin{equation*}
\gamma(z_1)=\frac{(r+\beta)\beta}{2}z_1^2+\frac{(r\beta-m)\beta-(r+\beta)^2\beta}{2}z_1^3+\mathcal{O}(z_1^4).
\end{equation*}
The system restricted to the local center manifold is given by
\begin{equation*}
\dot{z_1}=rz_1^2+\mathcal{O}(z_1^4),
\end{equation*}
thus the equilibrium $E_1$ is unstable.
\begin{remark}
Note that if $a=0$ (i.e. without harvesting effort), the model has a boundary equilibrium $E^*_1=(1,0)$, which is a hyperbolic saddle and always unstable. When $a\neq0$, it has three kinds of different cases, the boundary equilibrium $E_1=(1,0)$ is a unstable saddle if $0r$. Notice that if $a=0$ the system always has an internal equilibrium $E^*_2=(u^*_2,v^*_2)$, where $(u^*_2,v^*_2)=(1-\frac{\beta}{1+m},1-\frac{\beta}{1+m})$, $\beta<1+m$. One can get the following conclusions when $a\neq0$: the system has only one positive equilibrium when $00$. At the point $\beta =\beta_H$, the eigenvector corresponding to the eigenvalue $i\omega_0$ is $\xi=(i\omega_0+r_0,br_0)^{T}$. Let
\begin{eqnarray*}
P&=\left(
\begin{array}{cc}
\omega_0 & r_0 \\
0 & br_0
\end{array}
\right),
\end{eqnarray*}
and
\begin{eqnarray*}
P^{-1}&=\left(
\begin{array}{cc}
\frac{1}{\omega_0} & -\frac{1}{bw_0} \\
0 & \frac{1}{br_0}
\end{array}
\right).
\end{eqnarray*}
Take the transformation $\left(\begin{array}{cc}u \\ v\end{array}\right)=P\left(\begin{array}{cc}s \\ h\end{array}\right)$, the system can become
\begin{eqnarray*}
\left(
\begin{array}{*{20}{l}}
\frac{ds}{dt} \\[5pt]
\frac{dh}{dt}
\end{array}
\right)
=&&P^{-1}J_{E_*}P
\left(
\begin{array}{cc}
s \\
h
\end{array}
\right)
+
P^{-1}\left(
\begin{array}{cc}
f_2P(s,h) \\
g_2P(s,h)
\end{array}
\right)\\
=&&\left(
\begin{array}{cc}
0 & -\omega_0 \\
\omega_0 & 0
\end{array}
\right)
\left(
\begin{array}{cc}
s \\
h
\end{array}
\right)
+\left(
\begin{array}{cc}
f_3(s,h) \\
g_3(s,h)
\end{array}
\right),
\end{eqnarray*}
where
\begin{eqnarray*}
\left(
\begin{array}{cc}
f_3(s,h) \\
g_3(s,h)
\end{array}
\right)
=\left(
\begin{array}{cc}
\frac{1}{\omega_0}f_2(\omega_0s+r_0h,br_0h)-\frac{1}{b\omega_0}g_2(\omega_0s+r_0h,br_0h) \\[5pt]
\frac{1}{br_0}g_2(\omega_0s+r_0h,br_0h)
\end{array}
\right),
\end{eqnarray*}
with
\begin{eqnarray*}
f_2(\omega_0s+r_0h,br_0h)
&=&\omega_0^2f_{uu}s^2+r_0^2(f_{uu}+bf_{uv}+b^2f_{vv})h^2+\omega_0r_0(2f_{uu}+bf_{uv})sh+\omega_0^3f_{uuu}s^3\\
&&+\omega_0^2r_0(3f_{uuu}+bf_{uuv})s^2h+\omega_0r_0^2(3f_{uuu}+2bf_{uuv}+b^2f_{uvv})sh^2\\
&&+r_0^3(f_{uuu}+bf_{uuv}+b^2f_{uvv}+b^3f_{vvv})h^3,
\\
g_2(\omega_0s+r_0h,br_0h)
&=& \omega_0^2g_{uu}s^2+r_0^2(g_{uu}+bg_{uv}+b^2g_{vv})h^2+\omega_0r_0(2g_{uu}+bg_{uv})sh+\omega_0^3g_{uuu}s^3\\
&&+\omega_0^2r_0(3g_{uuu}+bg_{uuv})s^2h+\omega_0r_0^2(3g_{uuu}+2bg_{uuv}+b^2g_{uvv})sh^2\\
&&+r_0^3(g_{uuu}+bg_{uuv}+b^2g_{uvv}+b^3g_{vvv})h^3,
\end{eqnarray*}
hence, we get
\begin{eqnarray*}
&&f_{2ss}=\omega_0^2f_{uu}, \quad f_{2sh}=\omega_0r_0(2f_{uu}+bf_{uv}), \quad f_{2hh}=r_0^2(f_{uu}+bf_{uv}+b^2f_{vv}),\\
&&f_{2sss}=\omega_0^3f_{uuu}, \quad f_{2hhh}=r_0^3(f_{uuu}+bf_{uuv}+b^2f_{uvv}+b^3f_{vvv}),\\ &&f_{2ssh}=\omega_0^2r_0(3f_{uuu}+bf_{uuv}), \quad f_{2shh}=\omega_0r_0^2(3f_{uuu}+2bf_{uuv}+b^2f_{uvv}), \\
&&g_{2ss}=\omega_0^2g_{uu}, \quad g_{2sh}=\omega_0r_0(2g_{uu}+bg_{uv}), \quad g_{2hh}=r_0^2(g_{uu}+bg_{uv}+b^2g_{vv}),\\
&&g_{2sss}=\omega_0^3g_{uuu}, \quad g_{2hhh}=r_0^3(g_{uuu}+bg_{uuv}+b^2g_{uvv}+b^3g_{vvv}),\\
&&g_{2ssh}=\omega_0^2r_0(3g_{uuu}+bg_{uuv}), \quad g_{2shh}=\omega_0r_0^2(3g_{uuu}+2bg_{uuv}+b^2g_{uvv}).
\end{eqnarray*}
By \cite{30}, stability of the Hopf bifurcation is determined by the sign of
\begin{equation*}
\sigma=\frac{1}{16}(f_{3sss}+g_{3ssh}+f_{3shh}+g_{3hhh})\\
+\frac{1}{16\omega_0}(f_{3sh}(f_{3ss}+f_{3hh})-g_{3sh}(g_{3ss}+g_{3hh})-f_{3ss}g_{3ss}+f_{3ss}g_{3hh}),
\end{equation*}
where, all partial derivatives are evaluated at $(s,h,\beta)=(0,0,\beta_H)$,
\begin{eqnarray*}
&&f_{3ss}=\frac{\omega_0}{b}(bf_{uu}-g_{uu}), \quad f_{3sss}=\frac{\omega_0^2}{b}(bf_{uuu}-g_{uuu}), \quad f_{3sh}=\frac{r_0}{b}(2bf_{uu}+b^2f_{uv}-2g_{uu}-bg_{uv}), \\ &&f_{3ssh}=\frac{\omega_0r_0}{b}(3bf_{uuu}+b^2f_{uuv}-3g_{uuu}-bg_{uuv}),\quad f_{3hh}=\frac{r_0^2}{b\omega_0}(bf_{uu}+b^2f_{uv}+b^3f_{vv}-g_{uu}-bg_{uv}-b^2g_{vv}), \\
&&f_{3shh}=\frac{r_0^2}{b}(3bf_{uuu}+2b^2f_{uuv}+b^3f_{uvv}-3g_{uuu}-2bg_{uuv}-b^2g_{uvv}), \\
&&f_{3hhh}=\frac{r_0^3}{b\omega_0}(bf_{uuu}+b^2f_{uuv}+b^3f_{uvv}+b^4f_{vvv}-g_{uuu}-bg_{uuv}-b^2g_{uvv}),\\
&&g_{3ss}=\frac{\omega_0^2}{br_0}g_{uu}, \quad g_{3sh}=\frac{\omega_0}{b}(2g_{uu}+bg_{uv}), \quad g_{3hh}=\frac{r_0}{b}(g_{uu}+bg_{uv}+b^2g_{vv}),\\
&&g_{3sss}=\frac{\omega_0^3}{br_0}g_{uuu}, \quad g_{3hhh}=\frac{r_0^2}{b}(g_{uuu}+bg_{uuv}+b^2g_{uvv}),\\ &&g_{3shh}=\frac{\omega_0r_0}{b}(3g_{uuu}+2bg_{uuv}+b^2g_{uvv}), \quad g_{3ssh}=\frac{\omega_0^2}{b}(3g_{uuu}+bg_{uuv}).
\end{eqnarray*}
\begin{figure}[t!]
\centering
\centering
\begin{tabular}{ccc}
\subfigure[]{
\includegraphics[width=2.2in,height=1.5in]{Fig110.eps}}
\subfigure[]{
\includegraphics[width=2.2in,height=1.5in]{Fig111.eps}}
\subfigure[]{
\includegraphics[width=2.2in,height=1.5in]{Fig112.eps}}\\
\end{tabular}
\caption{The supercritical Hopf bifurcation occurs at $\beta=\beta_H$. ($a=0.4, m=0.8, r=0.89)$.}
\label{fig1}
\end{figure}
Now, we claim our main result as below.
\begin{theorem}
If $\sigma<0(resp. >0)$, then the direction of the Hopf bifurcation is supercritical (resp. subcritical) and the periodic solution bifurcated from the Hopf bifurcation is stable (resp. unstable).
\end{theorem}
\begin{remark}
In order to verify the above theoretical analysis in Theorem 2.4, setting $a=0.4, m=0.8, r=0.89$, and we take $\beta=\beta_H=2.089$, then we can get the positive equilibrium $E_*=(0.0743, 0.0409)$ is locally asymptotically stable with the evolution of time $t$, and $\sigma=-0.0788<0$, a supercritical Hopf bifurcation occurs and the periodic solution bifurcated from the Hopf bifurcation is stable.
\end{remark}
\subsection{Transcritical bifurcation}
\noindent
According to the results of Remark 2.1, when $a=r$, the system has no positive equilibrium and a boundary equilibrium $E_1$, $E_1$ is nonhyperbolic, then the determinant of corresponding Jacobian matrix at $E_1$ is equal to zero if $a=r$, there is chance for the existence of transcritical bifurcation. By using the Sotomayer’s theorem, we can verify the occurrence of transcritical bifurcation.
\newpage
\begin{theorem}
System (3) undergoes the transcritical bifurcation when $a=a_{TC}=r$.
\end{theorem}
\begin{proof}
Let $\bar{x}$ and $\bar{y}$ be eigenvectors corresponding to zero eigenvalues of matrices $J_{E_1}$ and $J^T_{E_1}$, respectively, where $J^T_{E_1}$ is the transposed matrix of $J_{E_1}$. They can be written as follows
\begin{eqnarray*}
\bar{x}=\left(
\begin{array}{c}
x_1 \\
x_2
\end{array}
\right)
=\left(
\begin{array}{c}
-\beta \\
1
\end{array}
\right),
\quad
\bar{y}=\left(
\begin{array}{c}
y_1 \\
y_2
\end{array}
\right)
=\left(
\begin{array}{c}
0 \\
1
\end{array}
\right).
\end{eqnarray*}
Additionally, we can get
\begin{eqnarray*}
F_a(E_1;a_{TC})=\left(
\begin{array}{c}
0 \\
-v
\end{array}
\right)_{(E_1;a_{TC})}
=\left(
\begin{array}{c}
0 \\
0
\end{array}
\right),
\end{eqnarray*}
\begin{eqnarray*}
DF_a(E_1;a_{TC})\bar{x}=\left(
\begin{array}{cc}
0&0 \\
0&-1
\end{array}
\right)
\left(
\begin{array}{c}
-\beta \\
1
\end{array}
\right)_{(E_1;a_{TC})}
=\left(
\begin{array}{c}
0 \\
-1
\end{array}
\right),
\end{eqnarray*}
\begin{eqnarray*}
D^2F_a(E_1;a_{TC})(\bar{x},\bar{x})=\left(
\begin{array}{c}
\frac{\partial^2F_1}{\partial u^2}x_1^2+2\frac{\partial^2F_1}{\partial u\partial v}x_1x_2+\frac{\partial^2F_1}{\partial v^2}x_2^2 \\[2mm]
\frac{\partial^2F_2}{\partial u^2}x_1^2+2\frac{\partial^2F_2}{\partial u\partial v}x_1x_2+\frac{\partial^2F_2}{\partial v^2}x_2^2
\end{array}
\right)_{(E_1;a_{TC})}
=\left(
\begin{array}{c}
0 \\
-2r
\end{array}
\right),
\end{eqnarray*}
so we get
\begin{eqnarray*}
&&\bar{y}^TF_a(E_1;a_{TC})=0,\\
&&\bar{y}^T[DF_a(E_1;a_{TC})\bar{x}]=-1\neq0,\\
&&\bar{y}^T[D^2F(E_1;a_{TC})(\bar{x},\bar{x})]=-2r\neq0.
\end{eqnarray*}
Thus, by the Sotomayor’s Theorem, then system (3) undergoes the transcritical bifurcation.
\end{proof}
\section{Turing instability }
\noindent
In order to find out the critical values of the Turing instability and get the Turing instability region, the linearized form of
Eq. (2) at $E_*$ is given as
\begin{equation}
\dot{U}=J_{E_*}U+D\Delta U,
\end{equation}
where
\begin{eqnarray*}
U=\left(
\begin{array}{c}
u-u_* \\
v-v_*
\end{array}
\right),
\quad
D=\left(
\begin{array}{cc}
d_1&0 \\
0&d_2
\end{array}
\right),
\end{eqnarray*}
and $J_{E_*}$ can be found as before.
The Fourier expansion of the solution to Eq. (11) can be described as
\begin{eqnarray}
U=\sum^{\infty}_{k=0}\left(
\begin{array}{c}
a_k \\
b_k
\end{array}
\right)e^{\lambda t+i\mathbf{k}\cdot\mathbf{\gamma}},
\end{eqnarray}
$\lambda $ is the growth rate of perturbation, and the wave number $\mathbf{k}=(k_x,k_y)$ satisfies $k=\mid \mathbf{k}\mid$, $\mathbf{\gamma}=(x,y)$ denotes the spatial vector in two-dimensional space, $i$ is the imaginary unit and $i^2=-1$, $a_k$ and $b_k$ represent constants.
we can plug Eq. (12) into Eq. (11), then
\begin{eqnarray*}
\lambda\left(
\begin{array}{c}
a_k \\
b_k
\end{array}
\right)
=\left(
\begin{array}{cc}
r_0-k^2d_1 &r_1 \\
br_3 &-r_3-k^2d_2
\end{array}
\right)
\left(
\begin{array}{c}
a_k \\
b_k
\end{array}
\right),
\end{eqnarray*}
\begin{figure}[t!]
\centering
\includegraphics[width=3.5in,height=2.5in]{Fig2111.eps}
\caption{Bifurcation diagram. $D_1$: Turing-Hopf instability region; $D_2$: Turing instability region; $D_3$: stable region; $D_4$: Turing instability region. ($m=0.25$, $r=1$, $d_1=0.1$, $d_2=2.5)$.}
\label{fig2}
\end{figure}
\noindent
then the characteristic equation can be represented as
\begin{equation}
\ell_k(\beta)=\lambda^2-T_k(\beta)\lambda+D_k(\beta)=0,
\end{equation}
where
\begin{eqnarray}
\left\{ {\begin{array}{*{20}{l}}
T_k(\beta)=-(d_1+d_2)k^2+r_0-r_3,\\
D_k(\beta)=d_1d_2k^4+(r_3d_1-r_0d_2)k^2-r_3(r_0+br_1),
\end{array}} \right.
\end{eqnarray}
then
\begin{equation}
\lambda_k=\frac{T_k(\beta)\pm \sqrt{(T_k(\beta))^2-4D_k(\beta)}}{2}.
\end{equation}
The Turing instability means that Eq. (3) is stable, while Eq. (2) is unstable at $E_*$, that is to say, under the conditions $T_0(\beta)<0$, $D_0(\beta)>0$, $T_k(\beta)<0$ and $D_k(\beta)<0$, for any $k\neq0$, the Turing instability can be obtained. Thus we obtain the condition for the marginal stability as follows
\begin{equation}
\min(D_{k_c^2}(\beta))=0,
\end{equation}
where $k_c^2=\sqrt{\frac{D_0(\beta)}{d_1d_2}}=\sqrt{\frac{-r_3(r_0+br_1)}{d_1d_2}}$.
On the basis of the Eq. (16), we can obtain the critical values of Turing instability,
\begin{equation}
\beta_T=\frac{(1+mb^2)^2}{2bd_2}(\sqrt{r_3d_1(1-mb^2)(2d_2-r_3d_1(1+mb^2))}+d_2-mb^2r_3d_1).
\end{equation}
\begin{figure}[t!]
\centering
\subfigure[]{
\includegraphics[width=3.0in]{Fig211.eps}}
\subfigure[]{
\includegraphics[width=3.4in]{Fig212.eps}}
\caption{(a) plot of $D(k^2)$ for different values of parameter $\beta$, (b) plot of $Re(\lambda_{k^2})$ for different values of parameter $\beta$. ($a=0.2$, $m=0.25$, $r=1$, $d_1=0.1$, $d_2=2.5$). }
\label{fig3}
\end{figure}
\begin{remark}
According to the above analysis, we know that there are two critical values of the Turing instability
\begin{eqnarray*}
&&\beta_+=\frac{(1+mb^2)^2}{2bd_2}(\sqrt{r_3d_1(1-mb^2)(2d_2-r_3d_1(1+mb^2))}+d_2-mb^2r_3d_1),\\ &&\beta_-=\frac{(1+mb^2)^2}{2bd_2}(-\sqrt{r_3d_1(1-mb^2)(2d_2-r_3d_1(1+mb^2))}+d_2-mb^2r_3d_1).
\end{eqnarray*}
It can be deduced that $E_*(u_*,v_*)$ is stable when $\beta\in(\beta_-,\beta_+)$, and it is unstable when $\beta\in(-\infty,\beta_-)\cup(\beta_+,+\infty)$. In this paper, we take the critical values of the Turing instability as $\beta_T=\beta_+$.
In Fig.2, the bifurcation diagram is plotted with $m=0.25$, $r=1$, $d_1=0.1$, $d_2=2.5$, so it is shown that the domain is divided into four regions, from top to bottom: Turing-Hopf instability region, Turing instability region, stable region and Turing instability region. The Turing instability region is the space between the Turing bifurcation curve and the Hopf bifurcation curve, so Turing instability occurs when $\beta \in (\beta_T, \beta_H).$
In Fig.3, take the parameters $a=0.2$, $m=0.25$, $r=1$, $d_1=0.1$, $d_2=2.5$, so we get $\beta_T=1.0299$, $\beta_H=1.5138$. It is derived that $D_k(\beta)=0$ when $\beta=1.0299$, then there is no Turing instability. When $D_k(\beta)<0$, $Re(\lambda_{k^2})>0$, the Turing instability occurs,
as the value of the parameter $\beta$ increases.
\end{remark}
\section{Weakly nonlinear analysis }
\subsection{Amplitude equation }
\noindent
In this subsection, the weakly nonlinear analysis is applied to derive the amplitude equation near the critical value of the Turing instability.
Consider the following form
\begin{eqnarray}
\frac{\partial U }{\partial t}=LU+F(U,U),
\end{eqnarray}
suppose $U=(u,v)^{T}$,
where
\begin{eqnarray*}
L=\left(
\begin{array}{cc}
r_0+d_1\Delta &r_1 \\
br_3 &-r_3+d_2\Delta
\end{array}
\right),
\end{eqnarray*}
with
\begin{eqnarray*}
F=\left(
\begin{array}{c}
f_{uu}u^2+f_{uv}uv+f_{vv}v^2+f_{uuu}u^3+f_{uuv}u^2v+f_{uvv}uv^2+f_{vvv}v^3\\
g_{uu}u^2+g_{uv}uv+g_{vv}v^2+g_{uuu}u^3+g_{uuv}u^2v+g_{uvv}uv^2+g_{vvv}v^3
\end{array}
\right).
\end{eqnarray*}
Based on the method of small parameter expansion, we introduce the parameter $\varepsilon$ to analyze the dynamical behaviors of the system when it is near the critical point $\beta=\beta_T$,
\begin{equation}
\beta-\beta_T=\varepsilon \beta_1+\varepsilon^2 \beta_2+\varepsilon^3 \beta_3+\mathcal{O}(4),
\end{equation}
further introduce time scales $T_0=t, T_1=\varepsilon t, T_2=\varepsilon^2 t, T_3=\varepsilon^3 t,$ then
\begin{equation}
\frac{\partial}{\partial t}=\frac{\partial}{\partial T_0}+\varepsilon \frac{\partial}{\partial T_1}+\varepsilon^2 \frac{\partial}{\partial T_2}+\varepsilon^3 \frac{\partial}{\partial T_3}+\mathcal{O}(4),
\end{equation}
the variable $U$ can be expanded in the following form with the small parameter $\varepsilon $
\begin{eqnarray}
U=\left(
\begin{array}{c}
u \\
v
\end{array}
\right)
=\varepsilon \left(
\begin{array}{c}
u_1 \\
v_1
\end{array}
\right)
+\varepsilon^2 \left(
\begin{array}{c}
u_2 \\
v_2
\end{array}
\right)
+\varepsilon^3 \left(
\begin{array}{c}
u_3 \\
v_3
\end{array}
\right)+\mathcal{O}(4),
\end{eqnarray}
the expansion form of the nonlinear term $F$ can be expressed as
\begin{equation}
F=\varepsilon^2 F_2+\varepsilon^3 F_3+\mathcal{O}(4),
\end{equation}
where
\begin{eqnarray*}
F_2=\left(
\begin{array}{c}
f_{uu}u_1^2+f_{uv}u_1v_1+f_{vv}v_1^2 \\
g_{uu}u_1^2+g_{uv}u_1v_1+g_{vv}v_1^2
\end{array}
\right),
\end{eqnarray*}
and
\begin{eqnarray*}
F_3=\left(
\begin{array}{c}
f_{uu}u_1u_2+f_{uv}(u_1v_2+u_2v_1)+f_{vv}v_1v_2+f_{uuu}u_1^3+f_{uuv}u_1^2v_1+f_{uvv}u_1v_1^2+f_{vvv}v_1^3 \\
g_{uu}u_1u_2+g_{uv}(u_1v_2+u_2v_1)+g_{vv}v_1v_2+g_{uuu}u_1^3+g_{uuv}u_1^2v_1+g_{uvv}u_1v_1^2+g_{vvv}v_1^3
\end{array}
\right).
\end{eqnarray*}
The linear operator L can be decomposed into
\begin{equation}
L=L_c+(\beta -\beta_T)M,
\end{equation}
with
\begin{eqnarray*}
L_c=\left(
\begin{array}{cc}
r_0(\beta_T)+d_1\Delta &r_1(\beta_T) \\
br_3(\beta_T) &-r_3(\beta_T)+d_2\Delta
\end{array}
\right),\quad
M=\left(
\begin{array}{cc}
\frac{2b}{(1+mb^2)^2} &\frac{(mb^2-1)}{(1+mb^2)^2} \\
0 &0
\end{array}
\right).
\end{eqnarray*}
Substituting formulas of decomposition Eq. (19)-(23) into Eq. (18), and according to different orders of $\varepsilon$, we have
$\mathcal{O}(\varepsilon)$:
\begin{eqnarray}
L_c\left(
\begin{array}{c}
u_1 \\
v_1
\end{array}
\right)=0.
\end{eqnarray}
$\mathcal{O}(\varepsilon^2)$:
\begin{eqnarray}
L_c\left(
\begin{array}{c}
u_2 \\
v_2
\end{array}
\right)
=\frac{\partial}{\partial T_1}
\left(
\begin{array}{c}
u_1 \\
v_1
\end{array}
\right)
-\beta_1M
\left(
\begin{array}{c}
u_1 \\
v_1
\end{array}
\right)-F_2.
\end{eqnarray}
$\mathcal{O}(\varepsilon^3)$:
\begin{eqnarray}
L_c\left(
\begin{array}{c}
u_3 \\
v_3
\end{array}
\right)
=&\dfrac{\partial}{\partial T_1}
\left(
\begin{array}{c}
u_2 \\
v_2
\end{array}
\right)
+\dfrac{\partial}{\partial T_2}
\left(
\begin{array}{c}
u_1 \\
v_1
\end{array}
\right)
-\beta_1M
\left(
\begin{array}{c}
u_2 \\
v_2
\end{array}
\right)
-\beta_2M
\left(
\begin{array}{c}
u_1 \\
v_1
\end{array}
\right)-F_3.
\end{eqnarray}
For Eq. (24), $L_c$ is the linear operator of the system at critical point $\beta=\beta_T$, then the solution of Eq. (24) is the linear
combination of the eigenvectors corresponding to the zero eigenvalue. The form of the solution is
\begin{eqnarray}
\left(
\begin{array}{c}
u_1 \\
v_1
\end{array}
\right)\mathbf{}
=\left(
\begin{array}{c}
\phi \\
1
\end{array}
\right)(\sum^3_{j=1}W_j\exp(i\mathbf{k_j}\cdot \gamma)+\sum^3_{j=1}\bar{W}_j\exp(-i\mathbf{k_j}\cdot \gamma)),
\end{eqnarray}
where $\phi=\frac{r_1}{d_1k_c^2-r_0}$, and $|\mathbf{k_j}|=k_c$, $W_j$ is the amplitude of the mode $\exp(i\mathbf{k_j}\cdot \gamma)$ under
the first perturbation, which is determined by high-order terms of the perturbation.
In the light of Fredholm solvability condition, we get the non-trivial solution of the Eq. (25), that is, the right-hand-side vector function of Eq. (25) must be orthogonal to the corresponding zero eigenvectors of operator $L_c^+$, which is the adjoint operator of linear operator $L_c$ here, so the zero eigenvector of operator $L_c^+$ can be written as
\begin{eqnarray}
\left(
\begin{array}{c}
1 \\
\varphi
\end{array}
\right)(\exp(-i\mathbf{k_j}\cdot \gamma))+c.c., \quad j=1,2,3,
\end{eqnarray}
where $\varphi=\frac{d_1k_c^2-r_0}{br_3}$, c.c. denotes the conjugate of the former terms. For Eq. (25), let
\begin{eqnarray}
L_c\left(
\begin{array}{c}
u_2 \\
v_2
\end{array}
\right)
=&&\frac{\partial}{\partial T_1}
\left(
\begin{array}{c}
u_1 \\
v_1
\end{array}
\right)
-\beta_1M
\left(
\begin{array}{c}
u_1 \\
v_1
\end{array}
\right)-F_2\nonumber \\
\triangleq && \left(
\begin{array}{c}
H_u \\
H_v
\end{array}
\right),
\end{eqnarray}
then the orthogonality condition is equivalent to
\begin{eqnarray}
(
1 , \varphi
)\exp(-i\mathbf{k_j}\cdot \gamma)
\left(
\begin{array}{c}
H_u^j \\
H_v^j \\
\end{array}
\right)=0, \quad j=1,2,3,
\end{eqnarray}
where $H_u^j$ and $H_v^j$ represent the corresponding coefficients of $\exp(i\mathbf{k_j}\cdot \gamma)$ in $H_u$ and $H_v $, respectively. The Fredholm solvability condition Eq. (30) is applied to Eq. (25). It is derived that
\begin{eqnarray}
&&(\phi+\varphi)\frac{\partial W_1}{\partial T_1}=\beta_1m_1W_1+2(h_1+\varphi h_2)\bar{W_2}\bar{W_3},\nonumber \\
&&(\phi+\varphi)\frac{\partial W_2}{\partial T_1}=\beta_1m_1W_2+2(h_1+\varphi h_2)\bar{W_1}\bar{W_3},\nonumber \\
&&(\phi+\varphi)\frac{\partial W_3}{\partial T_1}=\beta_1m_1W_3+2(h_1+\varphi h_2)\bar{W_1}\bar{W_2},
\end{eqnarray}
where
\begin{eqnarray*}
&&m_1=\frac{2\phi b+(mb^2-1)}{(1+mb^2)^2},\\
&&h_1=f_{uu}\phi^2+f_{uv}\phi+f_{vv},\\
&&h_2=g_{uu}\phi^2+g_{uv}\phi+g_{vv}.
\end{eqnarray*}
Eq. (31) is the amplitude equation of the first order perturbation. Substitute Eq. (27) into Eq. (25), then we can obtain
\begin{eqnarray}
\left(
\begin{array}{c}
u_2 \\
v_2
\end{array}
\right)
=&&\left(
\begin{array}{c}
U_0 \\
V_0
\end{array}
\right)
+\sum^3_{j=1}
\left(
\begin{array}{c}
U_j \\
V_j
\end{array}
\right)e^{i\mathbf{\mathbf{k_j}}\cdot \gamma}
+\sum^3_{j=1}
\left(
\begin{array}{c}
U_{jj} \\
V_{jj} \\
\end{array}
\right)e^{2i\mathbf{\mathbf{k_j}}\cdot \gamma}\nonumber\\
&&+\left(
\begin{array}{c}
U_{12} \\
V_{12}
\end{array}
\right)e^{i(\mathbf{k_1}-\mathbf{k_2})\cdot \gamma}
+\left(
\begin{array}{c}
U_{23} \\
V_{23}
\end{array}
\right)e^{i(\mathbf{k_2}-\mathbf{k_3})\cdot \gamma}
+\left(
\begin{array}{c}
U_{31} \\
V_{31}
\end{array}
\right)e^{i(\mathbf{k_3}-\mathbf{k_1})\cdot \gamma}+c.c.,
\end{eqnarray}
where
\begin{eqnarray*}
&&\left(
\begin{array}{c}
U_0 \\
V_0
\end{array}
\right)
=\left(
\begin{array}{c}
u_{00} \\
v_{00}
\end{array}
\right)(|W_1|^2+|W_2|^2+|W_3|^2),
\quad U_j=\phi V_j,\\
&&\left(
\begin{array}{c}
U_{jj} \\
U_{jj}
\end{array}
\right)
=\left(
\begin{array}{c}
u_{11} \\
v_{11}
\end{array}
\right)W_j^2,
\quad \left(
\begin{array}{c}
U_{jl} \\
U_{jl}
\end{array}
\right)
=\left(
\begin{array}{c}
u_c \\
v_c
\end{array}
\right)W_j\bar{W_l},
\end{eqnarray*}
and
\begin{equation*}
\left(
\begin{array}{c}
u_{00} \\
v_{00}
\end{array}
\right)
=\left(
\begin{array}{c}
\frac{2(r_1h_2+r_3h_1)}{-r_3(r_0+br_1)} \\[5pt]
\frac{2(br_3h_1-r_0h_2)}{-r_3(r_0+br_1)}
\end{array}
\right),
\end{equation*}
\begin{equation*}
\left(
\begin{array}{c}
u_{11} \\
v_{11}
\end{array}
\right)
=\left(
\begin{array}{c}
\frac{r_1h_2+(r_3+4d_2k_c^2)h_1}{(r_0-4d_1k_c^2)(-r_3-4d_2k_c^2)-br_1r_3} \\[5pt]
\frac{br_3h_1-(r_0-4d_1k_c^2)h_2}{(r_0-4d_1k_c^2)(-r_3-4d_2k_c^2)-br_1r_3}
\end{array}
\right),
\end{equation*}
\begin{equation*}
\left(
\begin{array}{c}
u_c \\
v_c
\end{array}
\right)
=\left(
\begin{array}{c}
\frac{2(r_1h_2+(r_3+3d_2k_c^2)h_1)}{(r_0-3d_1k_c^2)(-r_3-3d_2k_c^2)-br_1r_3} \\[5pt]
\frac{2(br_3h_1-(r_0-3d_1k_c^2)h_2)}{(r_0-3d_1k_c^2)(-r_3-3d_2k_c^2)-br_1r_3}
\end{array}
\right).
\end{equation*}
For convenience, $(R_u, R_v)^\mathrm{T}$ denotes the right side of Eq. (26). $R_u^j$ and $R_v^j$ stand for the corresponding coefficients of $\exp(i\mathbf{k_j}\cdot \gamma)$ in $R_u$ and $R_v $, respectively, where $j=1,2,3$. One has
\begin{eqnarray*}
R_u^1=&&\phi \frac{\partial V_1}{\partial T_1}+\phi \frac{\partial W_1}{\partial T_2}-\beta_1m_1V_1-\beta_2m_1W_1-2h_1(\bar{W_2}\bar{V_3}+\bar{W_3}\bar{V_2})\\
&&-(I_1|W_1|^2+I_2(|W_2|^2+|W_3|^2))W_1,\\
R_v^1=&&\frac{\partial V_1}{\partial T_1}+\frac{\partial W_1}{\partial T_2}-2h_2(\bar{W_2}\bar{V_3}+\bar{W_3}\bar{V_2})-(I_3|W_1|^2+I_4(|W_2|^2|W_3|^2))W_1,\\
R_u^2=&&\phi \frac{\partial V_2}{\partial T_1}+\phi \frac{\partial W_2}{\partial T_2}-\beta_1m_1V_2-\beta_2m_1W_2-2h_1(\bar{W_1}\bar{V_3}+\bar{W_3}\bar{V_1})\\
&&-(I_1|W_2|^2+I_2(|W_1|^2+|W_3|^2))W_2,\\
R_v^2=&&\frac{\partial V_2}{\partial T_1}+\frac{\partial W_2}{\partial T_2}-2h_2(\bar{W_1}\bar{V_3}+\bar{W_3}\bar{V_1})-(I_3|W_2|^2+I_4(|W_1|^2+|W_3|^2))W_2,\\
R_u^3=&&\phi \frac{\partial V_3}{\partial T_1}+\phi \frac{\partial W_3}{\partial T_2}-\beta_1m_1V_3-\beta_2m_1W_3-2h_1(\bar{W_1}\bar{V_2}+\bar{W_2}\bar{V_1})\\
&&-(I_1|W_3|^2+I_2(|W_1|^2+|W_2|^2))W_3,\\
R_v^3=&&\frac{\partial V_3}{\partial T_1}+\frac{\partial W_3}{\partial T_2}-2h_2(\bar{W_1}\bar{V_2}+\bar{W_2}\bar{V_1})-(I_3|W_3|^2+I_4(|W_1|^2+|W_2|^2))W_3,
\end{eqnarray*}
and
\begin{eqnarray*}
I_1=&&(2\phi f_{uu}+f_{uv})(u_{00}+u_{11})+(\phi f_{uv}+2f_{vv})(v_{00}+v_{11})\\
&&+3(\phi^3 f_{uuu}+\phi^2 f_{uuv}+\phi f_{uvv}+f_{vvv}),\\
I_2=&&(2\phi f_{uu}+f_{uv})(u_{00}+u_c)+(\phi f_{uv}+2f_{vv})(v_{00}+v_c)\\
&&+6(\phi^3 f_{uuu}+\phi^2 f_{uuv}+\phi f_{uvv}+f_{vvv}),\\
I_3=&&(2\phi g_{uu}+g_{uv})(u_{00}+u_{11})+(\phi g_{uv}+2g_{vv})(v_{00}+v_{11})\\
&&+3(\phi^3 g_{uuu}+\phi^2 g_{uuv}+\phi g_{uvv}+g_{vvv}),\\
I_4=&&(2\phi g_{uu}+g_{uv})(u_{00}+u_c)+(\phi g_{uv}+2g_{vv})(v_{00}+v_c)\\
&&+6(\phi^3 g_{uuu}+\phi^2 g_{uuv}+\phi g_{uvv}+g_{vvv}).
\end{eqnarray*}
$(R_u, R_v)^\mathrm{T}$ must be orthogonal with eigenvectors corresponding to zero eigenvalue of operator $L_c^+$. Then the Fredholm solvability condition implies that
\begin{eqnarray}
\left(
\begin{array}{cc}
1 , \varphi \\
\end{array}
\right)\exp(-i\mathbf{k_j}\cdot \gamma)
\left(
\begin{array}{c}
R_u^j \\
R_v^j \\
\end{array}
\right)=0, \quad j=1,2,3,
\end{eqnarray}
then the Fredholm solvability condition Eq. (33) is applied, it follows that
\begin{eqnarray}
(\phi+\varphi)(\frac{\partial V_1}{\partial T_1}+\frac{\partial W_1}{\partial T_2})&=&m_1(\beta_1V_1+\beta_2W_1)+2(h_1+\varphi h_2)(\bar{W_2}\bar{V_3}+
\bar{W_3}\bar{V_2})\nonumber\\
&&+(G_1|W_1|^2+G_2(|W_2|^2+|W_3|^2))W_1,\nonumber \\
(\phi+\varphi)(\frac{\partial V_2}{\partial T_1}+\frac{\partial W_2}{\partial T_2})&=&m_1(\beta_1V_2+\beta_2W_2)+2(h_1+\varphi h_2)(\bar{W_1}\bar{V_3}+\bar{W_3}\bar{V_1})\nonumber\\
&&+(G_1|W_2|^2+G_2(|W_1|^2+|W_3|^2))W_2,\nonumber \\
(\phi+\varphi)(\frac{\partial V_3}{\partial T_1}+\frac{\partial W_3}{\partial T_2})&=&m_1(\beta_1V_3+\beta_2W_3)+2(h_1+\varphi h_2)(\bar{W_2}\bar{V_1}+\bar{W_1}\bar{V_2})\nonumber\\
&&+(G_1|W_3|^2+G_2(|W_1|^2+|W_2|^2))W_3,
\end{eqnarray}
where
\begin{equation*}
G_1=I_1+\varphi I_3,\quad G_2=I_2+\varphi I_4.
\end{equation*}
To investigate the Turing patterns, the form of solutions for Eq. (2) can be written as
\begin{eqnarray}
\left(
\begin{array}{c}
u \\
v
\end{array}
\right)
=(\sum^3_{j=1}A_j\exp(i\mathbf{k_j}\cdot \gamma)+\sum^3_{j=1}\bar{A}_j\exp(-i\mathbf{k_j}\cdot \gamma)),
\end{eqnarray}
where $|\mathbf{k_j}|=k_c$, $A_j=(A_j^u, A_j^v)^T$ and $\bar{A_j}=(\bar{A_j^u}$, $\bar{A_j^v})^T$ are amplitudes for the modes of $i\mathbf{k_j}$ and $-i\mathbf{k_j}$, respectively.\\
According to Eq. (20), it follows that
\begin{equation}
\frac{\partial A_j}{\partial t}=\varepsilon \frac{\partial A_j}{\partial T_1}+\varepsilon^2 \frac{\partial A_j}{\partial T_2}+\varepsilon^3 \frac{\partial A_j}{\partial T_3}+\mathcal{O}(4),
\end{equation}
suppose
\begin{equation}
A_j=\varepsilon W_j+\varepsilon^2 V_j+\mathcal{O}(3).
\end{equation}
Using Eq. (31), Eq. (34) and Eq. (36), then amplitude equations corresponding to $A_1$, $A_2$ and $A_3$ are given by
\begin{eqnarray}
&&\tau_0\frac{\partial A_1}{\partial t}=\mu A_1+h_0\bar{A_2}\bar{A_3}-(g_1|A_1|^2+g_2(|A_2|^2+|A_3|^2))A_1,\nonumber\\
&&\tau_0\frac{\partial A_2}{\partial t}=\mu A_2+h_0\bar{A_1}\bar{A_3}-(g_1|A_2|^2+g_2(|A_1|^2+|A_3|^2))A_2,\nonumber\\
&&\tau_0\frac{\partial A_3}{\partial t}=\mu A_3+h_0\bar{A_1}\bar{A_2}-(g_1|A_3|^2+g_2(|A_1|^2+|A_2|^2))A_3,
\end{eqnarray}
where
\begin{eqnarray*}
&&\tau_0=\frac{\phi+\varphi}{m_1\beta_T},\quad \mu=\frac{\beta-\beta_T}{\beta_T},\quad h_0=\frac{2(h_1+\varphi h_2)}{m_1\beta_T},\\
&&g_1=-\frac{G_1}{m_1\beta_T},\quad g_2=-\frac{G_2}{m_1\beta_T}.
\end{eqnarray*}
Amplitude equations in Eq. (38) can be decomposed as
\begin{equation}
A_j=\rho_j \exp(i\theta_j),\quad j=1,2,3,
\end{equation}
where $\rho_j=|A_j|$ and $\theta_j$ are represented as mode and phase angle, respectively. Substituting Eq. (39) into Eq. (38), and separating the real and imaginary part of the Eq. (38), we get
\begin{eqnarray}
&&\tau_0\frac{\partial \theta}{\partial t}=-h_0\frac{\rho_1^2\rho_2^2+\rho_1^2\rho_3^2+\rho_2^2\rho_3^2}{\rho_1\rho_2\rho_3}\sin\theta,\\
&&\tau_0\frac{\partial \rho_1}{\partial t}=\mu \rho_1+h_0\rho_2\rho_3\cos\theta-g_1\rho_1^3-g_2(|\rho_2|^2+|\rho_3|^2)\rho_1,\nonumber\\
&&\tau_0\frac{\partial \rho_2}{\partial t}=\mu \rho_2+h_0\rho_1\rho_3\cos\theta-g_1\rho_2^3-g_2(|\rho_1|^2+|\rho_3|^2)\rho_2,\nonumber\\
&&\tau_0\frac{\partial \rho_3}{\partial t}=\mu \rho_3+h_0\rho_2\rho_1\cos\theta-g_1\rho_3^3-g_2(|\rho_2|^2+|\rho_1|^2)\rho_3,
\end{eqnarray}
where $\theta=\theta_1+\theta_2+\theta_3$, the first formula of the above equations shows that $h_0>0$, the solution of $\theta=0$ is stable
and $h_0<0$, the solution of $\theta=\pi$ is stable. Since we only consider the stable solution of the above equations, it follows that
\begin{eqnarray}
&&\tau_0\frac{\partial \rho_1}{\partial t}=\mu \rho_1+|h_0|\rho_2\rho_3-g_1\rho_1^3-g_2(|\rho_2|^2+|\rho_3|^2)\rho_1,\nonumber\\
&&\tau_0\frac{\partial \rho_2}{\partial t}=\mu \rho_2+|h_0|\rho_1\rho_3-g_1\rho_2^3-g_2(|\rho_1|^2+|\rho_3|^2)\rho_2,\nonumber\\
&&\tau_0\frac{\partial \rho_3}{\partial t}=\mu \rho_3+|h_0|\rho_2\rho_1-g_1\rho_3^3-g_2(|\rho_2|^2+|\rho_1|^2)\rho_3.
\end{eqnarray}
\subsection{Analysis of amplitude equation}
\noindent
\begin{theorem}
Different solutions of amplitude Eq. (42) can be classified as follows.
(i)~Homogeneous stationary state:
\begin{equation*}
\rho_1=\rho_2=\rho_3.
\end{equation*}
If $\mu<\mu_2=0$, then it is stable, and if $\mu>\mu_2=0$, then it is unstable.
(ii)~Stripe pattern:
\begin{equation*}
\rho_1=\sqrt{\frac{\mu}{g_1}}\neq0 ,\quad \rho_2=\rho_3=0.
\end{equation*}
If $\mu>\mu_3=\frac{h_0^2g_1}{(g_2-g_1)^2}$, then it is stable, and if $\mu<\mu_3=\frac{h_0^2g_1}{(g_2-g_1)^2}$, then it is unstable.
(iii)~Hexagonal pattern:
\begin{equation*}
\rho_1=\rho_2=\rho_3=\frac{|h_0|\pm\sqrt{h_0^2+4(g_1+2g_2)\mu}}{2(g_1+2g_2)}.
\end{equation*}
When the control parameter $\mu$ satisfies the condition of $\mu<\mu_4=\frac{(2g_1+g_2)h_0^2}{(g_2-g_1)^2}$, then the solution $\rho^-$ is always stable, and $\rho^+$ is always unstable, where
\begin{equation*}
\rho^-=\frac{|h_0|-\sqrt{h_0^2+4(g_1+2g_2)\mu}}{2(g_1+2g_2)},\quad \rho^+=\frac{|h_0|+\sqrt{h_0^2+4(g_1+2g_2)\mu}}{2(g_1+2g_2)}.
\end{equation*}
(iv)~Mixed structure state:
\begin{equation*}
\rho_1=\frac{|h_0|}{g_2-g_1}, \quad \rho_2=\rho_3=\sqrt{\frac{\mu-g_1\rho_1^2}{g_1+g_2}}.
\end{equation*}
When $\mu>\mu_3$ and $g_2>g_1$, then it is always unstable
\end{theorem}
Three stationary solutions are found based on stability analysis of the above amplitude equations. We can draw a conclusion about
the relationship between critical points, when $h_0\neq0$, $\mu_1<\mu_2<\mu_3<\mu_4$, where $\mu_1=\frac{-h_0^2}{4(g_1+2g_2)}$. When
parameter $\mu$ is increased to critical point $\mu_2=0$, instability of homogeneous stationary state will happen. Then a hexagon pattern begins to appear through nonequilibrium phase transition. When $h_0>0$, hexagonal pattern will appear. There exists a bistable region $\mu_1<\mu<\mu_2$, where homogeneous stationary state and hexagonal pattern are stable. Stripe pattern will appear at supercritical bifurcation when $\mu>\mu_3$, it is stable, and if $\mu>\mu_4$, hexagonal pattern is unstable. The system admits another bistable region $\mu_3<\mu<\mu_4$, where stripe pattern and hexagonal pattern are stable. When the control parameter is less than $\mu_3$, there is a transition from stripe pattern to hexagonal pattern, the system can transform hexagonal pattern to stripe pattern if control parameter is greater than $\mu_4$.
\section{Numerical simulations}
\noindent
In this section, we will perform numerical simulations to verify the results in Theorem 4.1 by using the mathematical software. Here, the system is studied on a rectangular grid of $200\times200$ points and the time step length is $\Delta t=0.02$ and the spatial step is set to be $\Delta h=1$. Initial conditions of Eq. (2) are chosen as
\begin{eqnarray*}
u(x,y,0)=u_*-0.002\times\zeta^*,\\
v(x,y,0)=v_*-0.002\times\zeta^*,
\end{eqnarray*}
where $\zeta^*$ is the random perturbation.
According to the previous theoretical analysis, the positive equilibrium $E_*$ is locally asymptotically stable when $\beta <\beta_H$
and $E_*$ is unstable if $\beta >\beta_H$. The Turing instability domain is set as $\beta_T<\beta <\beta_H$. Various patterns are presented with different values of parameter $\beta$ in the Turing instability domain and the values of parameter $a$. Difference is also exhibited between patterns with initial conditions under uniformly distributed random perturbation and those with symmetric initial conditions. Here only present the evolution of patterns of the prey, because patterns of the predators are similar.
\begin{table}[t!]
\doublerulesep 0.1pt
\tabcolsep 7.8mm
\centering
\caption{\rm Values of parameters are set in Figs.4-6. $(a=0.2, m=0.25, r=1, d_1=0.1, d_2=2.5)$.}
%
References
-
| [1] |
Hsu, S.B. and Huang, T.W. (1995), Global stability for a class of predator-prey systems,
SIAM Journal on Applied Mathematics, 55, 763-783.
|
-
| [2] |
Guan, X.N., Wang, W.M., and Cai, Y.L. (2011), Spatiotemporal dynamics of a Leslie-Gower predator-prey model incorporating a prey refuge,
Nonlinear Analysis. Real World Applications, 12, 2385-2395.
|
-
| [3] |
Camara, B.I. and Aziz-Alaoui, M.A. (2009), Turing and Hopf patterns formation in a predator-prey model with Leslie-Gower-type
functional response, Dynamics of Continuous, Discrete \& Impulsive System. Series B., 16, 479-488.
|
-
| [4] |
Holling, C.S. (1959), The components of predation as revealed by a study of small-mammal predation of the European pine sawfly,
Canadian Entomologist, 91, 293-320.
|
-
| [5] |
Holling, C.S. (1959), Some characteristics of simple types of predation and parasitism,
Canadian Entomologist, 91, 385-398.
|
-
| [6] |
Holling, C.S. (1965), The functional response of predator to prey density and its role in mimicry and population regulation,
Memoirs of the Entomological Society of Canada, 97, 5-60.
|
-
| [7] |
Ivlev V.S. (1961), Experimental Ecology of the Feeding of Fishes, Yale University Press: New Haven.
|
-
| [8] |
Hassell, M.P. and Varley, G.C. (1969), New inductive population model for insect parasites and its bearing on biological control,
Nature, 223, 1133-1137.
|
-
| [9] |
Arditi, R. and Saiah, H. (1992), Empirical evidence of the role of the heterogeneity in ratio-dependent consumption,
Ecology, 73, 1544-1551.
|
-
| [10] |
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.
|
-
| [11] |
Shi, H.B., Ruan, S., Su, Y., and Zhang, J.F. (2015), Spatiotemporal dynamics of a diffusive Leslie-Gower predator-prey model with ratio-dependent
functional response, International Journal of Bifurcation and Chaos, 25, 1530014.
|
-
| [12] |
Chang, X.Y. and Zhang, J.M. (2019), Dynamics of a diffusive Leslie-Gower predator-prey system with ratio-dependent Holling III functional response, Advances in Difference Equations, 76.
|
-
| [13] |
Hu, G.P. and Feng, Z.S. (2020), Turing instability and pattern formation in a strongly coupled diffusive predator-prey system,
International Journal of Bifurcation and Chaos, 30, 2030020.
|
-
| [14] |
Chen, M.X., Wu, R.C., and Chen, L.P. (2020), Spatiotemporal pattern induced by Turing and Turing-Hopf bifurcations in a predator-prey system,
Applied Mathematics and Compution, 380, 125300.
|
-
| [15] |
Abid, W., Yafia, R., Aziz-Alaoui, M.A., and Aghriche, A. (2018), Turing instability and Hopf bifurcation in a modified Leslie-Gower predator-prey model with cross-diffusion, International Journal of Bifurcation and Chaos, 28, 1850089.
|
-
| [16] |
Kar, T.K. and Ghosh, B. (2013), Impacts of maximum sustainable yield policy to prey-predator systems,
Ecological Modeling, 250, 134-142.
|
-
| [17] |
Huang, J.C., Gong, Y.J., and Ruan, S.G. (2013), Bifurcation analysis in a predator-prey model with constant-yield predator harvesting,
Discrete and Continuous Dynamical Systems. Series B, 18, 2101-2121.
|
-
| [18] |
Liu, C., Zhang, Q.L., and Zhang, Y. (2008), Bifurcation and control in a differential-algebraic harvested prey-predator model with stage
structure for predator, Interantional Journal of Bifurcation and Chaos, 18, 3159-3168.
|
-
| [19] |
Xiao, M. and Cao, J.D. (2009), Hopf bifurcation and non-hyperbolic equilibrium in a ratio-dependent predator-prey model with linear
harvesting rate: Analysis and computation, Mathematical and Computer Modelling, 50, 360-379.
|
-
| [20] |
Zhang, N., Chen, F., Su, Q., and Wu, T. (2011), Dynamic behaviors of a harvesting Leslie-Gower predator-prey model,
Discrete Dynamics in Nature and Society, 2011, 473949.
|
-
| [21] |
Zuo, W.Q., Ma, Zh.P., and Cheng, Zh.B. (2020), Spatiotemporal dynamics induced by Michaelis–Menten type prey harvesting in a diffusive
Leslie–Gower predator–prey model, International Journal of Bifurcation and Chaos, 30, 2050204.
|
-
| [22] |
Yavuz, M. and Sene, N. (2020), Stability analysis and numerical computation of the fractional predator–prey model with
the harvesting rate, Fractal and Fractional, 4, 35.
|
-
| [23] |
Turing, A.M. (1952), The chemical basis of morphogenesis,
Philosophical Transacations of the Royal Society of London. Series B, 237, 37-72.
|
-
| [24] |
Peng, Y.H. and Liu, Y.Y. (2016), Turing instability and Hopf bifurcation in a diffusive Leslie-Gower predator-prey model,
Mathematical Methods in Applied Sciences, 39, 4158-4170.
|
-
| [25] |
Han, R.J. and Dai, B.X. (2017), Cross-diffusion induced Turing instability and ampulitude equation for a
toxic-phytoplankton-zooplankton model with nonmonotonic functional response,
International Journal of Bifurcation and Chaos, 27(6), 1750088.
|
-
| [26] |
Chen, M.X., Wu, R.C., and Chen, L.P. (2020), Pattern dynamics in a diffusive Gierer-Meinhardt Model,
International Journal of Bifurcation and Chaos, 30, 2030035.
|
-
| [27] |
Ghorai, S. and Poria, S. (2016), Turing patterns induced by cross-diffusion in a predator-prey system in presence of habitat complexity,
Chaos, Soliton \& Fractals, 91, 421-429.
|
-
| [28] |
Song, Y.L. and Tang, X.S. (2017), Stability, steady-state bifurcations, and Turing patterns in a predator-prey model with herd behavior and prey-taxis, Studies in Applied Mathematics, 139, 371-404.
|
-
| [29] |
Wu, R., Shao, Y., Zhou, Y., and Chen, L. (2017), Turing and Hopf bifurcation of Gierer-Meinhardt activator-substrate model,
Electronic Journal of Differential Equations, 2017, 1-19.
|
-
| [30] |
Ouyang, Q. (2010), Nonlinear Science and Introduction of Pattern Dynamics, Beijing University Press: Beijing.
|
-
| [31] |
Wiggins, S. (2003), Introduction to Applied Nonlinear Dynamical Systems and Chaos, Springer-Verlag: New York.
|
