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.
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: where and represent prey and predator population densities at time , respectively. , , , , , and are positive constants. and are the intrinsic growth rate of prey and predator, respectively. 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, , is known as the Holling type-II functional response [11]. The parameter is the maximum number of the prey that can be eaten per predator per time, and is the saturation value that corresponds to the number of the preys necessary to achieve one half of the maximum rate .
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 indicates the population size, we assume that the growth function satisfies the following:(i)if , , is called weak Allee effect;(ii)if , , is called strong Allee effect.
The most common mathematical form describing this phenomenon for a single species is given by where or , which is named the multiplicative Allee effect; here, a threshold value is incorporated such that population growth is negative below . When , 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: where is the Allee threshold, and is an auxiliary parameter ( and ). The auxiliary parameter affects the overall shape of the per capita growth curve of the prey. When 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: For the sake of convenience, we still use variables and instead of and .
(H1) The basic model is a Holling-Tanner type as the form where , , and .
(H2) Following Boukal et al. [22], in Allee effect equation (3), we choose the auxiliary parameter , and is the Allee threshold. That is, prey has the population growth function
Obviously, we have the following:(i)if , , , the Allee effect (6) is the weak one;(ii)if , , , the Allee effect (6) is the strong one;(iii)if , the Allee effect will disappear.
And we can get the following model with the Allee effect on prey:
(H3) Assume that the individuals in populations and move randomly described as Brownian random motion [26]. We can get a simple spatial model corresponding to model (7) as follows: Here, the nonnegative constants and are the diffusion coefficients of and , respectively. is the Laplacian operator in two-dimensional space, which describes the random moving. The initial distribution of species and are continuous functions. And the boundary condition is assumed to be zero-flux one as follows: indicates the size of the model in the directions of and , respectively, and 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 2 and 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 are uniformly bounded.
Proof. Let and be any solution of model (7) with initial conditions such that , . From the first equation of model (7), we have
a standard comparison theorem shows that
Then, from the second equation of model (7), we get , which implies that
Define the function , differentiating both sides with respect to ; we get
Then,
Using the theory of differential inequality, for all , we have
Hence, . This completes the proof.
Remark 2. In fact, if , always holds, which means that the prey and predator will extinct. Hence, we will later only focus on the case of .
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., )
In this subsection, we consider the existence and stability of the equilibrium of model (7) with strong Allee effect; that is, .
We note that model (7) is not defined at the axis, particularly at the point , but both isoclines pass through this point, and in this case, it is a point of particular interest [34]. The character of can be obtained after rescaling the time in model (7) by as follows:
Lemma 3. The point of model (16) has a hyperbolic and a parabolic sector [20, 35] determined for the line . That is, there exists a separatrix curve in the phase plane that divides the behavior of trajectories; the point is then an attractor point for certain trajectories and a saddle point for others.
Proof. As the Jacobian matrix of the point for model (16) is the zero matrix, we follow the methodology used in [20, 35] given by the function . Then, we have that
and rescaling the time by , it becomes
Clearly, if , then . Moreover, .
The singularities of model (18) are and ; that is, a separatrix straight exists in the phase plane , given by . The Jacobian matrixes of and for model (18) are
Then, is a hyperbolic saddle point, and is an attractor point. Using the blowing down, the point is a saddle node in model (16), and the line 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 and . And the behavior of model (7) around and is found as follows.
The Jacobian matrix of model (7) at the equilibrium takes the form Hence, the equilibrium is an unstable node point (nodal source).
The Jacobian matrix of model (7) at the equilibrium takes the form Hence, the equilibrium is a saddle point.
And model (7) has a positive equilibrium , where satisfies For simplicity, we consider and ; then, the two roots of (22) are given by
Lemma 4. (i)Suppose that and .(a)If holds, model (7) has two positive equilibria and .(b)If holds, model (7) has a unique positive equilibrium . Note that in this case .(c)If , model (7) has no positive equilibrium.(ii)If , model (7) has no positive equilibrium.
Let be an arbitrary positive equilibrium. The Jacobian matrix of model (7) at the positive equilibrium takes the form Then, we can get
We can see that the sign of is determined by
Thus, we can obtain
Hence, we obtain , , and . And the positive equilibrium is a saddle point. The nature of the equilibrium point is dependent on the sign of the trace of the Jacobian matrix evaluated in this point. Whether is a node or a focus depends on the sign of .
In the following results, we study the stability of the positive equilibrium and the unique positive equilibria .
Theorem 5. Define (a) If , the positive equilibrium is a locally asymptotically stable point;(a1) if , then is a stable focus,(a2) if , then is a stable node point.(b) If , the positive equilibrium is an unstable point;(b1) if , then is an unstable focus surrounded by a stable limit cycle,(b2) if , then is an unstable node and the limit cycle disappears.(c) A Hopf bifurcation occurs at around the positive equilibrium . 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) holds,(ii)the characteristic equation is , whose roots are purely imaginary,(iii).
From the Poincaré-Andronov-Hopf Bifurcation Theorem [36], we know that model (7) undergoes a Hopf bifurcation at as passes through the value . The proof is completed.
Figure 1 illustrates the local stability of the positive equilibrium and the separatrix curve generated by the stable manifold of the positive equilibrium . The orbits initiating the right of the separatrix curve tend to , while the orbits initiating the left of the separatrix curve tend to that represents the extinction of the population. Figure 2 illustrates a Hopf bifurcation situation of the model around . The parameter values are given in the figures.
Theorem 6. The unique equilibrium point is(i)a nonhyperbolic attractor node, if and only if ;(ii)a nonhyperbolic repellor node, if and only if ;(iii)a cusp point, if and only if , and in this case, there exists a unique trajectory which attains the point . And in this case, the point is a global attractor.
Proof. We have
Hence (i) and (ii) hold.
Moreover, , if and only if . In this case, we obtain the Jacobian matrix of (16) as follows:
and the associate Jordan matrix is
Then, the singularity 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.
2.3. Equilibria Analysis in the Case of the Weak Allee Effect (i.e., )
In this subsection, we consider the stability of the equilibrium of model (7) with weak Allee effect ().
It is easy to verify that model (7) always has one boundary equilibrium which is a saddle point and a positive equilibrium , where
From (24), we have
Hence, we have the following results on the stability of the positive equilibrium .
Theorem 7. Define (a) If , the positive equilibrium is a locally asymptotically stable point, and(a1) if , then is a stable focus;(a2) if , then is a stable node point.(b) If , the positive equilibrium is an unstable point, and(b1) if , then is an unstable focus surrounded by a stable limit cycle;(b2) if , then is an unstable node and the limit cycle disappears.(c) A Hopf bifurcation occurs at around the positive equilibrium . 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 .
Theorem 8. If , the positive equilibrium is globally asymptotically stable.
Proof. Construct the following Lyapunov function:
where
Then,
Substituting the value of and from the model of (7), we obtained
Note that ; we obtain
Hence, if , which is equivalent to .
Hence, the positive equilibrium is globally asymptotically stable. This completes the proof.
Figure 4 demonstrates the global stability situation of model (7) around . Figure 5 illustrates a Hopf bifurcation situation of the model around . Figure 6 shows a stable limit cycle around which is an unstable focus. The parameter values are given in the figures.
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 in the case of weak Allee effect ().
3.1. Turing Instability
Mathematically speaking, an equilibrium is Turing instability (diffusion-driven instability) means that it is an asymptotically stable equilibrium 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 , , where . To study the effect of diffusion on the model, we have considered the linearized form of model as follows: where 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 where , and and . , and and are the corresponding wavenumbers.
Having substituted and into (40), we obtain where .
A general solution of (42) has the form , where the constants and are determined by the initial conditions (3), and the exponents and are the eigenvalues of the following matrix:
Correspondingly, and are the solutions of the following equation: where
Summarizing the previous discussions, we can get the following theorem immediately.
Theorem 9. (i) The positive equilibrium of model (8) is locally asymptotically stable if holds.
(ii) If the positive equilibrium of model (7) is globally asymptotically stable, then the corresponding steady state of model (8) is also globally asymptotically stable.
Proof. (i) Using Routh-Hurwitz criteria, we can know that the positive equilibrium is locally asymptotically stable, if and only if and . So, we obtain .
(ii) We select the Lyapunov function for model (8) as follows:
where is the same as defined in (35). So,
Using Green’s first identity in the plane,
Considering the zero-flux boundary conditions, one can show that
From the previous analysis, we note that is valid if is true. This implies that the equilibrium of both model (7) and model (8) is globally asymptotically stable if holds. This ends the proof.
On the other hand, Turing instability sets in when at least one of the conditions is either or . It is evident that the condition is not violated when the requirement is met [39]. Hence, only violation of condition gives rise to diffusion instability. Then, the condition for diffusive instability is given by
is quadratic in , and the graph of is a parabola. The minimum of occurs at , where The critical wave number of the first perturbations to grow is found by evaluating from (51).
Thus, a sufficient condition for Turing instability is that is negative. Therefore, Combination of (51) and (52) leads to the following final criterion for diffusive instability:
Summarizing the previous discussions, we can obtain the following theorem.
Theorem 10. If and hold, the criterion for Turing instability for model (8) emerges, and the critical wave number .
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 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 ; that is, there is no Turing instability. While below the curve, the solutions of model (8) are unstable for and diffusive instability emerges; that is, Turing patterns emerge. This domain is called the Turing space.
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 , with discretized through and , with . Other parameters are fixed as , , , , and .
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 . 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, either the model goes into a time-dependent state or to an essentially steady state solution (time independent).
More precisely, the concentrations at the moment at the mesh position are given by with the Laplacian defined by where , , 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.
Figure 8 shows the evolution of the spatial pattern of prey at , , , , with small random perturbation of the stationary solution of the spatially homogeneous systems when 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)).
(a)
(b)
(c)
(d)
Figure 9 shows stripe patterns are interlaced stripes of high and low population densities of prey for the parameter at . In Figure 10, with the parameter , the spot-stripe mixtures pattern is time independent with low prey densities.
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 is an attractor in addition to locally stable positive equilibrium 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 . An increase of , 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).