Complex Dynamics of a Diffusive Holling-Tanner Predator-Prey Model with the Allee Effect

Abstract

We investigate the complex dynamics of a diffusive Holling-Tanner predation model with the Allee effect on prey analytically and numerically. We examine the existence of the positive equilibria and the related dynamical behaviors of the model and find that when the model is with weak Allee effect, the solutions are local and global stability for some conditions around the positive equilibrium. In contrast, when the model is with strong Allee effect, this may lead to the phenomenon of bistability; that is to say, there is a separatrix curve that separates the behavior of trajectories of the system, implying that the model is highly sensitive to the initial conditions. Furthermore, we give the conditions of Turing instability and determine the Turing space in the parameters space. Based on these results, we perform a series of numerical simulations and find that the model exhibits complex pattern replication: spots, spots-stripes mixtures, and stripes patterns. The results show that the impact of the Allee effect essentially increases the models spatiotemporal complexity.

(ProQuest: ... denotes non-US-ASCII text omitted.)

Recommended by Lan Xu

Faculty of Science, Shaanxi University of Science and Technology, Xi'an, Shaanxi 710021, China

Received 2 August 2012; Revised 9 December 2012; Accepted 21 December 2012

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Recently, there has been a great interest in studying nonlinear difference/differential equations and systems [ 1- 6]. One of the reasons for this is a necessity for some techniques which can be used in investigating equations arising in mathematical models describing real-life situations in population biology, economy, probability theory, genetics, psychology, sociology, and so forth. And the bases for analyzing the dynamics of complex ecological systems are the interactions between two species, particularly the dynamical relationship between predators and their preys [ 7]. From the Lotka-Volterra model [ 8, 9], several alternatives for modeling continuous time consumer-resource interactions have been proposed. In recent years, one of the important predator-prey models is Holling-Tanner model, which was described by May [ 10]. This model reads as follows: [figure omitted; refer to PDF] where H and P represent prey and predator population densities at time t , respectively. r , K , _{c 1} , _{k 1} , _{s 1} , and δ are positive constants. _{r 1} and _{s 1} are the intrinsic growth rate of prey and predator, respectively. K is the carrying capacity of the prey, and δ takes on the role of the prey-dependent carrying capacity for the predator. The rate at which the predator consumes the prey, _{c 1} HP / ( _{k 1} +H ) , is known as the Holling type-II functional response [ 11]. The parameter _{c 1} is the maximum number of the prey that can be eaten per predator per time, and _{k 1} is the saturation value that corresponds to the number of the preys necessary to achieve one half of the maximum rate _{c 1} .

The dynamics of model ( 1) has been considered in many articles. For example, Hsu and Huang [ 12] obtained some results on the global stability of the positive equilibrium. More precise, under the conditions which local stability of the positive equilibrium implies its global stability. Gasull and coworkers [ 13] investigated the conditions of the asymptotic stability of the positive equilibrium which does not imply global stability. Sáez and González-Olivares [ 14] showed the asymptotic stability of a positive equilibrium and gave a qualitative description of the bifurcation curve.

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 conspecifics constitutes what is usually called an Allee effect [ 15- 24], starting with the pioneer work of Allee [ 25]. The outflux of prey to constant rate can be considered as Allee effect because a change on interaction dynamics is provoked, for instance, due to difficulty of encountering mates [ 17]. Nowadays, it is widely accepted that the Allee effect greatly increases the likelihood of local and global extinction [ 18] and can lead to a rich variety of dynamical effects.

From an ecological point of view, the Allee effect has been denominated in different ways [ 19- 22] and modeled into strong and weak ones [ 15, 16, 19], depending on the degree of positive density dependence. Mathematically speaking, if H =H (t ) indicates the population size, we assume that the growth function G (H ) satisfies the following:

(i) if G (0 ) =0 , ^{G [variant prime]} (0 ) >0 , G (H ) is called weak Allee effect;

(ii) if G (0 ) =0 , ^{G [variant prime]} (0 ) <0 , G (H ) is called strong Allee effect.

The most common mathematical form describing this phenomenon for a single species is given by [figure omitted; refer to PDF] where 0 <m <1 or -1 <m ...4;0 , which is named the multiplicative Allee effect; here, a threshold value m is incorporated such that population growth is negative below m . When m <H <1 , the per capita growth rate is positive.

Furthermore, Boukal et al. [ 22] proposed that the prey exhibits a demographic Allee effect at low population densities due to reasons other than predation by the focal predator as follows: [figure omitted; refer to PDF] where b is the Allee threshold, and c is an auxiliary parameter ( c >0 and b ...5; -c ). The auxiliary parameter c affects the overall shape of the per capita growth curve of the prey. When c is fixed, the unit growth rate of the species is only in connection with the Allee threshold.

For model ( 1), we make a change of variables as follows: [figure omitted; refer to PDF] For the sake of convenience, we still use variables H and P instead of H ~ and P ~ .

(H1) The basic model is a Holling-Tanner type as the form [figure omitted; refer to PDF] where a = _{c 1} /K , _{k 1} =K , and r = _{s 1} /r .

(H2) Following Boukal et al. [ 22], in Allee effect equation ( 3), we choose the auxiliary parameter c =1 , and b +1 ...1;m is the Allee threshold. That is, prey H has the population growth function [figure omitted; refer to PDF]

Obviously, we have the following:

(i) if 0 <m ...4;1 , G (0 ) =0 , ^{G [variant prime]} (0 ) >0 , the Allee effect ( 6) is the weak one;

(ii) if m >1 , G (0 ) =0 , ^{G [variant prime]} (0 ) <0 , the Allee effect ( 6) is the strong one;

(iii): if m =0 , the Allee effect will disappear.

And we can get the following model with the Allee effect on prey: [figure omitted; refer to PDF]

(H3) Assume that the individuals in populations H and P move randomly described as Brownian random motion [ 26]. We can get a simple spatial model corresponding to model ( 7) as follows: [figure omitted; refer to PDF] Here, the nonnegative constants _{d 1} and _{d 2} are the diffusion coefficients of H (t ) and P (t ) , respectively. ^{∇ 2} = ^{∂ 2} / ∂ ^{x 2} + ^{∂ 2} / ∂ ^{y 2} is the Laplacian operator in two-dimensional space, which describes the random moving. The initial distribution of species _{N 0} and _{P 0} are continuous functions. And the boundary condition is assumed to be zero-flux one as follows: [figure omitted; refer to PDF] L indicates the size of the model in the directions of x and y , respectively, and n is the outward unit normal vector of the boundary ∂ Ω . The main reason for choosing such boundary conditions is that we are interested in the self-organization of pattern, and the zero-flux boundary conditions imply that no external input is imposed form exterior [ 27].

There are some excellent works on a Holling-Tanner model considering the diffusion [ 28- 33] and the references therein. In [ 28], Guan and co-workers studied the spatiotemporal dynamics of a modified version of the Leslie-Gower predator-prey model incorporating a prey refuge and showed that the model dynamics exhibits complex Turing pattern replication: stripes, cold/hot spots-stripes coexistence, and cold/hot spots patterns. Without the Allee effect, Peng and Wang [ 29, 30] analyzed the global stability of the unique positive constant steady state and established some results for the existence and nonexistence of positive nonconstant steady states. Wang et al. [ 31] considered the Turing and Hopf bifurcations of the equilibrium solutions. Liu and Xue [ 32] investigated the pattern formation and found that spots, black-eye, and labyrinthine patterns can be observed in the model. Chen and Shi [ 33] proved global stability of the unique constant equilibrium.

However, the research about the influence of Allee effect on pattern formation of diffusive Holling-Tanner model seems rare. The main purpose of this paper is to study dynamical behaviors of a Holling-Tanner predator-prey model with the Allee effect. We will determine how the Allee effect affects the dynamics of the model and focus on the stability of the positive steady state and bifurcation mechanism and patterns formation analysis of the model.

The rest of the paper is organized as follows. In Sections 2and 3, we present our main results about the stability and bifurcation analysis of the nonspatial model ( 7) and the spatial model ( 8), respectively. Especially, in regards to the spatial model ( 8) in Section 3, we will give the conditions of the Turing instability and determine the Turing space, and by performing a series of numerical simulations, we illustrate the emergence of different patterns. Finally, in Section 4, some conclusions and remarks are given.

2. Dynamics Analysis of the Nonspatial Model ( 7)

2.1. Boundedness

Now, we prove that all solutions are eventually bounded.

Theorem 1.

All the solutions of model ( 7) which are initiated in ^{... + 2} are uniformly bounded.

Proof.

Let H (t ) and P (t ) be any solution of model ( 7) with initial conditions (H (0 ) ,P (0 ) ) = ( _{H 0} , _{P 0} ) such that _{H 0} >0 , _{P 0} >0 . From the first equation of model ( 7), we have [figure omitted; refer to PDF] a standard comparison theorem shows that [figure omitted; refer to PDF] Then, from the second equation of model ( 7), we get dP /dt ...4;rP (1 - (P / δ ) ) , which implies that [figure omitted; refer to PDF]

Define the function W (t ) =H (t ) +P (t ) , differentiating both sides with respect to t ; we get [figure omitted; refer to PDF] Then, [figure omitted; refer to PDF]

Using the theory of differential inequality, for all t ...5;T ...5;0 , we have [figure omitted; refer to PDF] Hence, _{limsup t [arrow right] ∞} (H (t ) +P (t ) ) ...4;M . This completes the proof.

Remark 2.

In fact, if m ...5;2 , dH / dt <0 always holds, which means that the prey and predator will extinct. Hence, we will later only focus on the case of 0 ...4;m <2 .

Next, we will investigate the existence of equilibria and their local and global stability with respect to model ( 7).

2.2. Equilibria Analysis in the Case of the Strong Allee Effect (i.e., 1 <m <2 )

In this subsection, we consider the existence and stability of the equilibrium of model ( 7) with strong Allee effect; that is, 1 <m <2 .

We note that model ( 7) is not defined at the P - axis, particularly at the point (0,0 ) , but both isoclines pass through this point, and in this case, it is a point of particular interest [ 34]. The character of (0,0 ) can be obtained after rescaling the time in model ( 7) by t = τH (1 +H ) as follows: [figure omitted; refer to PDF]

Lemma 3.

The point (0,0 ) of model ( 16) has a hyperbolic and a parabolic sector [ 20, 35] determined for the line P = ( δ (m -1 +r ) / r ) H . That is, there exists a separatrix curve in the phase plane that divides the behavior of trajectories; the point (0,0 ) is then an attractor point for certain trajectories and a saddle point for others.

Proof.

As the Jacobian matrix of the point (0,0 ) for model ( 16) is the zero matrix, we follow the methodology used in [ 20, 35] given by the function [straight phi] (u ,v ) = (uv ,v ) = (H ,P ) . Then, we have that [figure omitted; refer to PDF] and rescaling the time by T =v τ , it becomes [figure omitted; refer to PDF]

Clearly, if v =0 , then dv /dT =0 . Moreover, du / dT =u ( (1 -m )u -r ( u - ( 1 / δ ) ) ) .

The singularities of model ( 18) are (0,0 ) and ( r / ( δ (m -1 +r ) ) ,0 ) ; that is, a separatrix straight exists in the phase plane uv , given by u = r / ( δ (m -1 +r ) ) . The Jacobian matrixes of (0,0 ) and ( r / ( δ (m -1 +r ) ) ,0 ) for model ( 18) are [figure omitted; refer to PDF]

Then, (0,0 ) is a hyperbolic saddle point, and ( r / ( δ (m -1 +r ) ) ,0 ) is an attractor point. Using the blowing down, the point (0,0 ) is a saddle node in model ( 16), and the line P = ( ( δ (m -1 +r ) ) / r ) H divides the behavior of trajectories on the phase plane. The proof is completed.

Moreover, it is easy to verify that model ( 7) always has two boundary equilibria _{E 0} = (m -1,0 ) and _{E 1} = (1,0 ) . And the behavior of model ( 7) around _{E 0} and _{E 1} is found as follows.

The Jacobian matrix of model ( 7) at the equilibrium _{E 0} = (m -1,0 ) takes the form [figure omitted; refer to PDF] Hence, the equilibrium _{E 0} = (m -1,0 ) is an unstable node point (nodal source).

The Jacobian matrix of model ( 7) at the equilibrium _{E 0} = (1,0 ) takes the form [figure omitted; refer to PDF] Hence, the equilibrium _{E 1} = (1,0 ) is a saddle point.

And model ( 7) has a positive equilibrium E = (H , δH ) , where H satisfies [figure omitted; refer to PDF] For simplicity, we consider A =m -a δ and B = (m -a δ ^{) 2} +4 (1 -m ) ; then, the two roots of ( 22) are given by [figure omitted; refer to PDF]

Lemma 4.

(i) Suppose that m -a δ >0 and 1 <m <2 .

(a) If ^{B 2} >0 holds, model ( 7) has two positive equilibria _{E +} = ( _{H +} , δ _{H +} ) and _{E -} = ( _{H -} , δ _{H -} ) .

(b) If ^{B 2} =0 holds, model ( 7) has a unique positive equilibrium _{E e} = ( _{H e} , δ _{H e} ) . Note that in this case _{H e} = _{H +} = _{H -} =A /2 = m -1 .

(c) If ^{B 2} <0 , model ( 7) has no positive equilibrium.

(ii) If m -a δ ...4;0 , model ( 7) has no positive equilibrium.

Let E = (N ,P ) be an arbitrary positive equilibrium. The Jacobian matrix of model ( 7) at the positive equilibrium E = (H , δH ) takes the form [figure omitted; refer to PDF] Then, we can get [figure omitted; refer to PDF]

We can see that the sign of det ( _{J E} ) is determined by [figure omitted; refer to PDF]

Thus, we can obtain [figure omitted; refer to PDF]

Hence, we obtain det ( _{J E +} ) >0 , det ( _{J E -} ) <0 , and det ( _{J E e} ) =0 . And the positive equilibrium _{E -} = ( _{H -} , δ _{H -} ) is a saddle point. The nature of the equilibrium point _{E +} is dependent on the sign of the trace of the Jacobian matrix evaluated in this point. Whether _{E +} is a node or a focus depends on the sign of (tr ( _{J E +} ) ^{) 2} -4det ( _{J E +} ) .

In the following results, we study the stability of the positive equilibrium _{E +} and the unique positive equilibria _{E e} .

Theorem 5.

Define [figure omitted; refer to PDF]

(a) If r > _{r +} , the positive equilibrium _{E +} = ( _{H +} , δ _{H +} ) is a locally asymptotically stable point;

(a1) if ^{( _{r +} -r ) 2} <4det ( _{J E +} ) , then _{E +} is a stable focus,

(a2) if ^{( _{r +} -r ) 2} >4det ( _{J E +} ) , then _{E +} is a stable node point.

(b) If r < _{r +} , the positive equilibrium _{E +} = ( _{H +} , δ _{H +} ) is an unstable point;

(b1) if ( _{r +} -r ^{) 2} <4det ( _{J E +} ) , then _{E +} is an unstable focus surrounded by a stable limit cycle,

(b2) if ( _{r +} -r ^{) 2} >4det ( _{J E +} ) , then _{E +} is an unstable node and the limit cycle disappears.

(c) A Hopf bifurcation occurs at r = _{r +} around the positive equilibrium _{E +} = ( _{H +} , δ _{H +} ) . That is to say, model ( 7) has at least one positive periodic orbit.

Proof.

Here, we only give the proof of the existence of Hopf bifurcation. It is easy to see that

(i) tr (J ( _{E +} ) ) _{| r = r +} =0 holds,

(ii) the characteristic equation is ^{λ 2} +det (J ( _{E +} ) ) _{| r = r +} =0 , whose roots are purely imaginary,

(iii): ( d /dr ) [tr (J ( _{E +} ) ) _{] r = r +} = -1 ...0;0 .

From the Poincaré-Andronov-Hopf Bifurcation Theorem [ 36], we know that model ( 7) undergoes a Hopf bifurcation at _{E +} as r passes through the value _{r +} . The proof is completed.

Figure 1illustrates the local stability of the positive equilibrium _{E +} and the separatrix curve generated by the stable manifold of the positive equilibrium _{E -} . The orbits initiating the right of the separatrix curve tend to _{E +} , while the orbits initiating the left of the separatrix curve tend to (0,0 ) that represents the extinction of the population. Figure 2illustrates a Hopf bifurcation situation of the model around _{E +} . The parameter values are given in the figures.

Figure 1: The phase portrait of model ( 7) with the strong Allee effect. The parameters are taken as a =0.25 , δ =0.8 , m =1.3 and r =0.15 . In this case, _{E 0} = (0.3,0 ) is an unstable node point, _{E 1} = (1,0 ) and _{E -} = (0.5 ,0.4 ) are saddle points; the positive equilibrium _{E +} = (0.6 ,0.48 ) is local asymptotically stable. There exists a separatrix curve determined by the stable manifold of the equilibrium point _{E -} . The dotted curves are the nullclines.

[figure omitted; refer to PDF]

Figure 2: The phase portrait of model ( 7) with the strong Allee effect. The parameters are taken as a =0.25 , δ =0.8 , m =1.3 , and r =0.0375 . The model enters into a Hopf bifurcation around _{E +} = (0.6 ,0.48 ) at r = _{r +} .

[figure omitted; refer to PDF]

Theorem 6.

The unique equilibrium point _{E e} = ( m -1 , ( 1 / δ ) m -1 ) is

(i) a nonhyperbolic attractor node, if and only if r > m -1 (1 - m -1 ^{) 2} / ( 1 + m -1 ) ;

(ii) a nonhyperbolic repellor node, if and only if r > m -1 (1 - m -1 ^{) 2} / ( 1 + m -1 ) ;

(iii): a cusp point, if and only if r = m -1 (1 - m -1 ^{) 2} / ( 1 + m -1 ) , and in this case, there exists a unique trajectory which attains the point _{E e} . And in this case, the point (0,0 ) is a global attractor.

Proof.

We have [figure omitted; refer to PDF] Hence (i) and (ii) hold.

Moreover, tr ( _{J E e} ) =0 , if and only if r = m -1 (1 - m -1 ^{) 2} / ( 1 - m +1 ) . In this case, we obtain the Jacobian matrix of ( 16) as follows: [figure omitted; refer to PDF] and the associate Jordan matrix is [figure omitted; refer to PDF] Then, the singularity ( m -1 , ( 1 / δ ) m -1 ) is a cusp point, since it is a point of codimension 2, and we have a Bogdanov-Takens bifurcation [ 37].

The cusp point is shown in Figure 3.

Figure 3: The phase portrait of model ( 7) with the strong Allee effect. The parameters are taken as a =0.2556936062 , δ =0.8 , m =1.3 , and r =0.07238979895 . _{E 0} = (0.3,0 ) is an unstable node point; _{E 1} = (1,0 ) is saddle point; the positive equilibrium _{E 2} = (0.2 ,0.4 ) is a cusp point. In this case, the point (0 ,0 ) is globally asymptotically stable. The dotted curves are the nullclines.

[figure omitted; refer to PDF]

2.3. Equilibria Analysis in the Case of the Weak Allee Effect (i.e., 0 <m ...4;1 )

In this subsection, we consider the stability of the equilibrium of model ( 7) with weak Allee effect ( 0 <m ...4;1 ).

It is easy to verify that model ( 7) always has one boundary equilibrium _{E 0} = (1,0 ) which is a saddle point and a positive equilibrium ^{E *} = ( ^{H *} , δ ^{H *} ) , where [figure omitted; refer to PDF]

From ( 24), we have [figure omitted; refer to PDF]

Hence, we have the following results on the stability of the positive equilibrium ^{E *} = ( ^{H *} , δ ^{H *} ) .

Theorem 7.

Define [figure omitted; refer to PDF]

(a) If r > ^{r *} , the positive equilibrium ^{E *} = ( ^{H *} , δ ^{H *} ) is a locally asymptotically stable point, and

(a1) if ( ^{r *} -r ^{) 2} <4det ( _{J ^{E *}} ) , then ^{E *} is a stable focus;

(a2) if ( ^{r *} -r ^{) 2} >4det ( _{J ^{E *}} ) , then ^{E *} is a stable node point.

(b) If r < ^{r *} , the positive equilibrium ^{E *} = ( ^{H *} , δ ^{H *} ) is an unstable point, and

(b1) if ( ^{r *} -r ^{) 2} <4det ( _{J ^{E *}} ) , then ^{E *} is an unstable focus surrounded by a stable limit cycle;

(b2) if ( ^{r *} -r ^{) 2} >4det ( _{J ^{E *}} ) , then ^{E *} is an unstable node and the limit cycle disappears.

(c) A Hopf bifurcation occurs at r = ^{r *} around the positive equilibrium ^{E *} = ( ^{H *} , δ ^{H *} ) . That is to say, model ( 7) has at least one positive periodic orbit.

In the following theorem, we study the global behavior of the positive equilibrium ^{E *} .

Theorem 8.

If 0 <m < 1 / ( 1 +a δ ) , the positive equilibrium ^{E *} = ( ^{H *} , δ ^{H *} ) is globally asymptotically stable.

Proof.

Construct the following Lyapunov function: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] Then, [figure omitted; refer to PDF]

Substituting the value of dH /dt and dP /dt from the model of ( 7), we obtained [figure omitted; refer to PDF]

Note that a δ ^{H *} = ^{H *} ( 1 - ^{H *} ) ( 1 + ^{H *} -m ) ; we obtain [figure omitted; refer to PDF] Hence, if 0 <m <1 / (1 +a δ ) , ^{H *} -m >0 which is equivalent to dV /dt <0 .

Hence, the positive equilibrium ^{E *} = ( ^{H *} , δ ^{H *} ) is globally asymptotically stable. This completes the proof.

Figure 4demonstrates the global stability situation of model ( 7) around ^{E *} . Figure 5illustrates a Hopf bifurcation situation of the model around ^{E *} . Figure 6shows a stable limit cycle around ^{E *} which is an unstable focus. The parameter values are given in the figures.

Figure 4: The phase portrait of model ( 7) with the weak Allee effect. The parameters are taken as a =0.25 , δ =0.8 , m =1.3 , and r =0.1 . In this case, _{E 0} = (1,0 ) is saddle point; the positive equilibrium ^{E *} = (0.3232928050 ,0.4849392075 ) is globally asymptotically stable. The dotted curves are the nullclines.

[figure omitted; refer to PDF]

Figure 5: The phase portrait of model ( 7) with the weak Allee effect. The parameters are taken as a =0.25 , δ =0.8 , m =1.3 , and r =0.025265 . The model enters into a Hopf bifurcation around ^{E *} = (0.3232928050 ,0.4849392075 ) at r = ^{r *} .

[figure omitted; refer to PDF]

Figure 6: The phase portrait of model ( 7) with the weak Allee effect. The parameters are taken as a =0.8 , δ =1.5 , m =0.75 , and r =0.025 . The positive equilibrium ^{E *} = (0.3232928050 ,0.4849392075 ) is an unstable focus surrounded by a stable limit cycle.

[figure omitted; refer to PDF]

3. Dynamics of the Spatial Model ( 8)

In this section, we will investigate the dynamics of the spatial model ( 8). As an example, we only focus on the positive equilibrium point ^{E *} = ( ^{H *} , δ ^{H *} ) in the case of weak Allee effect ( 0 <m <1 ).

3.1. Turing Instability

Mathematically speaking, an equilibrium is Turing instability (diffusion-driven instability) means that it is an asymptotically stable equilibrium ^{E *} of model ( 7) but is unstable with respect to the solutions of reaction-diffusion model ( 8).

In the presence of diffusion, we will introduce small perturbations _{U 1} =H - ^{H *} , _{U 2} =P - δ ^{H *} , where | _{U 1} | , | _{U 2} | ...a;1 . To study the effect of diffusion on the model, we have considered the linearized form of model as follows: [figure omitted; refer to PDF] where ^{r *} is defined as ( 34).

Following Malchow et al. [ 38], we can know that any solution of model ( 40) can be expanded into a Fourier series so that [figure omitted; refer to PDF] where r = (x ,y ) , and 0 <x <L and 0 <y <L . k = ( _{k n} , _{k m} ) , and _{k n} =n π /L and _{k m} =m π /L are the corresponding wavenumbers.

Having substituted _{u nm} and _{v nm} into ( 40), we obtain [figure omitted; refer to PDF] where ^{k 2} = ^{k n 2} + ^{k m 2} .

A general solution of ( 42) has the form _{C 1} exp ( _{λ 1} t ) + _{C 2} exp ( _{λ 2} t ) , where the constants _{C 1} and _{C 2} are determined by the initial conditions ( 3), and the exponents _{λ 1} and _{λ 2} are the eigenvalues of the following matrix: [figure omitted; refer to PDF]

Correspondingly, _{λ 1} and _{λ 2} are the solutions of the following equation: [figure omitted; refer to PDF] where [figure omitted; refer to PDF]

Summarizing the previous discussions, we can get the following theorem immediately.

Theorem 9.

(i) The positive equilibrium ^{E *} of model ( 8) is locally asymptotically stable if r >max { ^{r *} , ^{r *} _{d 2} / _{d 1} } holds.

(ii) If the positive equilibrium ^{E *} of model ( 7) is globally asymptotically stable, then the corresponding steady state ^{E *} of model ( 8) is also globally asymptotically stable.

Proof.

(i) Using Routh-Hurwitz criteria, we can know that the positive equilibrium ^{E *} is locally asymptotically stable, if and only if tr ( D ~ ) <0 and det ( D ~ ) >0 . So, we obtain r >max { ^{r *} , ^{r *} _{d 2} / _{d 1} } .

(ii) We select the Lyapunov function for model ( 8) as follows: [figure omitted; refer to PDF] where V (H ,P ) is the same as defined in ( 35). So, [figure omitted; refer to PDF]

Using Green's first identity in the plane, [figure omitted; refer to PDF] Considering the zero-flux boundary conditions, one can show that [figure omitted; refer to PDF]

From the previous analysis, we note that d _{V 2} /dt <0 is valid if dV /dt <0 is true. This implies that the equilibrium ^{E *} of both model ( 7) and model ( 8) is globally asymptotically stable if 0 <m <1 / (1 +a δ ) holds. This ends the proof.

On the other hand, Turing instability sets in when at least one of the conditions is either tr ( D ~ ) <0 or det ( D ~ ) >0 . It is evident that the condition tr ( D ~ ) <0 is not violated when the requirement ^{r *} -r <0 is met [ 39]. Hence, only violation of condition det ( D ~ ) >0 gives rise to diffusion instability. Then, the condition for diffusive instability is given by [figure omitted; refer to PDF]

G ( · ) is quadratic in ^{k 2} , and the graph of G ( ^{k 2} ) =0 is a parabola. The minimum of G ( ^{k 2} ) occurs at ^{k 2} = ^{k m 2} , where [figure omitted; refer to PDF] The critical wave number _{k c} of the first perturbations to grow is found by evaluating _{k m} from ( 51).

Thus, a sufficient condition for Turing instability is that G ( ^{k m 2} ) is negative. Therefore, [figure omitted; refer to PDF] Combination of ( 51) and ( 52) leads to the following final criterion for diffusive instability: [figure omitted; refer to PDF]

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

Theorem 10.

If r _{d 1} / _{d 2} < ^{r *} <r and ^{r *} _{d 2} -r _{d 1} >2 _{d 1 d 2} det (J ( ^{E *} ) ) hold, the criterion for Turing instability for model ( 8) emerges, and the critical wave number _{k c} = ( ^{r *} _{d 2} -r _{d 1} ) / 2 _{d 1 d 2} .

The Turing instability (or bifurcation) breaks spatial symmetry, leading to the formation of patterns that are stationary in time and oscillatory in space [ 40, 41]. We adopt the intrinsic growth rates of predator r as the bifurcation parameter, and the linear stability analysis yields the bifurcation diagram shown in Figure 7. The Turing bifurcation curve separates the parametric space into two domains. Above the curve, the solutions of model ( 8) are stable for all pairs of (m ,r ) ; that is, there is no Turing instability. While below the curve, the solutions of model ( 8) are unstable for (m ,r ) and diffusive instability emerges; that is, Turing patterns emerge. This domain is called the Turing space.

Figure 7: Turing bifurcation diagram for model ( 8) using m and r as parameters. Other parameters are taken as a =0.8 , δ =1.75 , _{d 1} =0.01 , and _{d 2} =1 . Above the curve, the positive equilibrium ^{E *} is the only stable solution of model ( 3). Below the curve, the positive equilibrium ^{E *} loses its stability with respect to model ( 3), and Turing instability occurs; this domain is called the Turing space.

[figure omitted; refer to PDF]

3.2. Pattern Formation

In this subsection, we performed extensive numerical simulations of the spatially extended model ( 8) in two-dimension spaces, and the qualitative results are shown here. All our numerical simulations employ the zero-flux boundary conditions with a model size of L ×L , with L =100 discretized through x [arrow right] ( _{x 0} , _{x 1} , _{x 2} , ... , _{x n} ) and y [arrow right] ( _{y 0} , _{y 1} , _{y 2} , ... , _{y n} ) , with n =200 . Other parameters are fixed as a =0.8 , δ =1.75 , r =0.15 , _{d 1} =0.01 , and _{d 2} =1 .

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

More precisely, the concentrations ( ^{H i ,j n +1} , ^{P i ,j n +1} ) at the moment (n +1 ) τ at the mesh position (i ,j ) are given by [figure omitted; refer to PDF] with the Laplacian defined by [figure omitted; refer to PDF] where f (H ,P ) =H (1 -H ) ( 1 -m / (H +1 ) ) -aHP / (1 +H ) , g (H ,P ) =rP ( 1 -P / δH ) , and the space stepsize h =1 /3 .

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 H , for instance.

Figure 8shows the evolution of the spatial pattern of prey H at t =0 , 500 , 1000 , 2000 , with small random perturbation of the stationary solution ^{E *} of the spatially homogeneous systems when m is located in "Turing space." In this case, one can see that for model ( 8), the random initial distribution leads to the formation of a strongly irregular transient pattern in the domain. After the irregular pattern is formed (c.f., Figures 8(b)and 8(c)), it grows slightly and jumps alternately for a certain time, and finally spots patterns, which are isolated zones with low prey densities, prevail over the whole domain, and the dynamics of the model does not undergo any further changes (c.f., Figure 8(d)).

Spots pattern of H in model ( 8) for m =0.75 . Times: (a) 0, (b) 500, (c) 1000, and (d) 2000.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

(c) [figure omitted; refer to PDF]

(d) [figure omitted; refer to PDF]

Figure 9shows stripe patterns are interlaced stripes of high and low population densities of prey H for the parameter m =0.78 at t =1000 . In Figure 10, with the parameter m =0.85 , the spot-stripe mixtures pattern is time independent with low prey densities.

Figure 9: Stripes pattern of H in model ( 8) for m =0.78 at t =1000 .

[figure omitted; refer to PDF]

Figure 10: Spot-stripe mixtures pattern of H in model ( 8) for m =0.85 at t =1000 .

[figure omitted; refer to PDF]

4. Concluding Remarks

In this paper, we are concerned with the complex dynamics in a diffusive Holling-Tanner predator-prey model with the Allee effect on prey. The value of this study lies in two folds. First, the local asymptotic stability conditions for coexisting equilibrium and conditions for Hopf bifurcation are described briefly for the model with the weak and strong Allee effects. Second, it gives the analysis of Turing instability which determines the Turing space in the spatial domain and meanwhile illustrates the Turing pattern formation close to the onset Turing bifurcation via numerical simulations, which shows that the model dynamics exhibits complex pattern replication.

We note that in the analyzed models, a big difference between the dynamics of model with strong or weak Allee effect exists. In the case of strong Allee effect, two positive equilibria can coexist for a subset of parameters with a varied dynamics but different to other Holling-Tanner models analyzed earlier [ 12- 14]. We have shown that one of these equilibria is always a saddle point and proved the existence of a separatrix curve. And there is no global asymptotically stable positive equilibrium. In this case, the point (0,0 ) is an attractor in addition to locally stable positive equilibrium _{E +} for determined parameter values, which leads to the existence of a bistability phenomenon. The dynamics of the model is determined by the initial conditions; the predator and prey may be extinction or coexistence. This means that the strong Allee effect could easily lead to the risk of population extinction. Nevertheless, in the case of weak Allee effect, model ( 7) can only have one unique positive equilibrium, which is globally asymptotically stable under some conditions. Therefore, the predators and preys can coexist in stable conditions.

Furthermore, we have investigated the conditions for the predator-prey model which experiences spatial patterns through diffusion-driven instability. We have derived the conditions of Turing instability in terms of our model parameters analytically. In addition, to get a deeper insight into the model's dynamics behaviour, we select the different values of parameter m . An increase of m , from the numerical results, one can see that our model has rich and complex spatiotemporal behavior. We find three typical Turing patterns, that is, spots pattern, stripes pattern, and spots-stripes mixtures pattern. To the best of our knowledge, the Turing pattern we illustrate here is the first reported case to our model. And our complete analysis of the spatial model will give new suggestion to the models with the Allee effect.

Acknowledgments

The authors would like to thank the anonymous referee for very helpful suggestions and comments which led to improvements of their original paper. And this work is supported by the Cooperative Project of Yulin City (2011).

References

[1] J. D. Murray Mathematical Biology. II: Spatial Models and Biomedical Applications , vol. 18, of Interdisciplinary Applied Mathematics, Springer, New York, NY, USA, 2003., 3rd.

[2] S. Stevic, "A short proof of the Cushing-Henson conjecture," Discrete Dynamics in Nature and Society , 2006.

[3] S. Stevic, "Asymptotic behavior of a class of nonlinear difference equations," Discrete Dynamics in Nature and Society , 2006.

[4] M. De La Sen, S. Alonso-Quesada, "Model-matching-based control of the Beverton-Holt equation in ecology," Discrete Dynamics in Nature and Society , 2008.

[5] M. De la Sen, S. Alonso-Quesada, "A control theory point of view on Beverton-Holt equation in population dynamics and some of its generalizations," Applied Mathematics and Computation , vol. 199, no. 2, pp. 464-481, 2008.

[6] J. He, "Asymptotic methods for solitary solutions and compactons," Abstract and Applied Analysis , vol. 2012, 2012.

[7] Y. Li, D. Xiao, "Bifurcations of a predator-prey system of Holling and Leslie types," Chaos, Solitons & Fractals , vol. 34, no. 2, pp. 606-620, 2007.

[8] A. J. Lotka Elements of Physical Biology , Williams and Wilkins, Baltimore, Md, USA, 1925.

[9] V. Volterra, "Fluctuations in the abundance of a species considered mathematically," Nature , vol. 118, no. 2972, pp. 558-560, 1926.

[10] R. M. May Stability and Complexity in Model Ecosystems , Princeton University Press, Princeton, NJ, USA, 1974.

[11] C. S. Holling, "The components of predation as revealed by a study of small-mammal predation of the European Pine Sawfly," The Canadian Entomologist , vol. 91, no. 5, pp. 293-320, 1959.

[12] S. B. Hsu, T. W. Huang, "Global stability for a class of predator-prey systems," SIAM Journal on Applied Mathematics , vol. 55, no. 3, pp. 763-783, 1995.

[13] A. Gasull, R. E. Kooij, J. Torregrosa, "Limit cycles in the Holling-Tanner model," Publicacions Matematiques , vol. 41, no. 1, pp. 149-167, 1997.

[14] E. Sáez, E. González-Olivares, "Dynamics of a predator-prey model," SIAM Journal on Applied Mathematics , vol. 59, no. 5, pp. 1867-1878, 1999.

[15] G. Wang, X. G. Liang, F. Z. Wang, "The competitive dynamics of populations subject to an Allee effect," Ecological Modelling , vol. 124, no. 2-3, pp. 183-192, 1999.

[16] A. Kent, C. P. Doncaster, T. Sluckin, "Consequences for predators of rescue and Allee effects on prey," Ecological Modelling , vol. 162, no. 3, pp. 233-245, 2003.

[17] E. Gonzalez-Olivares, J. Mena-Lorca, H. Meneses-Alcay, B. Gonzalez-Yanez, J. D. Flores, "Allee effect, emigration and immigration in a class of predator-prey models," Biophysical Reviews and Letters , vol. 3, pp. 195-251, 2008.

[18] S. R. Zhou, Y. F. Liu, G. Wang, "The stability of predator-prey systems subject to the Allee effects," Theoretical Population Biology , vol. 67, no. 1, pp. 23-31, 2005.

[19] J. Wang, J. Shi, J. Wei, "Predator-prey system with strong Allee effect in prey," Journal of Mathematical Biology , vol. 62, no. 3, pp. 291-331, 2011.

[20] E. González-Olivares, J. Mena-Lorca, A. Rojas-Palma, J. D. Flores, "Dynamical complexities in the Leslie-Gower predator-prey model as consequences of the Allee effect on prey," Applied Mathematical Modelling , vol. 35, no. 1, pp. 366-381, 2011.

[21] P. Aguirre, E. González-Olivares, E. Sáez, "Three limit cycles in a Leslie-Gower predator-prey model with additive Allee effect," SIAM Journal on Applied Mathematics , vol. 69, no. 5, pp. 1244-1262, 2009.

[22] D. S. Boukal, M. W. Sabelis, L. Berec, "How predator functional responses and Allee effects in prey affect the paradox of enrichment and population collapses," Theoretical Population Biology , vol. 72, no. 1, pp. 136-147, 2007.

[23] Y. Cai, W. Wang, J. Wang, "Dynamics of a diffusive predator-prey model with additive Allee effect," International Journal of Biomathematics , vol. 5, no. 2, 2012.

[24] A. Morozov, S. Petrovskii, B. L. Li, "Bifurcations and chaos in a predator-prey system with the Allee effect," Proceedings of the Royal Society B , vol. 271, no. 1546, pp. 1407-1414, 2004.

[25] W. C. Allee Animal Aggregations: A Study in General Sociology , University of Chicago Press, Chicago, Ill, USA, 1931.

[26] R. S. Cantrell, C. Cosner Spatial Ecology via Reaction-Diffusion Equations , of Wiley Series in Mathematical and Computational Biology, John Wiley & Sons, Chichester, UK, 2003.

[27] J. D. Murray, "Discussion: turing's theory of morphogenesis-its influence on modelling biological pattern and form," Bulletin of Mathematical Biology , vol. 52, no. 1-2, pp. 119-132, 1990.

[28] X. Guan, W. Wang, Y. Cai, "Spatiotemporal dynamics of a Leslie-Gower predator-prey model incorporating a prey refuge," Nonlinear Analysis: Real World Applications , vol. 12, no. 4, pp. 2385-2395, 2011.

[29] R. Peng, M. Wang, "Positive steady states of the Holling-Tanner prey-predator model with diffusion," Proceedings of the Royal Society of Edinburgh A , vol. 135, no. 1, pp. 149-164, 2005.

[30] R. Peng, M. Wang, "Global stability of the equilibrium of a diffusive Holling-Tanner prey-predator model," Applied Mathematics Letters , vol. 20, no. 6, pp. 664-670, 2007.

[31] W. Wang, Y. Lin, L. Zhang, F. Rao, Y. Tan, "Complex patterns in a predator-prey model with self and cross-diffusion," Communications in Nonlinear Science and Numerical Simulation , vol. 16, no. 4, pp. 2006-2015, 2011.

[32] P. Liu, Y. Xue, "Spatiotemporal dynamics of a predator-prey model," Nonlinear Dynamics , vol. 69, no. 1-2, pp. 71-77, 2012.

[33] S. Chen, J. Shi, "Global stability in a diffusive Holling-Tanner predator-prey model," Applied Mathematics Letters , vol. 25, no. 3, pp. 614-618, 2012.

[34] J. Mena-Lorca, E. Gonzalez-Olivares, B. Gonzalez-Yanez, "The Leslie-Gower predator-prey model with Allee effect on prey: a simple model with a rich and interesting dynamiscs," in Proceedings of the International Symposium on Mathematical and Computational Biology, pp. 105-132, 2006.

[35] E. Sáez, E. González-Olivares, "Dynamics of a predator-prey model," SIAM Journal on Applied Mathematics , vol. 59, no. 5, pp. 1867-1878, 1999.

[36] S. Wiggins Introduction to Applied Nonlinear Dynamical Systems and Chaos , vol. 2, of Texts in Applied Mathematics, Springer, New York, NY, USA, 1990.

[37] F. Dumortier, J. Llibre, J. C. Artés Qualitative Theory of Planar Differential Systems , of Universitext, Springer, Berlin, Germany, 2006.

[38] H. Malchow, S. V. Petrovskii, E. Venturino Spatiotemporal Patterns in Ecology and Epidemiology: Theory, Models, and Simulation , of Chapman & Hall/CRC Mathematical and Computational Biology Series, Chapman & Hall/CRC, Boca Raton, Fla, USA, 2008.

[39] A. Okubo, S. A. Levin Diffusion and Ecological Problems: Modern Perspectives , vol. 14, of Interdisciplinary Applied Mathematics: Mathematical Biology, Springer, New York, NY, USA, 2001., 2nd.

[40] W. Wang, Q.-X. Liu, Z. Jin, "Spatiotemporal complexity of a ratio-dependent predator-prey system," Physical Review E , vol. 75, no. 5, 2007.

[41] W. Wang, L. Zhang, H. Wang, Z. Li, "Pattern formation of a predator-prey system with Ivlev-type functional response," Ecological Modelling , vol. 221, no. 2, pp. 131-140, 2010.

[42] M. R. Garvie, "Finite-difference schemes for reaction-diffusion equations modeling predator-prey interactions in MATLAB," Bulletin of Mathematical Biology , vol. 69, no. 3, pp. 931-956, 2007.

[]