Abstract

We investigate a diffusive Leslie-Gower predator-prey model with the additive Allee effect on prey subject to the zero-flux boundary conditions. Some results of solutions to this model and its corresponding steady-state problem are shown. More precisely, we give the stability of the positive constant steady-state solution, the refined a priori estimates of positive solution, and the nonexistence and existence of the positive nonconstant solutions. We carry out the analytical study for two-dimensional system in detail and find out the certain conditions for Turing instability. Furthermore, we perform numerical simulations and show that the model exhibits a transition from stripe-spot mixtures growth to isolated spots and also to stripes. These results show that the impact of the Allee effect essentially increases the model spatiotemporal complexity.

1. Introduction

The dynamics of a predator-prey model in a homogeneous environment can be described by the following reaction-diffusion equations: where and are the densities of prey and predator at time and position , respectively. The Laplace operator describes the spatial dispersal with passive diffusion; and are the diffusion coefficients corresponding to species and . describes the per-capita growth rate of the prey; is the functional response of the predator, which corresponds to the saturation of their appetites and reproductive capacity; , the so-called numerical response, is the per-capita growth rate of the predator [1ā€“4].

Functions , , and can be formulated in various specific situations. In general, is of the standard logistic growth: which was first created by Verhulst [5]. Here is the prey carrying capacity and is the intrinsic growth rate of prey.

Some conventional functional response functions include Holling types I, II, and III (see [6ā€“10]). Among many possible choices of , the Holling type-II functional response is most commonly used in the ecological literature, which is defined by [11]: where describes the maximum predation rate and measures the extent to which environment provides protection to prey . The Leslie-Gower type numerical response is given by which was first proposed by Leslie [12], and has been discussed by Leslie and Gower [13] and Pielou [14]. A modified version of Leslie-Gower functional response is given by Aziz-Alaoui et al. [15, 16]. Here, describes the growth rate of the predator ; has a similar meaning to ; takes on the role of the prey-dependent carrying capacity for the predator; is the extent to which environment provides protection to predator . Hence, we can rewrite model (1) as follows: The biological significance of all parameters in model (6) is as above.

For model (6), in the case of , Du et al. [17, 18] mainly focused attention on the steady-state problem and observed some quite interesting phenomena of pattern formation. In the case of , the so-called Holling-Tanner model, Peng and Wang [19, 20] analyzed the global stability of the unique positive constant steady-state and established the results for the existence and nonexistence of positive nonconstant steadystates; Shi and coworkers [21] studied the global attractor and persistence property, local and global asymptotic stability of the unique positive constant equilibrium, and the existence and nonexistence of nonconstant positive steady-states; Li et al. [22] considered the Turing and Hopf bifurcations of the equilibrium solutions; Liu and Xue [23] found the model exhibits the spotted, black-eye, and labyrinthine patterns. For model (6), that is, , Camara and Aziz-Alaoui [24ā€“26] paid more attention to pattern formation in the spatial domain and observed the labyrinth, chaos, and spiral wave patterns.

On the other hand, in population dynamics, any mechanism that can lead to a positive relationship between a component of individual fitness and either the number or density of conspecific can be termed an Allee effect [27ā€“30], starting with the pioneer work of ecologist Allee [31]. In particular, theoretical investigations have shown that an Allee effect can greatly increase the likelihood of local and global extinction [32] and can lead to a rich variety of dynamical effects. As a consequence, it is necessary to explore the influence of Allee effect in the growth of a population.

The Allee effect has been modeled in different ways [33ā€“37]. From an ecological viewpoint, the Allee effect has been modeled into strong and weak ones [33, 38ā€“42]. In a recent analytic approach by Wang and Kot [38], the Allee effect is ā€œstrongā€ if the sign of the growth function in the limit of law density is negative; that is, It is ā€œweakā€ if the sign of the growth function in the limit of law density is positive; that is, The strong Allee effect introduces a population threshold, and the population must surpass this threshold to grow. In contrast, the weak case has not any threshold [10, 35, 38, 42].

In particular, the growth function considering Allee effect is expressed by the equation: having an additive Allee effect, which was first deduced in [43] and applied in [34ā€“36]. Where is the term of additive Allee effect, and are the Allee effect constants. It should be noted that, if , the Allee effect in (9) is weak, while if , the Allee effect in (9) is strong.

Based on the above discussions, in this paper, we rigorously consider the spatiotemporal dynamics of the following modified Leslie-Gower predation model with the additive Allee effect on prey:

We make a change of variables: And for the sake of convenience, we still use variables , instead of , . Thus, the model to be studied is as follows: where Here, is the growth rate of the predator . is the term of additive Allee effect, and and are the Allee effect constants. is a bounded domain with smooth boundary , and is the outward unit normal vector on . The initial data and are continuous functions on , and the zero-flux boundary conditions mean that model (12) is self-contained and has no population flux across the boundary [44, 45].

By the standard theory for semilinear parabolic systems (see, [46]), we have model (12) that admits a unique classical solution for all time.

The stationary problem of model (12), which may display the dynamical behavior of solutions to model (12) as time goes to infinity, satisfies the following elliptic system: Unless otherwise specified, in this paper, we always assume that ; that is, we only focus on the case of weak Allee effect.

The rest of the paper is organized as follows. In Section 2, we investigate the stability of nonnegative constant steady-state solutions. In Section 3, we mainly give a priori upper and lower bounds for positive solutions of model (14). In Section 4 we discuss existence and nonexistence of nonconstant positive solutions, which might give us some suggestions on the conditions under which the patterns may or may not occur. In Section 5 we first use the method of linearized stability analysis to deduce the conditions under which the Turing instability might occur, and next we perform a series of numerical simulations to show the occurrence of different patterns. Finally, in the last section we make a summary to our results and give some concluding remarks.

2. Dynamics Analysis of Model (12)

2.1. The Existence of the Constant Steady-State Solution

It is easy to verify that model (12) has the following nonnegative constant steady-state solutions:(i)the trivial constant solution (extinction of two species);(ii)the semitrivial constant solution (extinction of the prey);(iii)the semitrivial constant solution (extinction of the predator);(iv)the unique positive constant solution (coexistence of two species), where , and is a real positive root of the cubic where .

By the transformation , (15) is reduced to where

It is worthy to note that if holds, (16) has a real positive root. Considering (17) and (18), from , one can determine , where

Hence, we have the following lemma regarding the existence of the positive constant steady-state solution of model (12).

Lemma 1. If either of the following inequalities holds: model (12) has a unique positive constant steady-state solution .

2.2. Stability of the Constant Steady-State Solution

In this subsection, we will analyze the asymptotical stability of the nonnegative constant solutions for the corresponding reaction-diffusion dynamics (12).

For sake of simplicity, we rewrite model (12) in the vectorial form: where , , and

Let be the eigenvalues of the operator- on with the zero-flux boundary conditions. And set being an orthonormal basis of , and ; then where .

Let be any arbitrary constant steady-state solution of model (12). And the linearization of model (12) at the constant steady-state solution can be expressed by where and

From [46], it is known that if all the eigenvalues of the operator have negative real parts, then is asymptotically stable; if there is an eigenvalue with positive real part, then is unstable; if all the eigenvalues have nonpositive real parts while some eigenvalues have zero real parts, then the stability of cannot be determined by the linearization [10].

For each , is invariant under the operator , and is an eigenvalue of if and only if is an eigenvalue of the matrix for some .

In the following, we denote evaluated at , and . So, the local stability of the constant steady-state solution can be analyzed as follows.

Theorem 2. For any positive parameters,(a) the trivial constant solution is unstable;(b) the semitrivial constant solution is(b1) locally asymptotically stable if or holds(b2) unstable if holds;(c) the semitrivial constant solution is unstable.

Proof. The stability of the constant steady-state solution is reduced to consider the characteristic equation: with (a), for ; the eigenvalues are and , so is unstable.(b)We can obtain .(b1)If or hold, then , so for , ā€‰Hence, is locally asymptotically stable.(b2)When , then . For , , which implies that (27) has at least one root with positive real part. Hence, is an unstable steady-state solution of model (12).(c)Since , for , one of the eigenvalues is , so is unstable.
The proof is complete.

Straightforward calculations show that where .

The determinant of is given bythen, the sign of depends on the factor : where is the same definition as (15), and

Therefore, we have the following.

Theorem 3. Assume that , , and the first eigenvalue subject to the zero-flux boundary conditions satisfies Then the positive constant steady-state solution of model (15) is uniformly asymptotically stable.

Proof. When , , then and . So, for , Note that for any , we have . Therefore, the eigenvalues of the matrix have negative real parts. It thus follows from the Routh-Hurwitz criterion that, for each , the two roots and of all have negative real parts.
In the following, we prove that there exists such that
Let ; then Since as , it follows that
By the Routh-Hurwitz criterion, it follows that the two roots , of all have negative real parts. Thus, there exists a positive constant , such that . By continuity, we see that there exists such that the two roots , of satisfy , . In turn, , . Let Then and (36) holds for .
Consequently, the spectrum of , which consists of eigenvalues, lies in . In the sense of [46], we obtain that the positive constant steady-state solution of model (12) is uniformly asymptotically stable. This ends the proof.

3. A Priori Estimates

In this section, we give a priori estimates for the steady-state solutions of model (14). To prove that we recall the following maximum principle [47].

Lemma 4 (see [47, maximum principle]). Let be a bounded Lipschitz domain in and .(a)Assume that and satisfies If , then .(b)Assume that and satisfies If , then .

Theorem 5. Suppose that . Any positive solution of model (14) satisfies where

Proof. Assume that is a positive solution of (14). Set Applying Lemma 4 to model (14), we obtain that By virtue of the definitions of (), it follows from (45) that and , and So, we have If or and hold, from (48), we get that .
If or and hold, from (47), we get that By simple computations, indicates that . So, if holds, we can obtain and . The proof is complete.

In order to obtain the desired bounds, we need to use the following Harnack inequality due to [48].

Lemma 6 (see [48, Harnack inequality]). Let be a positive solution to , where , satisfying zero-flux boundary conditions. Then there exists a positive constant , such that

Theorem 7. Let be an arbitrary fixed positive number. Then there exist positive , such that if and , any positive solution of model (14) satisfies

Proof. By Theorem 5, we note that and . And it suffices to verify the lower bounds of . We will verify the conclusion by a contradiction argument.
On the contrary, suppose that the conclusion is not true; then there exist sequences and with and the positive solution of model (14) corresponding to , such that and satisfies We observe that Lemma 4 guarantees for all by use of the second equation of (53). On the other hand, by the Harnack inequality (Lemma 6), we know that there is a positive constant independent of , such that for all . Consequently, which contradicts Theorem 5. We finish the proof.

4. Nonexistence and Existence of the Nonconstant Solutions

In this section, we apply the energy method and the topological degree theory [49] to establish some results on the nonexistence and existence of the positive nonconstant solutions of model (14), respectively.

4.1. Nonexistence of the Nonconstant Solutions
4.1.1. Nonexistence of Nonconstant Semitrivial Steady-State Solution

In the case that model (14) has the semitrivial steady-state solutions and , such and should satisfy, respectively,

Note that is the smallest positive eigenvalue of the operator- in subject to the zero-flux boundary conditions. Now, using the energy estimates, for the semitrivial solution, we can claim the following.

Theorem 8. (i) If , the only nonnegative solutions of model (56) are and .
(ii) If , the only nonnegative solution of model (57) is .

Proof. Let be a nonnegative solution of models (56) and (57), respectively. Denote and . Then .
Multiplying the equation in (56) by , we get Then with the PoincarƩ inequality we find that which implies that unless , and must be a constant solution.
In a similar manner, we multiply the equation in (57) by to have In view of , we have and must be a constant solution. This ends the proof.

4.1.2. Nonexistence of the Positive Nonconstant Solutions

This subsection is devoted to the consideration of the nonexistence for the positive nonconstant solutions of model (12), and in the below results, the diffusion coefficients do play a significant role.

Theorem 9. There exists a positive constant such that model (14) has no positive nonconstant solution provided that .

Proof. Let be any positive solution of model (14). Then, multiplying the first equation of model (14) by , integrating over and using Theorem 7, we have that In a similar manner, we multiply the second equation in model (14) by to have Using Theorem 7 again, we have that Similarly, From the well-known PoincarƩ inequality (62), (63), (64), and (65), we obtain that where , .
This shows that if then , which asserts our results.

4.2. Existence of the Nonconstant Positive Solutions

In this subsection, we will discuss the existence of the positive nonconstant solution of model (14).

Unless otherwise specified, in this subsection, we always require that Lemma 1 holds, which guarantees that model (14) has the unique positive constant solution . From now on, we denote .

Let be the space defined in (23) and let We write model (14) in the form: where Then is a positive solution of model (69) if and only if ā€‰ satisfies where is the inverse operator of subject to the zero-flux boundary condition. Then where

If is invertible, by Theorem of [50] the index of at is given by where is the multiplicity of negative eigenvalues of .

On the other hand, using the decomposition (24), we have that is an invariant space under and is an eigenvalue of in , if and only if is an eigenvalue of . Therefore, is invertible, if and only if for any the matrix is invertible.

Let be the multiplicity of . For the sake of convenience, we denote Then, if is invertible for any , with the same arguments as in [51], we can assert the following conclusion.

Lemma 10. Assume that, for all , the matrix is nonsingular; then

To compute , we have to consider the sign of . A straightforward computation yields where .

If , then has two positive solutions given by

Theorem 11. Assume that . If and for some , and is odd, then model (14) has at least one nonconstant solution.

Proof. By Theorem 9, we can fixed and such that(i)model (12) with diffusion coefficients and has no nonconstant solutions;(ii) for all .
By virtue of Theorem 7, there exists a positive constant such that .
Set and define by where It is clear that solving model (14) is equivalent to finding a fixed point of in . Further, by virtue of the definition of , we have that has no fixed point in for all .
Since is compact, the Leray-Schauder topological degree is well defined. From the invariance of Leray-Schauder degree at the homotopy, we deduce
Clearly, . Thus, if model (14) has no other solutions except the constant one , then Lemma 10 shows that
On the contrary, by the choice of and , we have that is the only solution of . Furthermore, by (ii) above, we have
From (83)ā€“(85), we get a contradiction, and the proof is completed.

Corollary 12. Let be fixed. If and all the eigenvalues have odd multiplicity, then, there exists a sequence of intervals with (as ) such that the steady-state model (14) has at least one nonconstant solution for all .

Corollary 13. Let be fixed. If and is odd, then there exists such that the steady-state model (14) has at least one nonconstant solution for any .

We omit the proofs of Corollaries 12 and 13 here and refer the reader to more detailed proofs in [52].

Remark 14. Our results suggest that if the parameters are properly chosen, both the general stationary pattern and more interesting Turing pattern can arise as a result of diffusion.

5. Turing Instability and Pattern Formation

5.1. Turing Instability

In this subsection, we mainly focus on the effects of reaction-diffusion on Turing instability for model (12).

Let us consider the spatially homogeneous system corresponding to model (12):

Mathematically speaking, a positive constant steady-state solution is Turing instability, which was emphasized by Turing in his pioneering work in 1952 [53], meaning that it is an asymptotically stable steady-state solution of model (86) but is unstable with respect to the solutions of spatial model (12). Since , then (the matrix ) is always true. Hence has an eigenvalue with a positive real part; then it must be a real value and the other eigenvalue must be a negative real one. As a consequence, we can obtain the following results.

Theorem 15. Assume that the following conditions are true:(a);(b);(c);
then the positive constant steady-state solution of model (12) is Turing instability if for some , where

Proof. Using the same approach as in Theorem 3, we have that the positive constant steady state solution of model (86) is asymptotically stable provided (a) and (b).
For the Turing instability, we must have for some . And we notice that achieves its minimum: at the critical value when .
We now discuss the conditions for Turing instability in the further calculation. The condition is equivalent to . In this case, has two positive roots and which satisfy Therefore, if we can find some such that , then . By [46] it follows that is uniformly asymptotically instable. This finishes the proof.

5.2. Pattern Formation

In this section, we perform extensive numerical simulations of the spatially extended model (12) in two-dimensional space, and the qualitative results are shown here. All our numerical simulations employ the zero-flux boundary conditions with a system size of .

The numerical integration of model (12) is performed by using a finite difference approximation for the spatial derivatives and an explicit Euler method for the time integration [54] with a time stepsize of . The initial condition is always a small amplitude random perturbation around the positive constant steady-state solution . After the initial period during which the perturbation spread, the model goes either into a time-dependent state, or to an essentially steady-state solution (time independent).

We use the standard five-point approximation [55] for the 2D Laplacian with the zero-flux boundary conditions. More precisely, the concentrations at the moment at the mesh position are given by with the Laplacian defined by where and , and the space stepsize .

In the numerical simulations, different types of dynamics are observed and it is found that the distributions of predator and prey are always of the same type. Consequently, we can restrict our analysis of pattern formation to one distribution. In this section, we show the distribution of prey , for instance. We have taken some snapshots with red (blue) corresponding to the high (low) value of prey .

Now, we show the Turing pattern for the different values of the parameters. Via numerical simulation, one can see that the model dynamics exhibits spatiotemporal complexity of pattern formation, including spots, stripes, and holes Turing patterns.

In Figure 1, we show the time process of spots pattern formation of the prey at for , , , , , , , and . In this case, there appears a competition between stripes and spots. The pattern takes a long time to settle down; starting with a homogeneous state (cf., Figure 1(a)), the random perturbations lead to the formation of stripes and spots (cf., Figure 1(b)), and the later random perturbations make these stripes decay and end with the time-independent regular spots pattern (cf., Figure 1(d)), which is isolated zones with low prey densities. Ecologically, spots pattern shows that the prey population is driven by predators to a very high level in those regions. The final result is the formation of patches of high prey density surrounded by areas of high prey densities; that is to say, under the control of these parameters, the prey is predominant in the area.

In Figure 2, we show the time process of stripes pattern formation of the prey at for . In this case, starting with a homogeneous state (cf., Figure 2(a)), the random perturbations lead to the formation of stripes and spots (cf., Figure 2(b)), and the later random perturbations make these stripes decay and end with the time-independent regular spots pattern (cf., Figure 2(d)).

In Figure 3, we show the time process of holes pattern formation of the prey at for , and . In this case, starting with a homogeneous state (cf., Figure 3(a)), the random perturbations lead to the formation of independent regular holes pattern (cf., Figure 1(d)), which is isolated zones with low prey densities.

Comparing Figure 1(d) with Figure 3(d), we find that they share similarities. Figure 1(d) consists of red (maximum density of ) spots on a blue (minimum density of ) background; that is, the preys are isolated zones with high population density. Figure 3(d) consists of blue (minimum density of ) spots on a red (maximum density of ) background; that is, the preys are isolated zones with low population density. For the sake of learning the pattern dynamics of model (12) further, we illustrate the three-dimensional patterns in the space . From Figure 4, one can realize the relations of the patterns (e.g., Figures 4(a) and 4(c)) with their corresponding numerical solutions (e.g., Figures 4(b) and 4(d)). In fact, the patterns (e.g., Figures 1ā€“3) are the projections in plane of the numerical solutions to the model (12). In the software MATLAB, imagesc(u) displays as an image; each numerical solution of corresponds to a rectangular area in the image; that is, the values of the numerical solutions to the model (12) are indices into the current colormap that determine the color of each patch.

6. Concluding Remarks

In the current investigation, we propose and analyze the dynamics of a reaction-diffusive Leslie-Gower predator-prey model with the additive Allee effect on prey. We are mainly concerned with the coexistence of the predator and prey and focus on the case of weak Allee effect (i.e., ). The value of this study is twofold. First, it shows the nonexistence and existence of the nonconstant positive steady-states, which guarantees the existence of Turing patterns. Second, it rigorously proves Turing instability by linear stability analysis and illustrates all three categories of Turing patterns close to the onset of Turing bifurcation via numerical simulations which indicates that the model dynamics exhibits complex pattern replication.

We summarize our findings as well as their related biological implications as follows.(1)Theorems 2 and 3 provide us with a full picture on the dynamics of the model with weak Allee effect. The dynamics of the model introduced can be very complicated due to Allee effects. In Theorem 2(a), we find that the trivial constant solution of the model which is subject to an Allee effect is unstable, which is exactly consistent with the model without Allee effect [24ā€“26]. This demonstrates that there is no extinction of a species in the present of Allee effects. By comparing them to their corresponding models without Allee effect, we can conclude that Allee effect can make the extinction of the prey although the maximum predation rate is small (see Theorem 2(b)). In this case, the predator species intends to change its food habits as predator has sufficient resources for alternative foods [56]. Furthermore, from Theorem 3, one can obtain that is locally uniformly asymptotically stable, which means that nonconstant positive solution (stationary pattern) of model (12) unlikely exists.(2)Theorems 9 and 11 indicate the existence and nonexistence of nonconstant steady-states with respect to various parameters. Roughly speaking, we can state that there is no pattern if the diffusion coefficients are suitably chosen. While pattern occurs provided that and all the eigenvalues have odd multiplicity (see Theorem 11). Hence, we can conclude that the multiplicity of patterns seems very interesting from the viewpoint of mathematics.(3)Theorem 15 and numerical simulations give us the existence of conditions of Turing instability and the types of Turing pattern. From the numerical results, one can see that our model has rich and complex spatiotemporal behavior. We find three typical Turing patterns, that is, stripes pattern, spot-stripe mixtures pattern, and spots pattern. That is to say, the effect of the Allee effect for pattern formation is tremendous. Therefore, the results of the present paper and [24, 25] show that the types of Turing pattern in the biological models depend on the effect of the Allee effect. In other words, the Allee effect may be one of the determining factors in producing spots and spot-stripe mixtures Turing patterns.

It is believed that our results made in this investigation related to predator-prey interactions due to the effect of Allee effect would certainly be of some help to theoretical mathematicians and the ecologists who are engaged in performing experimental work. Further studies are necessary to analyze the behaviour of a reaction-diffusion predator-prey model with the strong Allee effect (i.e., ).

Acknowledgments

The authors thank the anonymous referee and Professor Shangbin Cui for very helpful suggestions and comments which led to improvement of their original paper. This research was supported by the National Science Foundation of China (61373005, 11171357, and 11271290) and Zhejiang Provincial Natural Science Foundation (LY12A01014).