Spatiotemporal Pattern in a Self- and Cross-Diffusive Predation Model with the Allee Effect

Abstract

This paper proposes and analyzes a mathematical model for a predator-prey interaction with the Allee effect on prey species and with self- and cross-diffusion. The effect of diffusion which can drive the model with zero-flux boundary conditions to Turing instability is investigated. We present numerical evidence of time evolution of patterns controlled by self- and cross-diffusion in the model and find that the model dynamics exhibits a cross-diffusion controlled formation growth to spotted and striped-like coexisting and spotted pattern replication. Moreover, we discuss the effect of cross-diffusivity on the stability of the nontrivial equilibrium of the model, which depends upon the magnitudes of the self- and cross-diffusion coefficients. The obtained results show that cross-diffusion plays an important role in the pattern formation of the predator-prey model. It is also useful to apply the reaction-diffusion model to reveal the spatial predation in the real world.

1. Introduction

It is well known that predator-prey interaction is one of the basic interspecies relations for ecological and social models. The mathematical model of predator-prey type has played a major role in the studies of biological invasion of foreign species, epidemics spreading, and extinction/spread of flame balls in combustion or autocatalytic chemical reaction [1]. The first differential equation model of predator-prey type Lotka-Volterra equation was formulated by Lotka and Volterra in 1920s, when attempts were first made to find ecological laws of nature [2, 3]. In recent years, one of the important predator-prey models is the Holling type II model, originally due to Holling [4, 5], which has been studied in many articles; see, for example, [6–12].

In population dynamics, any mechanism that can lead to the positive relationship between population density and the per capita growth rate can be termed an Allee effect [13–15], starting with the pioneer work of Allee [16]. The Allee effect can be caused by a number of sources such as difficulties in finding mates at small densities, reproductive facilitation, predation, environment conditioning, and inbreeding depression. This effect usually saturates or disappears as populations get larger [17, 18]. The Allee effect’s strong potential impact on population dynamics has been attracting much attention recently. Detailed investigations relating to the Allee effect may be found in [19–33].

On the other hand, in nature, there is a tendency that the preys would keep away from predators and the escape velocity of the preys may be taken as proportional to the dispersive velocity of the predators. In the same manner, there is a tendency that the predators would get closer to the preys and the chase velocity of predators may be considered to be proportional to the dispersive velocity of the preys [34, 35]. In view of these, cross-diffusion arises, which was proposed first by Kerner [36] and first applied in competitive population system by Shigesada et al. [37] to describe the population pressure by other species. Cross-diffusion expresses the population fluxes of one species due to the presence of the other species. In the realistic situation, self-diffusion implies passive diffusion where the diffusing species moves along its concentration gradient, while cross-diffusion implies countertransport. The cross-diffusion coefficient may be positive or negative. When it is positive, one species moves along the concentration gradient of another species, and when it is negative, one species moves against the concentration gradient of another species [38, 39].

Among theoretical and mathematical biologists, there has been considerable interest to investigate the dynamical behavior of a predation system of interacting populations by taking into account the effect of self- and cross-diffusion [40–48]. But in the studies on spatiotemporal dynamics of a Holling type II predator-prey model with an Allee effect on prey, little attention has been paid to study the effect of self- and cross-diffusion. Based on the previous discussions, the main purpose of this paper is to examine out the stability behavior of the Holling-type-II predator-prey model with Allee effect on prey by taking into account self- and cross-diffusion.

The organization of this paper is as follows. In Section 2, we present a Holling-type-II predation model in the presence of Allee effect with self- and cross-diffusion. In Section 3, we analyze the model and derive the condition of Turing instability with respect to parameters. Furthermore, we use numerical simulations to reveal the emergence of different patterns and the influence of cross-diffusion on the dynamical behavior of the model. Finally, a brief discussion is given in Section 4.

2. Model

In this paper, we consider the following model of two-dimensional (i.e., 2𝐷) population dynamics: ̇(𝑁)(),̇()𝑁=𝐹−𝑓𝑁,𝑃𝑃=𝜇𝑓𝑁,𝑃−𝑀𝑃, (1) where 𝑁 and 𝑃 are the densities of prey and predator at time 𝑇>0, respectively. The function 𝐹(𝑁) describes prey multiplication, 𝑓(𝑁,𝑃) describes predation, and the term 𝑀𝑃 stands for predator mortality. 𝜇 is the food utilization coefficient. We assume that the predator response is of Holling-type-II which is usually parameterized as follows: 𝑓(𝑁,𝑃)=𝐴𝑁𝑃,𝑁+𝐵 (2) where 𝐴 is the predation rate and 𝐵 is the half-saturation density.

Referring to [18], we assume that Allee dynamics for the prey population, its growth rate can be parameterized as follows: 𝐹(𝑁)=4𝜔𝐾−𝑁₀²𝑁𝑁−𝑁₀(𝐾−𝑁), (3) where 𝐾 is the prey carrying capacity, 𝜔 is the maximum per capita growth rate, and 𝑁₀ is the “threshold” (𝑁₀<𝐾), so that for 𝑁<𝑁₀, the growth rate becomes negative. The value of 𝑁₀ can be considered as a measure of the intensity of the Allee effect: the less the value of 𝑁₀ is, the less prominent is the Allee effect [21]. In other words, the Allee effect is called “strong” if 0<𝑁₀<𝐾 (when the growth rate becomes negative for 𝑁<𝑁₀) and “weak” if −𝐾<𝑁₀≤0. For 𝑁₀≤−𝐾, the Allee effect is absent [49].

Based on the Holling-type-II predator-prey model with the Allee effect for the prey, in the following, we investigate the dynamical complexity induced by self- and cross-diffusion in the model. By choosing appropriate scales for the variables of (1), (2), and (3), the number of parameters can be lessened. Considering dimensionless variables as in [21] with the following scaling: 𝑁𝑢=𝐾𝑃,𝑣=𝜇𝐾,𝑡=𝜇𝐴𝐾𝑇𝐵, (4) we obtain and employ the corresponding diffusive model as follows: 𝜕𝑢()−𝜕𝑡=𝛾𝑢𝑢−𝛽1−𝑢𝑢𝑣1+𝛼𝑢+𝐷₁₁∇²𝑢+𝐷₁₂∇²𝑣,𝜕𝑣=𝜕𝑡𝑢𝑣1+𝛼𝑢−𝛿𝑣+𝐷₂₁∇²𝑢+𝐷₂₂∇²𝑣. (5) Here, 𝐾𝛼=𝐵𝑁,𝛽=⁰𝐾,𝛾=4𝜔𝐵𝐾𝜇𝐴𝐾−𝑁₀²,𝛿=𝑀𝐵.𝜇𝐴𝐾 (6)

In model (5), the nonnegative constants 𝐷₁₁ and 𝐷₂₂ are the self-diffusion coefficients of 𝑢 and 𝑣, respectively, which imply the movement of individuals from a higher to lower concentration region. The cross-diffusion coefficients 𝐷₁₂ and 𝐷₂₁ express the respective population fluxes of the preys and predators resulting from the presence of the other species, respectively. That is to say, 𝐷₁₂ in (5) represents the tendency of the prey 𝑢 to keep away from its predators 𝑣, and 𝐷₂₁ represents the tendency of the predator to chase its prey. The cross-diffusion coefficients, 𝐷₁₂ and 𝐷₂₁, may be positive or negative. Positive cross-diffusion coefficient denotes that one species tends to move in the direction of lower concentration of another species, and negative cross-diffusion expresses the population fluxes of one species in the direction of higher concentration of another species [39, 50]. Note that ∇²=𝜕²/𝜕𝑥²+𝜕²/𝜕𝑦² is the usual Laplacian operator in a two-dimensional space.

Model (5) is to be analyzed with the following nonzero initial conditions: 𝑢[]×𝑥,𝑦,𝑡>0,𝑣𝑥,𝑦,𝑡>0,𝑥,𝑦∈Ω=0,𝐿𝑥0,𝐿𝑦 (7) and the following zero-flux boundary conditions: 𝜕𝑢=𝜕𝐧𝜕𝑣𝜕𝐧=0,𝑥,𝑦∈𝜕Ω. (8)

In the previous part, 𝐿𝑥 and 𝐿𝑦 give the size of the model in the directions of 𝑥 and 𝑦, respectively. 𝐧 is the outward unit normal vector of the boundary 𝜕Ω which we will assume to be smooth. Zero-flux boundary conditions imply that no external input is imposed from outside [1, 51].

For model (5), without diffusion terms, there are four stationary states: 𝐸₀=(0,0), 𝐸₁=(1,0), 𝐸₂=(𝛽,0), and the coexistence state 𝐸^∗=(𝑢^∗,𝑣^∗), where 𝑢^∗=𝛿1−𝛼𝛿,𝑣^∗=𝛾𝑢^∗−𝛽1−𝑢^∗1+𝛼𝑢^∗. (9)

Using the standard linear stability analysis, it can be readily seen that(i) 𝐸₀=(0,0) is always a stable node and corresponds to extinction of both prey and predator;(ii) two uniform steady states 𝐸₁=(1,0) and 𝐸₂=(𝛽,0) are either nodes or saddle points, depending on the relation between the parameters, and correspond to extinction of the predator;(iii) the coexistence state 𝐸^∗=(𝑢^∗,𝑣^∗) lies in the biologically meaningful domain 𝑢>0, 𝑣>0 under condition 𝛽<𝑢^∗<1 (for these parameter values 𝐸₁=(1,0) and 𝐸₂=(𝛽,0) are saddle points) and it can be either a stable or unstable focus or a node.  For more details about model (5) without diffusions, we can refer to reference [29].

When we take the parameters as 𝛼=1.11, 𝛽=0.1, 𝛾=11.25, and 𝛿=0.36 for (5), 𝐸₀=(0,0) is a stable node, 𝐸₁=(1,0) and 𝐸₂=(0.1,0) are saddle points, and the positive equilibrium 𝐸^∗=(0.6,3.75) exists, which is unstable and shown in Figure 1.

Figure 1 

Phase portrait of model (5). The parameters are taken as 𝛼=1.11, 𝛽=0.1, 𝛾=11.25, and 𝛿=0.36. The horizontal axis is the prey population 𝑢, and the vertical axis is the predator population 𝑣. 𝐸0=(0,0) is a stable node, 𝐸1=(1,0) and 𝐸2=(0.1,0) are saddle points, and 𝐸∗=(0.6,3.75) is unstable.

3. Pattern Formation

To study patterns of the self- and cross-diffusion model given by (5), we must consider a spatially homogeneous system.

In this section, we perform extensive numerical simulations of the spatially extended model (5) in a two-dimensional space, and the qualitative results are shown here. All our numerical simulations employ the zero-flux boundary conditions (8) with a system size of 𝐿𝑥×𝐿𝑦 discretized through 𝑥→(𝑥₀,𝑥₁,𝑥₂,…,𝑥_𝑛) and 𝑦→(𝑦₀,𝑦₁,𝑦₂,…,𝑦_𝑛) with 𝑛=200, and the time step of 𝑡=1/100.

We use the standard five-point approximation for the 2𝐷 Laplacian with the zero-flux boundary conditions [52, 53]. That is, the concentrations (𝑢^𝑛+1_𝑖,𝑗,𝑣^𝑛+1_𝑖,𝑗) at the time (𝑛+1)𝑡 at the mesh position (𝑥_𝑖,𝑦_𝑗) are given by 𝑢^𝑛+1_𝑖,𝑗=𝑢^𝑛_𝑖,𝑗+𝑡𝐷₁₁Δ_ℎ𝑢^𝑛_𝑖,𝑗+𝑡𝐷₁₂Δ_ℎ𝑣^𝑛_𝑖,𝑗+𝑡𝑓𝑢^𝑛_𝑖,𝑗,𝑣^𝑛_𝑖,𝑗,𝑣^𝑛+1_𝑖,𝑗=𝑣^𝑛_𝑖,𝑗+𝑡𝐷₂₁Δ_ℎ𝑢^𝑛_𝑖,𝑗+𝑡𝐷₂₂Δ_ℎ𝑣^𝑛_𝑖,𝑗+𝑡𝑔𝑢^𝑛_𝑖,𝑗,𝑣^𝑛_𝑖,𝑗, (10) with the Laplacian defined by Δ_ℎ𝑢^𝑛_𝑖,𝑗=𝑢^𝑛_𝑖+1,𝑗+𝑢^𝑛_{𝑖−1,𝑗}+𝑢^𝑛_𝑖,𝑗+1+𝑢^𝑛_{𝑖,𝑗−1}−4𝑢^𝑛_𝑖,𝑗ℎ², (11) where the space stepsize ℎ=1/3.

We are mainly interested in investigating the behavior of model (5) around the interior equilibrium point; so we will put emphasis on the positive equilibrium point 𝐸^∗=(𝑢^∗,𝑣^∗). The entire system is initially placed in the stationary state (𝑢^∗,𝑣^∗), and the propagation velocity of the initial perturbation is thus of the order of 5×10⁻⁴ space units per time unit. After the initial period during which the perturbation spreads, either the system goes into a time-dependent state or goes into a time-independent steady state [54].

To study the effect of diffusion on model (5), we consider the linearized form of the model about 𝐸^∗=(𝑢^∗,𝑣^∗) as follows: 𝜕𝑢𝜕𝑡=𝐽₁₁𝑢+𝐽₁₂𝑣̂+𝐷₁₁∇²𝑢+𝐷₁₂∇²̂𝜕𝑣̂𝑣,𝜕𝑡=𝐽₂₁𝑢+𝐽₂₂𝑣̂+𝐷₂₁∇²𝑢+𝐷₂₂∇²𝑣,̂ (12) where 𝑢=𝑢^∗+𝑢, 𝑣=𝑣^∗+𝑣̂ and 𝐽₁₁=−𝑑(𝑏−𝑎−𝑎𝑏)𝑛²+𝑎²𝑑()+()𝑑2𝑛+𝑑𝑎𝑏+𝑎+𝑏²×()−2𝑏𝑑𝑛𝑏𝑛𝑛−𝑑²⁻¹,𝐽₁₂=−𝑑𝑚𝑛,𝐽₂₁()=−𝑎𝑑+𝑑−𝑛)(𝑎𝑑+𝑏𝑑−𝑏𝑛(),𝐽𝑏𝑚𝑛−𝑑₂₂=0. (13) We note that (𝑢,𝑣)̂ are small perturbations of (𝑢,𝑣) about the equilibrium point (𝑢^∗,𝑣^∗) and the stationary state can become unstable to a nonuniform perturbation 𝑢(),𝑣̂().𝐫,𝑡∼exp𝜆𝑡+i𝑘⋅𝐫𝐫,𝑡∼exp𝜆𝑡+i𝑘⋅𝐫,𝐫=𝑥,𝑦 (14) Any solution of model (12) can be expanded into a Fourier series so that 𝑢(𝐫,𝑡)=^∞_𝑖,𝑗=0𝑢_𝑖𝑗(𝐫,𝑡)=^∞_𝑖,𝑗=0𝛼_𝑖𝑗(𝑡),𝑣̂()=sin𝑘𝐫𝐫,𝑡^∞_𝑖,𝑗=0𝑣_𝑖𝑗(𝐫,𝑡)=^∞_𝑖,𝑗=0𝛽_𝑖𝑗(𝑡)sin𝑘𝐫, (15) where 0<𝑥<𝐿𝑥, 0<𝑦<𝐿𝑦. Note that 𝑘=(𝑘_𝑖,𝑘_𝑗), 𝑘_𝑖=𝑖𝜋/𝐿𝑥, and 𝑘_𝑗=𝑗𝜋/𝐿𝑦 are the corresponding wavenumbers.

Substituting 𝑢_𝑖𝑗 and 𝑣_𝑖𝑗 into (12), we obtain d𝛼_𝑖𝑗=𝐽d𝑡¹¹−𝐷₁₁𝑘²𝛼_𝑖𝑗+𝐽₁₂−𝐷₁₂𝑘²𝛽_𝑖𝑗,d𝛽_𝑖𝑗=𝐽d𝑡²¹−𝐷₂₁𝑘²𝛼_𝑖𝑗−𝐷₂₂𝑘²𝛽_𝑖𝑗, (16) where 𝑘²=𝑘²_𝑖+𝑘²_𝑗. The characteristic equation of the linearized system is given by 𝜆²−tr(𝒥)𝜆+det(𝒥)=0, (17) where 𝐽𝒥=¹¹−𝐷₁₁𝑘²𝐽₁₂−𝐷₁₂𝑘²𝐽₂₁−𝐷₂₁𝑘²−𝐷₂₂𝑘²,tr(𝒥)=𝐽₁₁−𝐷₁₁+𝐷₂₂𝑘²,det(𝒥)=det(𝐷)𝑘⁴+𝐷₁₂𝐽₂₁+𝐷₂₁𝐽₁₂−𝐷₂₂𝐽₁₁𝑘²−𝐽₁₂𝐽₂₁. (18) Here, the diffusive matrix is 𝐷𝐷=¹¹𝐷₁₂𝐷₂₁𝐷₂₂, (19) and the determinant of 𝐷 is det(𝐷)=𝐷₁₁𝐷₂₂−𝐷₁₂𝐷₂₁. Then the roots of (17) yield the dispersion relation 𝜆_1,2=12tr(𝒥)±(tr(𝒥))²−4det(𝒥). (20)

In order to have diffusion-induced instability (i.e., Turing instability), which means that it is an asymptotically stable equilibrium of the diffusion-less model (5) but is unstable with respect to solutions of the spatial model (5), at least one of the conditions tr(𝒥)<0 and det(𝒥)>0 must be violated by Routh-Hurwitz criterion. It is evident that the first condition tr(𝒥)<0 is not violated when the condition 𝐽₁₁<0 is met. Hence, in order to obtain diffusion instability, we need det(𝒥)<0. Thus, the condition is given by 𝐺𝑘²=det(𝐷)𝑘⁴+𝐷₁₂𝐽₂₁+𝐷₂₁𝐽₁₂−𝐷₂₂𝐽₁₁𝑘²−𝐽₁₂𝐽₂₁<0. (21) If the minimum of 𝐺(𝑘²) is negative, then a range of modes 𝑘²₋<𝑘²<𝑘²₊ will grow, where 𝑘²₋=1(𝐷)×𝐷2det₂₂𝐽₁₁−𝐷₂₁𝐽₁₂−𝐷₁₂𝐽₂₁−𝐷²₂₂𝐽²₁₁+𝐷₂₁𝐽₁₂+𝐷₁₂𝐽₂₁²+4𝐽₁₂𝐽₂₁det(𝐷)−2𝐷₂₂𝐽₁₁𝐷₂₁𝐽₁₂+𝐷₁₂𝐽₂₁^1/2,𝑘²₊=1(𝐷)×𝐷2det₂₂𝐽₁₁−𝐷₂₁𝐽₁₂−𝐷₁₂𝐽₂₁+𝐷²₂₂𝐽²₁₁+𝐷₂₁𝐽₁₂+𝐷₁₂𝐽₂₁²+4𝐽₁₂𝐽₂₁det(𝐷)−2𝐷₂₂𝐽₁₁𝐷₂₁𝐽₁₂+𝐷₁₂𝐽₂₁^1/2. (22)

In addition, when 𝑘²=𝑘²_𝑐=𝐷₂₂𝐽₁₁−𝐷₁₂𝐽₂₁−𝐷₂₁𝐽₁₂2det(𝐷)>0, (23) we obtain 𝐺_min=−𝐽₁₂𝐽₂₁−𝐷₂₂𝐽₁₁−𝐷₁₂𝐽₂₁−𝐷₂₁𝐽₁₂²(𝐷).4det (24) Then we have 𝐺𝑘²<0⟺𝑘²₋<𝑘²<𝑘²₊. (25) If both 𝑘²₋ and 𝑘²₊ exist and have positive values, they limit the range of instability for a locally stable equilibrium (see Figure 2).

Figure 2 

Model (5) with wavenumbers lying between the zeros of 𝐺(𝑘2), 𝑘2−, and 𝑘2+ growing in the Turing instability. Parameters are taken as 𝛼=1.11, 𝛽=0.1, 𝛾=11.25, 𝛿=0.36, 𝐷11=0.2, 𝐷22=2, 𝐷12=−0.5, and 𝐷21=−0.4; hence, 𝑘2−=0.427433403, 𝑘2+=2.807586678, and 𝐺min=−0.283261256.

Summarizing the above discussions, we can obtain the following theorem.

Theorem 1. The criterion for Turing instability for model (5) satisfies the following condition: 𝐷₂₂𝐽₁₁−𝐷₁₂𝐽₂₁−𝐷₂₁𝐽₁₂>2−𝐽₁₂𝐽₂₁det(𝐷)>0. (26)

In Figure 3, we plot the dispersion relation corresponding to several values of parameter 𝐷₁₂ while keeping the others fixed as 𝐷𝛼=1.11,𝛽=0.1,𝛾=11.25,𝛿=0.36,₁₁=0.2,𝐷₂₂=2,𝐷₂₁=−0.4. (27)

Figure 3 

Variation of dispersion relation of model (5) for different 𝐷12: (A) 𝐷12=−1.6<𝐷T12, the Turing instability occurs; (B) 𝐷12=𝐷T12=−0.31, the critical Turing value; and (C) 𝐷12=0.4>𝐷T12, the Turing instability decays. Other parameters values are fixed as (27) in the text.

Curve (B) in Figure 3 corresponds to the critical Turing value 𝐷^T₁₂=−0.31. When the parameter 𝐷₁₂=−1.6 below the 𝐷^T₁₂ (i.e., curve (A) in Figure 3), the Turing instability occurs, while 𝐷₁₂=0.4 above the 𝐷^T₁₂ (i.e., curve (C) in Figure 3), the Turing instability decays. That is to say, when 𝐷₁₂>𝐷^T₁₂, the steady state is the only stable solution of model (5). The spatial patterns are generated when 𝐷₁₂ passes through the critical Turing bifurcation point 𝐷^T₁₂.

Figure 4 shows the spatial patterns of prey 𝑢 and predator 𝑣 and the corresponding phase portraits of model (5). Figures 4(a1) and 4(a2) show stationary spotted and striped-like patterns emerge mixed in the distribution of prey 𝑢 and predator 𝑣 population density, respectively. In Figures 4(a1) and 4(a2), when 𝐷₁₂=−0.365<𝐷^T₁₂=−0.31, starting with a homogeneous state (0.6,3.75), the random perturbations make the formations end with spotted and striped-like coexistence in the spatial domain at 𝑡=10000 (i.e., 1000000 iterations), and the mixed patterns prevail over the whole domain at last. In Figure 4(a3), one can see that as 𝑡→∞, 𝑢 and 𝑣 coordinates fix to (0.745,3.428) which is the steady state.

Figure 4 

Dynamical behavior of model (5). 𝐷12=−0.365 (the first row); 𝐷12=−0.5 (the second row); and other parameters are fixed as (27) in the text. (a1), (b1) Two categories of Turing patterns of prey 𝑢 at 𝑡=10000; (a2), (b2) two categories of Turing patterns of predator 𝑣 at 𝑡=10000; and (a3), (b3) phase portraits.

Figures 4(b1) and 4(b2) show spotted pattern formations of prey 𝑢 and predator 𝑣, respectively. In Figures 4(b1) and 4(b2), 𝐷₁₂=−0.5<𝐷^T₁₂=−0.31, starting with a homogeneous state (0.6,3.75), the random perturbations make the formations end with a time-independent spotted-patterns at 𝑡=10000 (i.e., 1000000 iterations). From Figures 4(b1) and 4(b2), we can see that the spotted patterns prevail over the whole domain and the dynamics of model (5) do not undergo any further changes. In Figure 4(b3), which is a phase portrait, the steady state (0.725,3.061) exists and is asymptotically stable. From the figure, we can see that (𝑢,𝑣)→(0.725,3.061) with time.

In Figure 4, two different snapshots during the temporal evolution of model (5) are presented in a two-dimensional space under different 𝐷₁₂ values. These results indicate that the effect of cross-diffusion 𝐷₁₂ for pattern formation is tremendous, and model (5) is conditionally stable at unique positive equilibrium.

In order to investigate quantitatively the evolution of model (5), we depict the typical emergence of Turing pattern in the spacetime plots of Figure 5. The method of spacetime plots is to let 𝑦 be a constant, choose the line 𝑦=100 from each pattern snapshot, and pile these lines in time order [55–59]. Figures 5(a) and 5(b) show the evolution process of prey 𝑢 throughout time 𝑡 and space 𝑥 with different 𝐷₁₂, respectively, where time increases from bottom to top, and the horizontal axis represents the spatial location. From Figures 5(a) and 5(b), the single parameter of model (5) namely 𝐷₁₂, which is the diffusion coefficient of the prey, can lead to dramatic changes in the qualitative dynamics of solutions. From a biological perspective, the cross term plays a constructive role in the pattern formation of the prey, that is, changes the Turing pattern into regular spatial pattern with time.

(a)

(b)

Figure 5 

Spacetime plots of prey 𝑢 of model (5). The time interval shown is 1500. (a) 𝐷12=−0.365 and (b) 𝐷12=−0.5. Other parameters are fixed as (27) in the text.

According to the works of many researchers [53, 60–62], the choice of initial conditions can affect the spatiotemporal dynamics of a reaction-diffusion model in ecosystems. In Figure 6, when 𝐷₁₂=−0.5 and other parameters are fixed as (27), we choose the special initial condition of model (5) to investigate the evolution of the spatial pattern of prey 𝑢 at a different evolution time. The initial condition is introduced as the following form: 𝑢𝑥,𝑦,0=𝑢^∗,𝑣𝑥,𝑦,0=𝑣^∗,90<𝑥,𝑦<110, (28) which is a square. Initially, the entire system is placed in the stationary state (𝑢^∗,𝑣^∗)=(0.6,3.75), and the propagation velocity of the initial random perturbations is thus of the order of 5×10⁻⁴ space units per time unit. Also the system is then integrated for 10⁵ time steps, and some images are saved.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 6 

The process of pattern formation of prey 𝑢 in model (5) with 𝐷12=−0.5; other parameters are fixed as (27) in the text. Time steps: (a) 𝑡=10; (b) 𝑡=105; (c) 𝑡=151; (d) 𝑡=179; (e) 𝑡=300; and (f) 𝑡=1000.

From Figure 6(a), we can see that a target wave pattern emerges after perturbation of the stationary solutions 𝑢^∗ and 𝑣^∗ of model (5). Some iterations later, it breaks towards the core, and a spot pattern (interior) with target wave pattern (exterior) occurs (see Figure 6(b)). Then one can see that after the symmetrical target-spotted mixture patterns (see Figures 6(c) and 6(d)), it grows slightly and the spots increase with time (see Figures 6(e) and 6(f)). From Figure 6, we observe that when the iteration time is large enough, the time-independent striped-like and target-spotted mixed patterns take a long time to settle down and stay stable. It finally prevails over the whole domain at 𝑡=1000.

Figure 6 illustrates the formation of the spatial patterns behavior of prey 𝑢 in two dimensions with self- and cross-diffusion, where the initial condition (at 𝑡=0) adds random and nonuniform small perturbations to the equilibrium values.

4. Conclusions and Remarks

In this paper, we investigate the pattern formation of a Holling-type-II predator-prey interacting model with an Allee effect in prey population and with self- and cross-diffusion under the zero-flux boundary conditions. We show that cross-diffusion in both predator and prey species tends to stabilize the predation model. That is to say, the unstable equilibrium (𝑢^∗,𝑣^∗) of model (5) becomes stable under certain conditions. Moreover, the numerical results show that model (5) dynamics exhibits a cross-diffusion controlled formation growth not only to spotted-striped, but also to spotted pattern replication.

By a series of analysis and numerical simulations, we find that the evolutionary process depends on several parameters. In this paper, we fix the values of 𝛼, 𝛽, 𝛾, 𝛿, 𝐷₁₁, 𝐷₂₂, and 𝐷₂₁ as (27) and perform a large number of computer simulations corresponding to different 𝐷₁₂. The results show that model (5) has rich spatiotemporal patterns and different patterns appear on account of the change of cross-diffusion coefficient 𝐷₁₂. Biological speaking, cross-diffusion may have effect on the distribution of the species. It is also shown that if the equilibrium state (𝑢^∗,𝑣^∗) of model (5) with no diffusion is unstable, then the corresponding uniform steady state of the model with self- and cross-diffusion can be made stable by changing cross-diffusion coefficients. Therefore, we hope that the results presented here will be useful for studying the dynamical complexity of ecosystems or physical systems.

Acknowledgment

The author is grateful to the reviewers and the editor of this paper for their helpful comments and suggestions.