Turing Instability and Spatiotemporal Pattern Formation Induced by Nonlinear Reaction Cross-Diffusion in a Predator–Prey System with Allee Effect

Abstract

The Allee effect is widespread among endangered plants and animals in ecosystems, suggesting that a minimum population density or size is necessary for population survival. This paper investigates the stability and pattern formation of a predator–prey model with nonlinear reactive cross-diffusion under Neumann boundary conditions, which introduces the Allee effect. Firstly, the ODE system is asymptotically stable for its positive equilibrium solution. In a reaction system with self-diffusion, the Allee effect can destabilize the system. Then, in a reaction system with cross-diffusion, through a linear stability analysis, the cross-diffusion coefficient is used as a bifurcation parameter, and instability conditions driven by the cross-diffusion are obtained. Furthermore, we show that the system (5) has at least one inhomogeneous stationary solution. Finally, our theoretical results are illustrated with numerical simulations.

1. Introduction

Starting with the seminal research of Lotka and Volterra, the analysis of predator–prey model has become an important field of mathematical biology because of its practical and theoretical significance. Most articles about the simulation of predator–prey interaction use the logistic growth function. Meanwhile, some examples verify the existence of the Allee effect, which Allee proposed [1]. The Allee effect has many possible origins, e.g., difficulty looking for a partner, reproductive convenience, environmental regulation, or close relative decline. In general, for a population that is initially low density, if the per capita growth rate is an increasing function, then the species has an Allee effect. The Allee effect is divided into strong and weak; if the growth rate per capita is negative in a low-density limit, this is a strong Allee effect, while a weak Allee effect is expressed at zero density when the per capita growth rate is positive. A strong Allee effect creates a population threshold beyond which the population could grow. On the contrary, a group with a weak Allee effect has no threshold [2,3,4]. Prey is influenced by strong Allee effects. Numerous studies have revealed rich and diverse pattern-formation scenarios. Especially in a predator–prey system with a strong Allee effect on the population of the prey, space is vital for existence of species. Pattern formation can be observed both inside and outside the Turing field. Especially in reaction-diffusion systems, the function of space is significant for the existence of a species. On the other hand, numerous studies have shown that pattern formation is altered with the addition of a weak Allee effect. Currently, the Allee effect receives great attention from ecologist and mathematicians, because it can affect the population growth of a species [5,6,7,8,9,10,11,12].

In 1952, Turing proposed that space-diffusion can destroy interacting chemicals in a homogeneous and stable state, which is occurs through two coupled reaction–diffusion equations. This instability is known as diffusion-driven instability or Turing instability [13,14]. In recent years, reaction–diffusion systems have attracted increasing attention from mathematical biologists seeking insights into patterns occurring in ecosystems [15,16,17,18,19,20,21,22,23,24,25,26]. Many researchers have studied the emergence of patterns in population dynamics systems with self-diffusion [27,28,29,30,31]. While the reaction–diffusion system is not sterically stable, it exhibits Turing modes in space. With the broad growth of biomathematics and theoretical biology, the most prominent theme in the scopes of the dynamic system is the influence of cross-diffusion on mode appearance [32,33,34,35,36,37,38]. Biologically, integrated defenses against prey tend to move predators to areas with lower prey concentrations, which means that this can give rise to a spatially inhomogeneous population distribution. This can be described by reaction–diffusion systems with cross-diffusion. It has been seen that a appropriate cross-diffusion coefficient is able to generate this pattern [39,40,41].

In view of the above introduction about the Allee effect and the diffusion term, some scholars have studied the impact of the Allee effect on the stability and the pattern formation of the predator–prey system, and some scholars have studied the impact of diffusion on the stability and the pattern formation of the system. By studying the temporal and spatial distribution and dynamics of populations, human beings can better protect and adjust populations and rationally develop and utilize resources, which has far-reaching significance. Due to the aforementioned motivations, in this paper, we consider the effects of the Allee effect and nonlinear cross-diffusion terms on system stability and spatial mode formation in our model. The structure of the text is as follows: in Section 2, in terms of its mathematical and biological significance, we introduce the model. In Section 3, we discuss the influence of the Allee effect constant on the stability of the coexistence equilibrium of the predator–prey system with self-diffusion and ODE systems. Additionally, the conditions of cross-diffusion-driven Turing instability are obtained. In Section 4, we prove the existence of a nonconstant positive solution. In Section 5, our theoretical results are illustrated with numerical simulations. Finally, in Section 6, the conclusion is given.

2. The Mathematical Model

In the section, firstly, we introduce the Allee effect into the prey growth function.

In 1983, Verhulst–Pearl proposed the Verhulst–Pearl equation, which is also called the Logistic model:

$$\frac{dU}{dt}=r_{m}U(1-\frac{U}{K}),$$

(1)

$U(t)$ is species density, $K$ is carrying capacity, also known as environmental capacity, and it is a constant. $r_m$ is defined as the intrinsic growth rate. At this time, the average growth rate is $r=r_m(1-\frac{U}{K})$. When the population density reaches $K$, the birth rate is equal to the death rate.

The following model is used to describe the growth of a single population:

$$\frac{dU}{dt}\frac{1}{U}=e(1-\frac{U}{P})-g,$$

(2)

$e$ is the birth rate per capita when undisturbed by other individuals; when the actual birth rate per capita is zero, the density of the population is $P$. $g$ is the per capita mortality rate for adults. Equation (2) is equivalent to the logistic equation, and environmental carrying capacity is $K=P(e-g)/e$ [42]. Introducing the Allee effect to Equation (2):

$$\frac{dU}{dt}\frac{1}{U}=e(1-\frac{U}{P})(\frac{U}{U+m})-g,$$

(3)

$m$ is the Allee effect constant, and $U/U+m$ is known as the Allee effect, where the larger the $m$, the better the Allee effect. Solving equation $e(1-\frac{U}{P})(\frac{U}{U+m})-g=0$, the equivalent equation is $e{U}^2-P(e-g)U+Pgm=0$. Assuming ${\Delta }_0={[P(e-g)]}^2-4Pegm>0$, let $l=\frac{P(e-g)-\sqrt{{\Delta }_0}}{2e},$ $k=\frac{P(e-g)+\sqrt{{\Delta }_0}}{2e}$. Additionally, Equation (3) can also be expressed as:

$$\frac{dU}{dt}\frac{1}{U}=e(1-\frac{U}{P})(\frac{U}{U+m})-g=\frac{e(U-l)(k-U)}{P(U+m)},$$

(4)

Equation (4) can be simplified to: $\frac{du}{dt}=\frac{\omega u(1-u)(u-\beta )}{\sigma +u}$, where $u=\frac{U}{K},\omega =\frac{Ke}{P},\beta =\frac{l}{K},\sigma =\frac{m}{K}$.

Let $u_1, u_2$ denote the density of prey and predator. Let $y={\left(u_1, u_2\right)}^T$, $D(y)={(D_{ij}(y))}_{2\times 2}$. $y_t=\Delta [D(y)]+H(y)$ represents the strongly coupled reaction–diffusion system. Concerning the nonlinear cross-diffusion coefficient, the authors of [43,44] researched the systems for $D(y)=dig(c_1+\frac{k}{\epsilon +{u_3}^2},c_2,c_3)y$, $D(y)=dig(c_1+\frac{k}{\epsilon +{u_3}^2},c_2,c_3+c_3{c}_4{u}_1)y$. Additionally, here, we study the system for $D(y)=dig(c_1+\frac{k}{\eta +{u_2}^2},c_2+c_4{u}_1)y$. It can be simplified to $D(y)=dig(c_1+\frac{c_3}{1+{c_5{u}_2}^2},c_2+c_4{u}_1)y$, where $c_3=\frac{k}{\eta },c_5=\frac{1}{\eta }$.

In this article, we study a predator–prey model with the Allee effect and nonlinear cross-diffusion.

$$\{\begin{array}{l}\frac{\partial u_1}{\partial t}-\Delta (c_1+\frac{c_3}{1+{c_5{u}_2}^2})u_1=\frac{\omega u_1(1-u_1)(u_1-\beta )}{\sigma +u_1}-u_1{u}_2,\hspace{1em} t>0, x\in \Omega ,\\ \frac{\partial u_2}{\partial t}-c_2\Delta (1+c_4{u}_1)u_2=u_1{u}_2-\delta u_2,\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} t>0, x\in \Omega ,\\ \frac{\partial u_1}{\partial n}=\frac{\partial u_2}{\partial n}=0,\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} t>0, x\in \partial \Omega ,\\ u_i(x,0)=u_{i0}(x)\ge 0,\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}i=1,2, x\in \Omega ,\end{array}$$

(5)

$u_1(x)$ and $u_2(x)$ are the population densities of prey and predator, respectively. $\omega$ is the intrinsic prey growth rate, and δ represents the relative mortality of the predator. $\sigma$ is the Allee effect constant and $\beta$ is the Allee threshold. $c_1,c_2,c_3,c_4$ and $c_5$ are positive constants, the constants $c_m(m=1, 2)$ are termed self-diffusion coefficients, and $c_m(m=3, 4)$ are termed cross-diffusion coefficients. $\Omega$ is a bounded domain in $R^N$ with a smooth boundary $\partial \Omega$. The Homogeneous Neumann boundary condition states the population flux across the boundary is zero.

The diffusive flux of $u_1$ is

$$J=-\nabla \left(c_1{u}_1+\frac{c_3{u}_1}{1+c_5{u}_2^2}\right)=-\left(c_1+\frac{c_3}{1+c_5{u}_2^2}\right)\nabla u_1+\frac{2c_3{c}_5{u}_1{u}_2}{{(1+c_5{u}_2^2)}^2}\nabla u_2,$$

as $2c_3{u}_1{u}_2/{(1+c_5{u}_2^2)}^2\ge 0$, the $2c_3{u}_1{u}_2/{(1+c_5{u}_2^2)}^2\nabla u_2$ part of the flux is in the direction of the increasing population density of $u_2$.

The diffusive flux of $u_2$ is

$$J=-\nabla (c_2{u}_2+c_2{c}_4{u}_1{u}_2)=-c_2{c}_4{u}_2\nabla u_1-(c_2+c_2{c}_4{u}_1)\nabla u_2,$$

as $-c_2{c}_4{u}_2<0$, $-c_2{c}_4{u}_2\nabla u_1$ of the flux is in the direction of the decreasing population density of $u_1$.

As predators hunt their prey, the flux should point in the direction of increasing prey population density. However, in certain ecosystems, large numbers of prey species come together to form a large group to protect themselves from predators.

3. Stability of the Coexisting Equilibrium

In this section, we discuss the effect of the Allee effect constant on the stability of the coexistence equilibrium point of ODE systems and predator–prey systems with self-diffusion. When $(1-\delta )(\delta -\beta )>0$, $\tilde{y}=\left({\tilde{u}}_1, {\tilde{u}}_2\right)$ is the unique positive coexistence equilibrium point in system (5), where

$${\tilde{u}}_1=\delta ,{\tilde{u}}_2=\frac{\omega (1-\delta )(\delta -\beta )}{\sigma +\delta }.$$

(6)

Consider the ODE system of system (5) as follows:

$$\{\begin{array}{l}\frac{d{u}_1}{dt}=\frac{\omega u_1(1-u_1)(u_1-\beta )}{\sigma +u_1}-u_1{u}_2, t>0,\\ \frac{d{u}_2}{dt}=u_1{u}_2-\delta u_2, t>0.\end{array}$$

(7)

Theorem 1.

If the parameters of (5) satisfy

$$\sigma >\frac{\beta -{\delta }^2}{2\delta -\beta -1}$$

(8)

$\tilde{u}$ is locally asymptotically stable in system (7).

Proof. 

Indicate

$$\begin{array}{l}H(y)=\left(\begin{array}{c}{H}_1(y)\\ H_2(y)\end{array}\right)=\left(\begin{array}{c}{u}_1{g}_1(y)\triangleq u_1\left(\frac{\omega (1-u_1)(u_1-\beta )}{\sigma +u_1}-u_2\right)\\ u_2{g}_2(y)\triangleq u_2(u_1-\delta )\end{array}\right)\\ =\left(\begin{array}{c}{u}_1\left(\frac{\omega (-{u_1}^2+(\beta +1)u_1-\beta )}{\sigma +u_1}-u_2\right)\\ u_2(u_1-\delta )\end{array}\right).\end{array}$$

(9)

Noticing $H(\tilde{y})=0$,

$$H_u(\tilde{y})=\left(\begin{array}{cc}{\tilde{u}}_1\left(\frac{-\omega {{\tilde{u}}_1}^2-2\omega \sigma {\tilde{u}}_1+\omega \sigma \beta +\omega \sigma +\omega \beta }{{(\sigma +{\tilde{u}}_1)}^2}\right)& -{\tilde{u}}_1\\ {\tilde{u}}_2& 0\end{array}\right).$$

(10)

The characteristic polynomial of $H_y(\tilde{y})$ is

$$\chi (\eta )={\eta }^2+A_1\eta +A_2,$$

where $A_1=Trace(H_y(\tilde{y})),$ $A_2=Det(H_y(\tilde{y}))$.

If (8) holds, $A_1$ and $A_2$ are positive, which is easy to verify. So, due to the Routh–Hurwitz criterion, a conclusion is established. □

Next, the following system (5) with self-diffusion can be studied:

$$\{\begin{array}{l}\frac{\partial u_1}{\partial t}-c_1\Delta u_1=\frac{\omega u_1(1-u_1)(u_1-\beta )}{\sigma +u_1}-u_1{u}_2, t>0, x\in \Omega ,\\ \frac{\partial u_2}{\partial t}-c_2\Delta u_2=u_1{u}_2-\delta u_2,\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}t>0, x\in \Omega ,\\ \frac{\partial u_1}{\partial n}=\frac{\partial u_2}{\partial n}=0,\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} t>0, x\in \partial \Omega ,\\ u_i(x,0)=u_{i0}(x)\ge 0,\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}i=1,2, x\in \Omega .\end{array}$$

(11)

To research the local asymptotic stability of Equation (11), similar to [45], we give Notation 1.

Notation 1.

Let$0={\pi }_1<{\pi }_2<\cdots \to \infty$be the eigenvalues of$-\Delta$on$\Omega$ under no-flux boundary conditions, and let$E({\pi }_m)$be the space of eigenfunctions corresponding to ${\pi }_m$. Then,

(i)

${\mathbf{X}}_{mn}:=\left\{\mathbf{c}\cdot {\varphi }_{mn}:\mathbf{c}\in R^3\right\}$,$\left\{{\varphi }_{mn}\right\}$are the orthonormal bases of$E({\pi }_m)$,$n=1,2\dots \dim E({\pi }_m)$.

(ii)

$\mathbf{X}:=\left\{\mathbf{u}\in {\left[C^1(\overline{\Omega })\right]}^3:{\partial }_n{u}_1={\partial }_n{u}_2=0 on \partial \Omega \right\}$, then$\mathbf{X}={\oplus }_{m=1}^{\infty }{\mathbf{X}}_m$,${\mathbf{X}}_m={\oplus }_{n=1}^{\dim E({\pi }_m)}{\mathbf{X}}_{mn}$.

Theorem 2.

If (8) holds, $\tilde{y}$ is locally asymptotically stable in system (11).

Proof. 

Linearizing system (11), then

$$y_t=(D\Delta +H_y(\tilde{y}))y$$

(12)

where $D=\left(\begin{array}{cc}{c}_1& 0\\ 0& c_2\end{array}\right)$.

Due to Notation 1, if and only if ${\mathbf{X}}_m$ is an eigenvalue of the matrix $-{\pi }_{m}D+H_y(\tilde{y})$, it is invariant under operator $D\Delta +H_y(\tilde{y})$, when this operator on ${\mathbf{X}}_m$, $\lambda$ is an eigenvalue. ${\psi }_m(\lambda )=\left|\lambda E-(-{\pi }_{m}D+H_y(\pi ))\right|={\lambda }^2+B_1\lambda +B_2$ is the characteristic polynomial of $-{\pi }_{m}D+H_y(\tilde{y})$, where

$$B_1={\pi }_m{c}_1+{\pi }_m{c}_2+{\tilde{u}}_1(\frac{\omega {\tilde{u}}_1{}^2+2\omega \sigma {\tilde{u}}_1-(\omega \beta \sigma +\omega \sigma +\omega \beta )}{{(\sigma +{\tilde{u}}_1)}^2}),$$

$$B_2={\tilde{u}}_1{\tilde{u}}_2+{\pi }_m{c}_2({\pi }_m{c}_1+{\tilde{u}}_1(\frac{\omega {\tilde{u}}_1{}^2+2\omega \sigma {\tilde{u}}_1-(\omega \beta \sigma +\omega \sigma +\omega \beta )}{{(\sigma +{\tilde{u}}_1)}^2})).$$

If (8) holds, we are able to show that $B_1$ and $B_2$ are positive. Then, for each $i\ge 1$, because of the Routh–Hurwitz criterion, the two roots ${\lambda }_{i1}, {\lambda }_{i2}\ge 0$ of ${\psi }_i\left(\lambda \right)=0$ have negative real parts. So, the conclusion holds. □

Theorem 3.

Assume that

$$\frac{2c_5{\tilde{u}}_2}{(c_2+c_2{c}_4{\tilde{u}}_1)(1+c_5{\tilde{u}}_2{}^2)+2c_5{c}_2{c}_4{\tilde{u}}_2}>0$$

(13)

and (8) hold, if${\pi }_2<\tilde{\mu },$where${\pi }_2$can be seen in Notation 1, $\tilde{\mu }$ will be seen in (14). There must be a positive constant$c_3{}^{\ast },$when the $c_3\ge c_3{}^{\ast },$coexistence equilibrium solution$\tilde{y}$of system (5) is unstable.

Proof. 

We indicate that

$$\vartheta (y)={\left(c_1{u}_1+\frac{c_3{u}_1}{1+c_5{u}_2{}^2},c_2{u}_2+c_2{c}_4{u}_1{u}_2\right)}^T$$

Linearizing system (5) at $\tilde{y}$, then we have $y_t=({\vartheta }_y\Delta +H_y(\tilde{y}))y$, where

$${\vartheta }_y=\left(\begin{array}{l}\begin{array}{cc}{c}_1+\frac{c_3}{1+c_5{\tilde{u}}_2{}^2}& \frac{-2c_3{c}_5{\tilde{u}}_1{\tilde{u}}_2{}^2}{{(1+c_5{\tilde{u}}_2{}^2)}^2}\end{array}\\ \begin{array}{cc}{c}_2{c}_4{\tilde{u}}_2& \begin{array}{cc}\begin{array}{cc}& \end{array}& c_2+c_2{c}_4{\tilde{u}}_1\end{array}\end{array}\end{array}\right),$$

the characteristic polynomial of $-{\pi }_m{\vartheta }_{\mathrm{y}}+H_y(\tilde{y})$ is ${\varsigma }_m(\lambda )={\lambda }^2+{\overline{B}}_1\lambda +{\overline{B}}_2,$

$$\begin{array}{l}{\overline{B}}_1={\pi }_m(c_1+\frac{ c_3}{1+c_5{\tilde{u}}_2{}^2}+c_2+c_2{c}_4{\tilde{u}}_1)-{\tilde{u}}_1(\frac{-\omega {\tilde{u}}_1{}^2-2\omega \sigma {\tilde{u}}_1+\omega \sigma \beta +\omega \sigma +\omega \beta }{{(\sigma +{\tilde{u}}_1)}^2}),\\ {\overline{B}}_2={\pi }_m{}^2((c_1+\frac{c_3}{1+c_5{\tilde{u}}_2{}^2})\left(c_2+c_2{c}_4{\tilde{u}}_1\right)+\frac{2c_5{c}_2{c}_3{c}_4{\tilde{u}}_2}{(1+c_5{\tilde{u}}_2{}^2)})\\ \hspace{1em}\hspace{1em}+{\pi }_m(-{\tilde{u}}_1(\frac{-\omega {\tilde{u}}_1{}^2-2\omega \sigma {\tilde{u}}_1+\omega \sigma \beta +\omega \sigma +\omega \beta }{{(\sigma +{\tilde{u}}_1)}^2})\left(c_2+c_2{c}_4{\tilde{u}}_1\right)-c_2{c}_4{\tilde{u}}_1{\tilde{u}}_2-\frac{2c_3{c}_5{\tilde{u}}_2}{(1+c_5{\tilde{u}}_2{}^2)})+{\tilde{u}}_1{\tilde{u}}_2.\end{array}$$

Let ${\lambda }_1({\pi }_m)$, ${\lambda }_2({\pi }_m)$ be the two roots of ${\varsigma }_m(\lambda )=0$, and ${\lambda }_1({\pi }_m)\cdot {\lambda }_2\left({\pi }_m\right)=-{\overline{B}}_2$. To have at least one $\mathrm{Re}({\lambda }_m)>0$, ${\overline{B}}_2<0$ is enough. Next, we find the condition for ${\overline{B}}_2<0$. ${\overline{B}}_2=\det ({\pi }_m{\Phi }_y-H_y(\tilde{y}))$ and ${\overline{B}}_2=Q_2{\pi }_m{}^2+Q_1{\pi }_m+Q_0$, where

$$Q_2=(c_1+\frac{c_3}{1+c_5{\tilde{u}}_2{}^2})(c_2+c_2{c}_4{\tilde{u}}_1)+\frac{2c_5{c}_2{c}_3{c}_4{\tilde{u}}_2}{(1+c_5{\tilde{u}}_2{}^2)},$$

$$Q_1=-{\tilde{u}}_1(\frac{-\omega {\tilde{u}}_1{}^2-2\omega \sigma {\tilde{u}}_1+\omega \sigma \beta +\omega \sigma +\omega \beta }{{(\sigma +{\tilde{u}}_1)}^2})\left(c_2+c_2{c}_4{\tilde{u}}_1\right)-c_2{c}_4{\tilde{u}}_1{\tilde{u}}_2-\frac{2c_3{c}_5{\tilde{u}}_2}{(1+c_5{\tilde{u}}_2{}^2)},$$

$$Q_0={\tilde{u}}_1{\tilde{u}}_2,\det (H_y)={\tilde{u}}_1{\tilde{u}}_2.$$

Let $\tilde{Q}(\pi )=Q_2{\pi }^2+Q_1\pi +\det (H_y)$ and ${\tilde{\pi }}_1$, ${\tilde{\pi }}_2$, two roots of $\tilde{Q}(\pi )=0$, satisfy $\mathrm{Re}({\tilde{\pi }}_1)\le \mathrm{Re}({\tilde{\pi }}_2)$. Then, $\det (H_y)={\tilde{\pi }}_1{\tilde{\pi }}_2>0.$

$$\underset{c_3\to \infty }{\lim }\frac{Q_2}{c_3}=(\frac{1}{1+c_5{\tilde{u}}_2{}^2})\left(c_2+c_2{c}_4{\tilde{u}}_1\right)+\frac{2c_5{c}_2{c}_4{\tilde{u}}_2}{(1+c_5{\tilde{u}}_2{}^2)},\underset{c_3\to \infty }{\lim }\frac{Q_1}{c_3}=-\frac{2c_5{\tilde{u}}_2}{(1+c_5{\tilde{u}}_2{}^2)}.$$

Note that

$$\begin{gathered} \underset{c_3\to \infty }{\lim }\frac{\tilde{Q}(\pi )}{c_3}=\underset{c_3\to \infty }{\lim }\frac{Q_2{\pi }^2+Q_1\pi +\det \left(H_u\right)}{c_3} \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=((\frac{1}{1+c_5{u}_2{}^2})\left(c_2+c_2{c}_4{\tilde{u}}_1\right)+\frac{2c_5{c}_2{c}_4{\tilde{u}}_2}{(1+c_5{\tilde{u}}_2{}^2)}){\pi }_m{}^2+(-\frac{2c_5{\tilde{u}}_2}{(1+c_5{\tilde{u}}_2{}^2)}){\pi }_m. \end{gathered}$$

Under condition (13), coherence contention is demonstrated; when $c_3$ is big enough, $\mathrm{Re}({\tilde{\pi }}_1)>0$. Moreover, as $\mathrm{Re}({\tilde{\pi }}_1)\mathrm{Re}({\tilde{\pi }}_2)>0$, $\mathrm{Re}({\tilde{\pi }}_2)>0$, and

$$\begin{array}{l}\underset{c_3\to \infty }{\lim }{\tilde{\pi }}_1=\frac{2c_5{\tilde{u}}_2}{(c_2+c_2{c}_4{\tilde{u}}_1)(1+c_5{\tilde{u}}_2{}^2)+2c_5{c}_2{c}_4{\tilde{u}}_2}\triangleq \tilde{\mu }\\ \underset{c_3\to \infty }{\lim }{\tilde{\pi }}_2=0,\end{array}$$

(14)

So, there must be a constant ${c_3}^{\ast }$, when $c_3>{c_3}^{\ast }$, the following holds: $\tilde{Q}(\pi )<0, \pi \in ({\tilde{\pi }}_2,{\tilde{\pi }}_1).$ Since that $0<{\pi }_2<\tilde{\pi },$ that is

$$\frac{2c_5{\tilde{u}}_2}{(d_2+d_2{d}_4{\tilde{u}}_1)(1+c_5{\tilde{u}}_2{}^2)+2c_5{d}_2{d}_4{\tilde{u}}_2}>0.$$

We have ${\pi }_2\in \left({\tilde{\pi }}_2,{\tilde{\pi }}_1\right),$ so it follows that $\tilde{Q}({\pi }_2)<0.$ Thus, we show that ${\overline{B}}_2<0,$ and the proof is completed. □

It can be seen from the above theorem that cross-diffusion makes the positive equilibrium solution unstable.

4. Nonhomogeneous Steady States

In this section, our research concerns inhomogeneous positive solutions from the Leray–Schauder theory [45].

Firstly, we study the following form of the corresponding steady state of system (5):

$$\{\begin{array}{l}-\Delta \left(c_1{u}_1+\frac{c_3{u}_1}{1+c_5{u}_2{}^2}\right)=\frac{\omega u_1(1-u_1)(1-\beta )}{\sigma +u_1}-u_1{u}_2, x\in \Omega ,\\ -c_2\Delta \left(u_2+c_4{u}_1{u}_2\right)=u_1{u}_2-\delta u_2,\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} x\in \Omega .\end{array}$$

(15)

Domain Ω determines the generic constants $S^{\ast }, \stackrel{\_}{\mathrm{S}}, \underset{\_}{\mathrm{S}}$. A collection is defined: $B(S)=\left\{y={(u_1,u_2)}^T:\underset{\_}{S}<u_i<\stackrel{\_}{S},i=1,2\right\},$ then in the set $B(S)$, we can obtain nonhomogeneous positive solutions.

Let $\vartheta (y)={(c_1{u}_1+\frac{c_3{u}_1}{1+c_5{u}_2^2},c_2{u}_2+c_4{u}_1{u}_2)}^T$. Next, (15) is equivalent to

$$\{\begin{array}{l}-\Delta \vartheta (y)=H(y), t>0, x\in \Omega ,\\ {\partial }_{n}y=0, x\in \partial \Omega .\end{array}$$

(16)

For $y$, ${\vartheta }_y(y)$ exists and ${\vartheta }_y(y)>0$, the determinant of ${\vartheta }_y(y)>0$. So, $y>0$ to (22) when

$$R(y):={(I-\Delta )}^{-1}\left\{{\vartheta }_y^{-1}\left[H(y)+\nabla y\cdot {\Phi }_{yy}(y)\cdot \nabla y\right]+y\right\}=y,$$

${(I-\Delta )}^{-1}$ is the inverse of $I-\Delta$ under no-flux boundary constraints. Observe

$$R_y(\tilde{y})={(I-\Delta )}^{-1}\left\{{\Phi }_y^{-1}(\tilde{y})H_{\mathbf{y}}(\tilde{y})+I\right\},$$

For the problem of eigenvalue

$$\{\begin{array}{l}-(I-R_{\mathbf{y}}(\tilde{y}))\Psi =\lambda \Psi , x\in \Omega ,\\ {\partial }_{\mathrm{n}}\Psi =0, x\in \partial \Omega ,\end{array}$$

(17)

$\Psi ={({\Psi }_1,{\Psi }_2,{\Psi }_3)}^T$. Since $R(\cdot )$ has no fixed points in $\partial B(S)$, the Leray–Schauder degree is well defined. Utilizing the Leray–Schauder theorem,

$$index(I-R(\cdot ),\tilde{y})={(-1)}^{\tau }, \tau =\underset{\xi >0}{\Sigma }{n}_{\xi },$$

(18)

If zero is not the eigenvalue of (17), $n_{\xi }$ is the multiplicity of $\xi$.

If and only if $\xi$ is an eigenvalue of the matrix $-I+\frac{1}{1+{\pi }_m}({\vartheta }_y^{-1}\tilde{y}{H}_{\mathbf{y}}(\tilde{y})+I)=\frac{1}{1+{\pi }_m}(-{\pi }_{m}I+{\vartheta }_y^{-1}(\tilde{y})H_{\mathbf{y}}(\tilde{y})),$ $m \ge 1$, $1\le n\le \dim E({\pi }_m)$, ${\mathbf{X}}_{mn}$ is invariant under $-I+D_{\mathbf{y}}f({\mathbf{y}}^{*})$. The number of positive eigenvalues $\xi$ in ${\mathbf{X}}_{mn}$ is odd if and only if $\det (-{\pi }_{m}I+{\vartheta }_y^{-1}(\tilde{y})H_{\mathbf{y}}(\tilde{y}))>0$. Set

$$L(\pi )=L(\tilde{y},\pi )\triangleq \det (-\pi I+{\vartheta }_y^{-1}(\tilde{y})H_{\mathbf{y}}(\tilde{y})),$$

(19)

and next, we obtain some consequences.

Proposition 1.

Assume that all $m\ge 1,$matrix$-{\pi }_{m}I+{\pi }_y^{-1}(\tilde{y})H_{\mathbf{y}}(\tilde{y})$is nonsingular. So,$index(I-R(\cdot ),\tilde{y})={(-1)}^{\tau },$ $\tau =\underset{m\ge 1,H({\pi }_m)>0}{\Sigma }\dim E({\pi }_m)$.

For our convenience to compute$index(I-R(\cdot ),\tilde{y})$, next, we study the sign of$L({\pi }_m)$, fix$c_1$,$c_2$,$c_4$, and consider the dependence of$L({\pi }_m)$on$c_3$. Then, we note

$$L(\pi )=\det \{{\vartheta }_y^{-1}(\tilde{y}))\det (\pi {\vartheta }_{\mathbf{y}}\tilde{y}-H_{\mathbf{y}}(\tilde{y})\}.$$

Now that we have determined$\det \{\pi {\vartheta }_y(\tilde{y})\}>0$, we just need to study$\det \{\pi {\Phi }_{\mathbf{y}}\tilde{y}{H}_{\mathbf{y}}(\tilde{y})\}$. Actually, the value of obtained fromSection 3is equal to$\det \{\pi {\vartheta }_{\mathbf{y}}\tilde{y}-H_{\mathbf{y}}(\tilde{y})\}$. We obtain enough conditions to${\overline{B}}_2<0$during the proof of Theorem 3.

Proposition 2.

Suppose (3.8) holds. There must be a constant$c_3^{*}$, for all$c_3\ge c_3^{*}$, two roots${\tilde{\pi }}_1$,${\tilde{\pi }}_2$of$\det \{\pi {\vartheta }_{\mathbf{y}}(\tilde{y})-H_{\mathbf{y}}(\tilde{y})\}=0$are real and satisfy (3.9). Furthermore,

$$\{\begin{array}{l}0<{\tilde{\pi }}_1<{\tilde{\pi }}_2,\\ \det \{\pi {\vartheta }_y^{}\tilde{y}-H_{\mathbf{y}}(\tilde{y})\}>0, \pi \in (-\infty ,{\tilde{\pi }}_2)\cup ({\tilde{\pi }}_1,\infty ),\\ \det \{\pi {\vartheta }_y^{}\tilde{y}-H_{\mathbf{y}}(\tilde{y})\}<0, \pi \in ({\tilde{\pi }}_2,{\tilde{\pi }}_1).\end{array}$$

(20)

We study the influence of cross-diffusion on the positive solution of (16); firstly, we conclude that there is no nonconstant equilibrium solution without cross-diffusion.

Theorem 4.

Assume that

$$c_2\ge \frac{-\delta }{{\pi }_2},c_3=0,c_4=0,$$

(21)

which${\pi }_2$can see in Notation 1. Next, there must be a positive constant${\overline{D}}_1$; when$c_1\ge {\overline{D}}_1$, (15) has no nonconstant positive solution; moreover,

$$u_i\equiv {\overline{u}}_i, i=1,2 \mathrm{where} {\overline{u}}_i\triangleq \frac{1}{|\Omega |}{\int }_{\Omega }{u}_{i}dx.$$

Proof. 

Let ${\overline{u}}_i\triangleq \frac{1}{|\Omega |}{\int }_{\Omega }{u}_{i}dx$, $i=1,2$. Multiplying the ith equation of (15) through $u_i-{\overline{u}}_i$, then integrating the results over $\Omega$ by parts yields

$$\begin{array}{l}{c}_1{\int }_{\Omega }{\left|\nabla u_1\right|}^2={\int }_{\Omega }(u_1-{\overline{u}}_1)(u_1{g}_1(u_1,u_2)-{\overline{u}}_1{g}_1({\overline{u}}_1,{\overline{u}}_2))\mathrm{d}x\\ ={\int }_{\Omega }(u_1-{\overline{u}}_1)(u_1(\frac{\omega (1-u_1)(u_1-\beta )}{\sigma +u_1}-u_2)-{\overline{u}}_1(\frac{\omega (1-{\overline{u}}_1)({\overline{u}}_1-\beta )}{\sigma +{\overline{u}}_1}-{\overline{u}}_2))\mathrm{d}x\\ ={\int }_{\Omega }{(u_1-{\overline{u}}_1)}^2[\frac{(\omega \sigma +\omega \sigma \beta )(u_1+{\overline{u}}_1)+(\omega \beta -\omega +\omega (u_1+{\overline{u}}_1)u_1\cdot {\overline{u}}_1)-\omega \sigma \beta -\omega \sigma ({\overline{u}}_1{}^2+{\overline{u}}_1{u}_1+u_1^2)}{(\sigma +u_1)(\sigma +{\overline{u}}_1)}]\mathrm{d}x\\ +{\int }_{\Omega }(u_1-{\overline{u}}_1)({\overline{u}}_1{\overline{u}}_2-u_1{u}_2)\mathrm{d}x,\end{array}$$

(22)

$$\begin{array}{l}{c}_2{\int }_{\Omega }{\left|\nabla u_2\right|}^2={\int }_{\Omega }(u_2-{\overline{u}}_2)(u_2{g}_2(u_1,u_2)-{\overline{u}}_2{g}_2({\overline{u}}_1,{\overline{u}}_2))\mathrm{d}x\\ ={\int }_{\Omega }(u_2-{\overline{u}}_2)(u_2(u_1-\delta )-{\overline{u}}_2({\overline{u}}_1-\delta ))\mathrm{d}x\\ ={\int }_{\Omega }{(u_1-{\overline{u}}_1)}^2(-\delta )dx+{\int }_{\Omega }(u_2-{\overline{u}}_2)(u_2{u}_1-{\overline{u}}_2{\overline{u}}_1)\mathrm{d}x. \end{array}$$

(23)

$$\begin{array}{l}{ \mathrm{c}}_1{\int }_{\Omega }{\left|\nabla u_1\right|}^2+c_2{\int }_{\Omega }{\left|\nabla u_2\right|}^2\\ ={\int }_{\Omega }{(u_1-{\overline{u}}_1)}^2[\frac{(\omega \sigma +\omega \sigma \beta )(u_1+{\overline{u}}_1)+(\omega \beta -\omega +\omega (u_1+{\overline{u}}_1)u_1{\overline{u}}_1)-\omega \sigma \beta -\omega \sigma ({{\overline{u}}_1}^2+{\overline{u}}_1{u}_1+u_1^2)}{(\sigma +u_1)(\sigma +{\overline{u}}_1)}]\mathrm{d}x\\ +{\int }_{\Omega }{(u_1-{\overline{u}}_1)}^2(-\delta )dx+{\int }_{\Omega }(u_1-{\overline{u}}_1)(u_2-{\overline{u}}_2)(\frac{{\overline{u}}_1{\overline{u}}_2-u_1{u}_2}{u_2-{\overline{u}}_2}+\frac{u_1{u}_2-{\overline{u}}_1{\overline{u}}_2}{u_1-{\overline{u}}_1})\mathrm{d}x.\end{array}$$

(24)

By the Cauchy inequality,

$$\begin{array}{l}{ \mathrm{c}}_1{\int }_{\Omega }{\left|\nabla u_1\right|}^2+c_2{\int }_{\Omega }{\left|\nabla u_2\right|}^2\\ \le {\int }_{\Omega }{(u_1-{\overline{u}}_1)}^2[\frac{(\omega \sigma +\omega \sigma \beta )(u_1+{\overline{u}}_1)+(\omega \beta -\omega +\omega (u_1+{\overline{u}}_1)u_1{\overline{u}}_1)-\omega \sigma \beta -\omega \sigma ({{\overline{u}}_1}^2+{\overline{u}}_1{u}_1+u_1^2)}{(\sigma +u_1)(\sigma +{\overline{u}}_1)}]\mathrm{d}x\\ +{\int }_{\Omega }{(u_1-{\overline{u}}_1)}^2(-\delta +\theta )dx+{\int }_{\Omega }{(u_1-{\overline{u}}_1)}^2\frac{1}{4\theta }{[(\frac{{\overline{u}}_1{\overline{u}}_2-u_1{u}_2}{u_2-{\overline{u}}_2}+\frac{u_1{u}_2-{\overline{u}}_1{\overline{u}}_2}{u_1-{\overline{u}}_1})]}^2\mathrm{d}x,\end{array}$$

(25)

where $\theta$ is an arbitrary positive constant.

Due to Poincare Inequality,

$$c_1{\int }_{\Omega }{|\nabla u_1|}^2+c_2{\int }_{\Omega }{|\nabla u_2|}^2\ge {\int }_{\Omega }{c}_1{\mu }_2{(u_1-{\overline{u}}_1)}^2\mathrm{d}x+{\int }_{\Omega }{c}_2{\mu }_2{(u_2-{\overline{u}}_2)}^2\mathrm{d}x,$$

(26)

So there is a ${\theta }_0$ such that

$$\begin{array}{l}{\mu }_2{d}_2\ge -\delta +{\theta }_0,\\ {\mu }_2{D}_1\ge [\frac{(\omega \sigma +\omega \sigma \beta )(u_1+{\overline{u}}_1)+(\omega \beta -\omega +\omega (u_1+{\overline{u}}_1)u_1{\overline{u}}_1)-\omega \sigma \beta -\omega \sigma ({{\overline{u}}_1}^2+{\overline{u}}_1{u}_1+u_1^2)}{(\sigma +u_1)(\sigma +{\overline{u}}_1)}]\mathrm{d}x\\ +\frac{1}{4\theta }{[(\frac{{\overline{u}}_1{\overline{u}}_2-u_1{u}_2}{u_2-{\overline{u}}_2}+\frac{u_1{u}_2-{\overline{u}}_1{\overline{u}}_2}{u_1-{\overline{u}}_1})]}^2.\end{array}$$

(27)

We can conclude that $u_i={\overline{u}}_i, i=1,2$. □

Theorem 5.

Let $c_1,c_2,c_4$be fixed and satisfy (8) and (11), and let$\tilde{\mu }$be defined in (12). If$\tilde{\mu }\in ({\mu }_n,{\mu }_{n+1}),n\ge 1$and${\sigma }_n={\displaystyle \sum _{m=2}^n\dim E({\pi }_m)}$are odd, there must be a positive number$c_3^{*}$; if$c_3\ge c_3^{*}$, system (5) has at least a nonhomogeneous positive steady-state solution.

Proof. 

Let ${\overline{D}}_1,D_1,D_2,D_3,D_4$ be fixed and meet (14); in Proposition 2, there must be a positive constant $c_3^{*}$; when $c_3\ge c_3^{*}$, (20) holds and

$$0={\pi }_1<{\tilde{\pi }}_2<{\pi }_2,{\tilde{\pi }}_1\in ({\pi }_n, {\pi }_{n+1}) .$$

(28)

We show when $c_3\ge c_3^{*}$, (15) has at least an inhomogeneous positive solution. It is obtained because of the homotopy invariance of topological degree. Instead, we assume this protestation is not correct. Set

$$\vartheta (t;\mathbf{u})={(t(c_1-{\widehat{c}}_1)(c_1{u}_1+\frac{c_3{u}_1}{1+c_5{u}_2}),{\widehat{c}}_2+t(c_2-{\widehat{c}}_2))(u_2+c_4{u}_1{u}_2))}^T,t\in [0,1]$$

where ${\widehat{c}}_1={\overline{D}}_1, {\widehat{c}}_2=-\frac{\delta }{{\pi }_2}$, and the study issues

$$\{\begin{array}{l}-\Delta \vartheta (t;y)=H(y), t\in [0,1], x\in \Omega ,\\ {\partial }_{n}y=0, x\in \partial \Omega .\end{array}$$

(29)

$\tilde{y}$ is the unique constant positive solution of (29). Next, $y$ is a positive solution of (29) when

$$R(t;y):={(I-\Delta )}^{-1}\left\{{\vartheta }_y^{-1}\left[H(y)+\nabla y\cdot {\Phi }_{yy}(t;y)\cdot \nabla y\right]+y\right\}=y.$$

By the definition of $R(y),R(1;y)=R(y).$ Theorem 4 indicates that $R(0;y)=y$ has only one positive constant solution in $B(S)$. By a direct computation,

$$R_y(t;\tilde{y})={(I-\Delta )}^{-1}\left\{{\vartheta }_y(t;\tilde{y})H_y(\tilde{y})+I\right\}.$$

Especially

$$\begin{array}{l}{R}_y(0;\tilde{y})={(I-\Delta )}^{-1}\left\{D_{}^{-1}{H}_y(\tilde{y})+I\right\},\\ R_y(1;\tilde{y})={(I-\Delta )}^{-1}\left\{{\vartheta }_y^{-1}(\tilde{y})H_y(\tilde{y})+I\right\}{=\mathrm{R}}_y(\tilde{y}),\end{array}$$

$D=diag({\widehat{c}}_1,{\widehat{c}}_2,{\widehat{c}}_4)$. Given Propositions 4.2 and (28), the following conclusions can be drawn:

$$\{\begin{array}{l}L({\pi }_1)=L(0)>0,\\ L({\pi }_m)<0, 2\le i\le n,\\ L({\pi }_m)>0, i\ge n+1.\end{array}$$

So, for all $m\ge 1$, zero is not an eigenvalue of the matrix ${\pi }_{m}I-{\vartheta }_m^{-1}(\tilde{y})H_y(\tilde{y})$ and ${\sigma }_n={\displaystyle \sum _{m=2}^n\dim E({\pi }_m)}={\displaystyle \sum _{m\ge 1,H({\pi }_m)<0}^n\dim E({\pi }_m)}$, which is odd. Because of Proposition 1, we have

$$index(I-R(1;\cdot ),\tilde{y})={(-1)}^{\tau }={(-1)}^{{\sigma }_n}=-1.$$

(30)

Likewise, we could obtain

$$index(I-R(1;\cdot ),\tilde{y})={(-1)}^0=1.$$

(31)

Now, for all $t\in [0,1]$, there must be a positive constant $S$, and the positive solution of (29) meets $S^{-1}\le u_1,u_2\le S$. So, for all $t\in [0,1]$, $R(t;y)\ne y$ on $\partial B$. Using the homotopy invariance of the topological degree,

$$\mathrm{deg}(I-R(1;\cdot ),0,B(S))=\mathrm{deg}(I-R(0;\cdot ),0,B(S)).$$

(32)

Meanwhile, according to a previous assumption, the two equations $R(1;\mathbf{u})=0$ and $R(0;\mathbf{u})=0$ only have a positive solution $\tilde{u}$ in $B(S)$, so,

$$\begin{array}{l}\mathrm{deg}(I-R(1;\cdot ),0,B(S))=index(I-R(1;\cdot ),\tilde{y})=-1,\\ \mathrm{deg}(I-R(0;\cdot ),0,B(S))=index(I-R(0;\cdot ),\tilde{y})=1,\end{array}$$

This contradicts the assumption that the proof ends. □

5. Numerical Simulation

In this section, we use numerical simulations to validate our theoretical simulations, including the stability of the coexistence equilibrium points in Section 3, the existence of constant and nonconstant equilibrium solutions in Section 4, that is, system (11) with self-diffusion does not form a space pattern, and in system (5) with cross-diffusion, the system forms a spatial stationary pattern.

Firstly, we use numerical simulations to illustrate the stability of the coexistence equilibrium points in the ODE system. We choose appropriate parameters in ODE system (7), and the two populations tend to be stable and coexist, as shown in Figure 1. For PDE system (5) with cross-diffusion, Figure 2 and Figure 3 show the relation between the real part of the eigenvalue $\lambda$ and cross-diffusion coefficient $c_3$ or Allee effect constant $\sigma$. Due to Theorem 3, we state that when the real part of the maximum eigenvalue is greater than zero, the system equilibrium point will become unstable. Figure 2 shows the real part of the eigenvalues $\lambda$ as a function of $c_3$, and the cross-diffusion coefficient can affect the stability of the coexistence equilibrium points. When $\sigma$ is 0.3, 0.5, or 0.8, the critical value of $c_3$ is 173.4, 208.4, or 277.3, respectively. For different $\sigma$, the critical value of $c_3$ is different. It can be concluded that the Allee effect constant also affects the stability of the coexistence equilibrium points. Figure 3 shows that when the parameter $c_3$ is fixed, there is also a critical value for $\sigma$ to induce the instability of coexistence equilibrium points of the system (5). When $c_3$ is 200, 300, or 500, the critical value of $\sigma$ is 0.455, 0.882, or 1.385, respectively.

Next, the numerical simulations show that the system (11) with self-diffusion does not form a space pattern, and in the PDE system (5) with cross-diffusion, the system will form a spatial stationary pattern. We suppose the area of system (5) is a rectangular area $\Omega =[0,L]\times [0,L]\subset R^2$ with $L=200.$ Additionally, we study the system (5) in a domain with 200 × 200 stations using the simple Euler method with the Laplacian operator in a discrete grid of lattice points denoted by $(p,q), p=1,2\cdots N,q=1,2\cdots N$. The form is

$$\Delta k{|}_{(i,j)}=\frac{1}{s^2}(k(p-1,q)+k(p+1,q))+k(p,q-1)+k(p,q+1)-4\mathrm{k}(p,q) ,$$

(33)

$s$ is the lattice constant. Let $s=1$. When $(p,q)$ is on the left boundary, we define $k(p,0)\equiv k(p,1)$. When $(p,q)$ is on the right boundary, we define $k(p,N+1)\equiv k(p,N)$. When $(p,q)$ is on the upper boundary, we define $k(0,q)\equiv k(1,q)$. When $(p,q)$ is on the lower boundary, we define $k(N+1,q)\equiv k(N,q)$. This indicates that the population flux across the border is zero.

Then, we choose the Allee effect constant that makes system (11) with self-diffusion stable. Figure 4 illustrates changes in the population density of prey and predators in system (11). It can be seen that system (11) does not form a space pattern. This verifies the result of Theorem 4 in Section 4, and it shows that the system with self-diffusion cannot induce a space pattern. Next, we select suitable parameters to generate the spatial patterns. They satisfy condition (8) and (13). Figure 5 shows the spatial pattern for system (5). When the time reaches 30,000, the spatial pattern is unchanged. So, as time increases, it shows that the system forms a spatial stationary pattern. Numerical simulations illustrate our theoretical results of Theorem 5.

6. Concluding Remarks

In our article, we studied a spatial predator–prey model with the Allee effect and nonlinear reaction cross-diffusion to research the Turing instability and pattern formation caused by cross-diffusion. Moreover, the Allee effect produces considerable influence on the stability of the system.

For ODE system (7) and system (11) with self-diffusion, if the Allee effect constant is within the stable range, the coexistence equilibrium of the model is asymptotically stable (Theorem 1 and Theorem 2), which means that two populations will tend to be stable and coexist in ODE system (7) (see Figure 1). For system (5) with nonlinear cross-diffusion, the system will induce Turing instability when the cross-diffusion coefficient $c_3$ is big enough and satisfies conditions (8) and (13) (Theorem 3). Figure 2 shows the relationship between the real part of the maximum eigenvalue $\lambda$ and the cross-diffusion coefficient $c_3$. When $c_3$ is big enough, the real part of the maximum eigenvalue is positive, that is, the system is unstable. Figure 3 shows the relationship between the real part of the maximum eigenvalue $\lambda$ and the Allee effect constant. The critical value of $c_3$ increases with the value of the Allee effect constant. Thus, the Allee effect also affects the stability of system (5).

Furthermore, in Theorem 4, if the Allee effect constant is within the stable range, we prove that the model has a constant solution for system (11) with self-diffusion. Choosing appropriate parameters, numerical simulation shows the spatial patterns. It shows predator and prey densities do not change over time, and finally, it reaches the stationary uniform solution. The system cannot generate spatial patterns (see Figure 4). When the cross-diffusion is present, we prove reaction system (5) could produce nonconstant, positive, steady-state solutions (Theorem 5). Figure 5 shows stationary patterns in system (5). It shows that a predator–prey ecosystem with self-diffusion does not form spatial patterns, while in an ecosystem with cross-diffusion, the population forms spatial patterns.

In our model, the Allee effect and the nonlinear cross-diffusion term were introduced together to study the stability and spatial patterns of the system. Biologically, due to the fact that the impacts of the Allee effect on different predator–prey systems are quite different, under the Allee effect, a system with self-diffusion will change its stability at the equilibrium point. Meanwhile, the mobility of a population is largely influenced by the abundance of the presence or absence of another population. The uniform steady state of a self-diffusion system is stable under nonuniform perturbations, but because of cross-diffusion, it loses its stability and produces various modes. This explains the destabilizing effect of nonlinear cross-diffusion.