Pattern Dynamics in a Spatial Predator-Prey System with Allee Effect

and Applied Analysis 3 It is found from direct calculations that E 3 is unstable and E 4 is stable. Denote E 4 = (u ∗ , V). Following the standard linear analysis of the reactiondiffusion equation [22], we consider a perturbation near the steady state: u ( ⃗ r, t) = u ∗ + u (r, t) , V ( ⃗ r, t) = V ∗ + V (r, t) , (11) where u(r, t) ≪ u, V(r, t) ≪ V, and r = (X, Y). Assume that ( u (r, t) V (r, t) ) = ( α 1 α 2 ) e λt e i(κXX+κYY), (12) where λ is the growth rate of perturbation in time t, α 1 and α 2 represent the amplitudes, and κ X and κ Y are the wave number of the solutions. The characteristic equation of the systems (4a) and (4b) is (A − λI) ( u V ) = 0, (13)


Introduction
The Allee effect, named after the ecologist Warder Clyde Allee, has been recognized as an important phenomenon of positive density dependence in low-density population [1][2][3][4][5].Allee effect can occur whenever fitness of an individual in a small or sparse population decreases as the population size or density also declines [6,7].Since the outstanding work of Allee [1], the Allee effect has been regarded as one of the central and highly important issues in the population and community ecology.And its critical importance has widely been realized in the conservation biology that Allee effect is most likely to increase the extinction risk of low-density populations.As a result, studies on Allee effect have received more and more attention from both mathematicians and ecologists.
Long time series of the density of both prey and predator is needed, so it is difficult to analyse their dynamics.As a result, it may provide useful information by constructing mathematical models to investigate the dynamical behaviors of predator-prey systems.There have been a large group of papers on predator-prey systems with Allee effect [8][9][10][11][12][13].
However, these previous works did not take into account the effect of space.
There are also some works done on spatial predatorprey systems with Allee effect [14][15][16].Petrovskii et al. found that the deterministic system with Allee effect can induce patch invasion [14].Morozov et al. found that the temporal population oscillations can exhibit chaotic dynamics even when the distribution of the species in the space was regular [15].Moreover, they found that the chaos accompanied with patch invasion even though the environments were heterogeneous [16].However, their results were obtained by choosing particular initial conditions.Then, it is natural to ask what kind of patterns can be obtained in predator-prey systems with Allee effect by using other initial conditions.To understand that mechanism well, we will investigate a predator-prey system with Allee effect.
Because of the insightful work of many scientists over recent years, we can make research on pattern selection by using the standard multiple scale analysis [17,18], in which the control parameters and the derivatives are expanded in terms of a small enough parameter.In the neighborhood of the bifurcation points (Hopf and Turing bifurcation points), the critical amplitudes follow the normal forms, and thus their general forms can be obtained from the methods of symmetry-breaking bifurcations.
The paper is organized as follows.In Section 2, we present a predator-prey system with Allee effect and give Turing region in parameters space.In Section 3, by using multiple scale analysis, we obtain amplitude equations.In Section 4, we show the spatial patterns by a series of numerical simulations.Finally, conclusions and discussions are presented in Section 5.
The function () represents the intrinsic prey growth, (, ) = () represents predation term,  is the food utilization coefficient,  1 and  2 are diffusion coefficients, and () describes predator mortality.
It is assumed that the predation term is a bilinear form of prey and predator density and predator mortality is a nonlinear function of predator density.As a result, we choose (, ) =  and () =  2 [20].
When the prey population obeys Allee dynamics, its growth rate can be parameterized as follows [14,15,21]: where  is the prey-carrying capacity,  is the maximum per capita growth rate, and  0 quantifies the intensity of the Allee effect.If 0 <  0 < , () is a strong Allee effect; if − <  0 < 0, () is a weak Allee effect; if  0 ≤ −, the Allee effect is absent.
In order to minimize the number of parameters involved in the model system, it is extremely useful to write the system in a nondimensionalized form.Although there is no unique method of doing this, it is often a good idea to relate the variables to some key relevant parameters.Introducing dimensionless variables we obtain the following equations: where ( First of all, we need to investigate the dynamics of nonspatial model of systems (4a) and (4b) Systems (6a) and (6b) have three boundary equilibrium named  0 = (0, 0),  1 = (1, 0), and  2 = (, 0) and two interior equilibriums named  3 and  4 , where ) , (7a) where  = () 2 − 2() 2  − 2 + 1.
From a biological point of view, we are concerned with the dynamics of  3 and  4 .The Jacobian matrix corresponding to the equilibrium point is that where Diffusion-driven instability requires the stable, homogeneous steady state is driven unstable by the interaction of the dynamics and diffusion of the species; and therefore It is found from direct calculations that  3 is unstable and  4 is stable.Denote  4 = ( * , V * ).
Following the standard linear analysis of the reactiondiffusion equation [22], we consider a perturbation near the steady state: where (, ) ≪  * , V(, ) ≪ V * , and  = (, ).Assume that where  is the growth rate of perturbation in time ,  1 and  2 represent the amplitudes, and   and   are the wave number of the solutions.The characteristic equation of the systems (4a) and ( 4b) is where ) .
As a result, we have characteristic polynomial: = where  2 =  2  +  2  .The roots of ( 15) can be obtained by the following form: When Im(  ) ̸ = 0 and Re(  ) = 0, Hopf bifurcation will emerge.Then, we have that the critical value of Hopf bifurcation parameter- equals When  2 = (  ) 2 = √Δ  / and Im(  ) = 0, Re(  ) = 0, Turing bifurcation will occur.Denote   as the critical value of  as Turing instability occurs.Since the expression is complicated, we omit it here.
In Figure 1, we show the two critical lines in the parameter space spanned by  and .The equilibria that can be found in the region, marked by  (Turing space), are stable with respect to the homogeneous perturbations, but they lose their stability with respect to the perturbations of specific wave numbers .In this region, stationary patterns can be observed.To see the effect of parameter  well, we plot in Figure 2 the dispersion relation corresponding to several values of  while keeping the other parameters fixed.We see that the available Turing modes shift to higher wave numbers when  decreases.

Spatial Dynamics of Systems (4a) and (4b)
In the following, we use multiple scale analysis to determine the amplitude equations when || =   .Denote  as the controlled parameters.When the controlled parameter is larger than the critical value of Turing point, the solutions of systems (4a) and (4b) can be expanded as with || =   .  and the conjugate   are the amplitudes associated with the modes   and −  .Close to onset  =   , one has that Based on the center manifold near the Turing bifurcation point, it can be concluded that amplitude   satisfies In order to obtain the amplitude equations, we first need to investigate the linearized form of systems (4a) and (4b) at the equilibrium point  4 .By setting  =  * + and V = V * +, we have the following equations: Close to onset  =   , the solutions of systems (4a) and (4b) can be expanded as series form: System ( 19) can be expanded as where  0 = (( * 11  +  * 11 )/(2 * 21 ), 1)  is the eigenvector of the linearized operator.
From the standard multiple scale analysis, up to the third order in the perturbations, the spatiotemporal evolution of the amplitudes can be described as Due to spatial translational symmetry, we have the following equation: Comparing ( 25) with (26), one can find that the two equations hold only if   =   + ⋅ ⋅ ⋅ +   .From the center manifold theory, we know that amplitude equations do not include the amplitude with unstable mode.As a result, we have the following equations: where  = (  − )/  and  0 is a typical relaxation time.
In the following part, we will give the expressions of  0 , ℎ,  1 , and  2 .Let Then systems (4a) and (4b) can be written as: where  = ( We need to investigate the dynamical behavior when  is close to   , and thus we expand  as: where  is a small enough parameter.We expand  and  as the series form of : ) . ( Linear operator  can be expanded as where ) .(34) Let and   is a dependent variable.For the derivation of time, we have that The solutions of systems (4a) and (4b) have the following form: This expression implies that the bases of the solutions have nothing to do with time and the amplitude  is a variable that changes slowly.As a result, one has the following equation: Substituting the above equations into (29) and expanding (29) according to different orders of , we can obtain three equations as follows: where We first consider the case of the first order of .Since   is the linear operator of the system close to the onset, ( 1 ,  1 )  is the linear combination of the eigenvectors that corresponds to the eigenvalue zero.Since that we have that where |  | =  *  and   is the amplitude of the mode exp(  ).Now, we consider the case of the second order of .Note that According to the Fredholm solubility condition, the vector function of the right hand of the above equation must be orthogonal with the zero eigenvectors of operator L +  .And the zero eigenvectors of operator L +  are It can be found from the orthogonality condition that where    and    represent the coefficients corresponding to exp(  ) in   and   .
By investigating exp( 1 ⃗ ), one has ) . ( It can be obtained from the orthogonality condition that By using the same methods, one has where Using the Fredholm solubility condition, we can obtain where By using the same methods, we can obtain the other two equations.The amplitude   can be expanded as As a result, we have The other two equations can be obtained through the transformation of the subscript of .By calculations, we obtain the expressions of the coefficients of  0 , ℎ,  1 , and  2 as follows: By using substitutions, we have where  =  1 +  2 +  3 .In order to see the relationships between different parameters, we give the values of coefficients for different parameter sets in Table 1. The dynamical systems (4a) and (4b) possess five kinds of solutions [23] as follows.

Spatial Pattern of Systems (4a) and (4b)
In this section, we perform extensive numerical simulations of the spatially extended systems (4a) and (4b) in twodimensional spaces.All our numerical simulations employ the zero-flux boundary conditions with a system size of 200 × 200.The space step is Δ = 1, and the time step is Δ = 0.00001.
In Figure 3, we show the spatial pattern of prey population at different time.In the parameter set,  = 1.5,  = 0.15, and  = 0.92, we find that  ∈ ( 3 ,  4 ), which means that there is coexistence of spotted and stripe patterns.As shown in this figure, our theoretical results are consistent with the numerical results.
By setting  = 1.5,  = 0.15, and  = 0.96, one can obtain that  >  4 .In Figure 4, we show the spatial pattern of prey population when  equals 0, 50, 100, 200, 500, and 1000.At the initial time, the prey population shows patched invasion.As time increases, stripe pattern appears and the structure does not change a lot.While keeping other parameters fixed and increasing , we find that stripe pattern will occupy the whole space.However, some stripe patterns connect with each other and cause the emergence of spotted patterns which are shown in Figure 5.
Figure 6 shows the evolution of the spatial pattern of prey population at  = 0, 100, 300, 500, 1000, and 2000 iterations, with small random perturbation of the stationary solution of the spatially homogeneous systems (4a) and (4b).The corresponding parameters values are  = 1.5,  = 0.15, and  = 1.04.By the amplitude equations, we can conclude that there are spotted patterns of prey population for this parameter set.In this case, one can see that for the systems (4a) and (4b), the random initial distribution leads to the formation of an irregular transient pattern in the domain.After these forms, it grows slightly and spotted patterns emerge.When the time is large enough, the spotted patterns prevail over the two-dimensional space.As time further increases, the pattern structures of the prey population do not undergo any further changes.

Conclusion and Discussion
Allee effect has been paid much attention due to its strong potential impact on population dynamics [24].In this paper, we investigated the pattern dynamics of a spatial predatorprey systems with Allee effect.Based on the bifurcation analysis, exact Turing pattern region is obtained.By using amplitude equations, the Turing pattern selection of the predator-prey system is well presented.It is found that the predator-prey systems with Allee effect have rich spatial dynamics by performing a series of numerical simulations.
It should be noted that our results were obtained under the assumption that predation is modeled by the bilinear function of the prey and predator densities.However, this function has limitations to describe many realistic phenomena in the biology.By numerical simulations, we find that the system exhibits similar behaviour when the functional response is of other types, such as Holling-II and Holling-III forms.
To compare the spatial dynamics for different parameters, we give the spatial patterns of population  when the parameter values are out of the domain of Turing space.For this parameter set, systems (4a) and (4b) have Hopf bifurcation, and spiral waves occupy the whole domain instead of stationary patterns, which is shown in Figure 7.The stability of spiral wave can be done by using the spectrum theory analysis [25,26].In the further study, we will use the spectrum theory to show the stability of spiral wave.
In [15], they found that a spatial predator-prey model with Allee effect and linear death rate could increase the system's complexity and enhance chaos in population dynamics.However, in this paper, we showed that a spatial population model with Allee effect and nonlinear death rate can induce stationary patterns, which is different from the previous results.
From a biological point of view, our results show that predator mortality plays an important role in the spatial invasions of populations.More specifically, low predator mortality will induce stationary patterns (cf.Figures 3-6), and high predator mortality corresponds to travelling patterns (cf. Figure 7).When the populations exhibit wave distribution in space, the dynamics of populations may be accompanied with chaotic properties [27,28].If the chaotic behavior occurs, it may lead to the extinction of the population, or the population may be out of control [29,30].In that case, we need to find out the best way to control the chaos or change the chaotic behavior.

Figure 1 :
Figure 1: Bifurcation diagram for the systems (4a) and (4b).The green one is the Hopf bifurcation critical line and the red one, Turing bifurcation critical line.The figure shows the Turing space which is marked by .Parameters values:  = 1.5 and  = 0.15.

Table 1 :
Coefficients for different parameter sets.