Self-Organized Patterns Induced by Neimark-Sacker, Flip and Turing Bifurcations in a Discrete Predator-Prey Model with Lesie-Gower Functional Response

Abstract

The formation of self-organized patterns in predator-prey models has been a very hot topic recently. The dynamics of these models, bifurcations and pattern formations are so complex that studies are urgently needed. In this research, we transformed a continuous predator-prey model with Lesie-Gower functional response into a discrete model. Fixed points and stability analyses were studied. Around the stable fixed point, bifurcation analyses including: flip, Neimark-Sacker and Turing bifurcation were done and bifurcation conditions were obtained. Based on these bifurcation conditions, parameters values were selected to carry out numerical simulations on pattern formation. The simulation results showed that Neimark-Sacker bifurcation induced spots, spirals and transitional patterns from spots to spirals. Turing bifurcation induced labyrinth patterns and spirals coupled with mosaic patterns, while flip bifurcation induced many irregular complex patterns. Compared with former studies on continuous predator-prey model with Lesie-Gower functional response, our research on the discrete model demonstrated more complex dynamics and varieties of self-organized patterns.

1. Introduction

Predator-prey systems are some of the essential ecological systems in Nature. The dynamic behaviors of the predator-prey system have captured the interest of both biologists and ecologists [1,2,3,4,5]. There are a substantial number of predator-prey system models. Recently, the formation of patterns has become a very hot topic [6,7,8,9]. This is because the formation process demonstrates self-organization of spatial heterogeneity, and shows system complexity directly and visibly. This visible complexity matches well what has been found in real ecosystems, however, the dynamics of predator-prey systems are so complex that more studies are still needed to explore the mechanism of pattern formation.

Reaction-diffusion models are commonly used to investigate spatially extended predator-prey systems [6,7,8,10,11,12,13,14]. Through Turing bifurcation, reaction-diffusion models have successfully revealed the pattern formation mechanism in many different ecosystems [2,12,13,15,16,17]. Simulations have shown a variety of patterns such as spots, stripes, labyrinth, spirals, gaps, and so on [6,10,13]. Different reaction-diffusion models focus on different functional responses that describe different predation relationships. There are a variety of responses [18], such as Lesie-Gower functional responses [19,20], Beddington-DeAngelis functional responses [14,21,22,23,24] and so on. Among these functional responses, the Leslie-Gower formulation is based on the assumption that a reduction in a predator population has a reciprocal relationship with per capita availability of its preferred food. Indeed, Leslie [19] introduced a predator-prey model where the carrying capacity of the predator is proportional to the number of its prey. This interesting formulation for predator dynamics has been discussed by Leslie and Gower [20] and Pielou [25].

In expressing predator-prey models, two main types of models are considered—continuous models and discrete models—with the latter showing more advantages in revealing complex nonlinear dynamics. A classic example is the discrete logistic model that exhibits period-doubling cascade and a route to chaos [21]. On the contrary, the continuous logistic model shows a simple “S” form curve, but never demonstrates the above dynamic complexity. Through a variety of bifurcations, discrete models can generate periodic orbits, invariant circles, periodic windows, chaotic behavior and so on. More importantly, in predator-prey systems, births and deaths of bionts are discrete events, which means continuous models only make sense for very large populations. Many researchers think that discrete models can reveal the discontinuous properties (such as a patchy environment or a fragmented habitat) of predator-prey systems [26]. Besides, discrete models may exhibit new dynamic behaviors. For example, Han et al. found that Turing instability and Turing patterns can occur in a simple discrete competitive Lotka-Volterra system rather than the continuous one [27].

In this research, we will transform a well-recognized continuous predator-prey model into a discrete model. The continuous model incorporates the Holling-type-II and the modified Lesie-Gower functional responses, that can be shown as:

$$\{\begin{array}{c}\frac{\partial H}{\partial T}=(a_1-b_{1}H-\frac{c_{1}P}{H+k_1})H+D_1(\frac{{\partial }^{2}H}{\partial X^2}+\frac{{\partial }^{2}H}{\partial Y^2})\\ \frac{\partial P}{\partial T}=(a_2-\frac{c_{2}P}{H+k_2})P+D_2(\frac{{\partial }^{2}P}{\partial X^2}+\frac{{\partial }^{2}P}{\partial Y^2})\end{array}$$

(1)

where H and P are the densities of prey and predators, respectively. (X,Y) is the spatial position of species when they move in a two dimensional space. D₁ and D₂ are the diffusion coefficients of prey and predator respectively. a₁ is the growth rate of prey H. a₂ describes the growth rate of predator P. b₁ measures the strength of competition among individuals of species H. c₁ is the maximum value of the per capita reduction of H due to P. c₂ has a similar meaning to c₁. k₁ measures the extent to which environment provides protection to prey H. k₂ has a similar meaning to k₁ relatively to the predator P. With scaling transformations:

$$\begin{array}{c}t=a_{1}T,u(t)=\frac{b_1}{a_1}H(T),v(t)=\frac{c_2{b}_1}{a_1{a}_2}P(T),x=X{(\frac{a_1}{D_1})}^{\frac{1}{2}},y=Y{(\frac{a_1}{D_1})}^{\frac{1}{2}}\\ a=\frac{a_2{c}_1}{a_1{c}_2},b=\frac{a_2}{a_1},e_1=\frac{b_1{k}_1}{a_1},e_2=\frac{b_1{k}_2}{r_1},\delta =\frac{D_2}{D_1},\end{array}$$

(2)

the continuous model can be expressed as:

$$\{\begin{array}{c}\frac{\partial u}{\partial t}=\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}+u(1-u)-\frac{auv}{u+e_1}=\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}+f(u,v)\\ \frac{\partial v}{\partial t}=\delta (\frac{{\partial }^{2}v}{\partial x^2}+\frac{{\partial }^{2}v}{\partial y^2})+b(1-\frac{v}{u+e_2})v=\delta (\frac{{\partial }^{2}v}{\partial x^2}+\frac{{\partial }^{2}v}{\partial y^2})+g(u,v)\end{array}.$$

(3)

The local dynamics of this continuous model (3) has been studied in [28,29], and the global stability of a similar model has been investigated in [30]. A similar model with delay is studied in [31,32], and a three dimensional similar model with the same functional responses is studied in [33,34,35]. In [36], several Turing and Hopf bifurcation patterns were obtained by the continuous model (3). We can see that this model is well recognized.

A few researchers found that more complex dynamics could be generated through discretizing the continuous model [26,37,38,39,40]. Based on the approach in prior studies [26,37,38,39], here we investigated the complex dynamics via transforming the continuous model (3) to a discrete model. In the analysis of the discrete model, fixed points and stability analysis were studied. Three types of bifurcation analysis were conducted, including flip bifurcation, Neimark-Sacker bifurcation and Turing bifurcation. Then numerical simulations were carried out under these bifurcation conditions, to show the complex dynamics and the formation self-organized patterns with this discrete model. Discussions focused on the types of self-organized patterns, and the relations between bifurcations and pattern types.

2. Model and Stability Analysis

2.1. A Discrete Predator-Prey Model

In this research, the above continuous model Equation (3) will be transformed to a discrete model. We consider the model on a N × N lattice, and the two variables can be expressed as u_(i,j,t) and v_(i,j,t) (i,j ∈ {1,2,3,…N} and t ∈ Z⁺), that represent the prey density and the predator density in lattice (i,j) at time t, respectively. According to the former research works of [26,37,38,39], there are two stages, reaction stage and diffusion stage, when we discretize the continuous model (3). The diffusion stage is considered firstly as:

$$\{\begin{array}{c}u{\prime }_{(i,j,t)}=u_{(i,j,t)}+\frac{\tau }{h^2}{\nabla }_d^2{u}_{(i,j,t)}\\ v{\prime }_{(i,j,t)}=v_{(i,j,t)}+\frac{\tau }{h^2}\delta {\nabla }_d^2{v}_{(i,j,t)}\end{array},$$

(4)

where τ and h are the time interval and space interval. ${\nabla }_d^2$ denotes the discrete form of the Laplacian operator. Note that the whole research is using periodic boundary conditions. Then we consider the reaction stage:

$$\{\begin{array}{c}{u}_{(i,j,t+1)}=f_1(u{\prime }_{(i,j,t)},v{\prime }_{(i,j,t)})\\ v_{(i,j,t+1)}=g_1(u{\prime }_{(i,j,t)},v{\prime }_{(i,j,t)})\end{array},$$

(5)

in which:

$$\{\begin{array}{c}{f}_1(u,v)=u+\tau f(u,v)\\ g_1(u,v)=v+\tau g(u,v)\end{array}.$$

(6)

Equations (4)–(6) including both diffusion and reaction stages are defined as our discrete model.

2.2. Fixed Points and Stability

As we obtain a new discrete model, fixed points and stability analysis need to be investigated. The fixed points should satisfy:

$${\nabla }_d^2{u}_{(i,j,t)}={\nabla }_d^2{v}_{(i,j,t)}=0.$$

(7)

Then we get:

$$\{\begin{array}{c}u{\prime }_{(i,j,t)}=u_{(i,j,t)}\\ v{\prime }_{(i,j,t)}=v_{(i,j,t)}\end{array},$$

(8)

and:

$$\{\begin{array}{c}{u}_{(i,j,t+1)}=f_1(u_{(i,j,t)},v_{(i,j,t)})\\ v_{(i,j,t+1)}=g_1(u_{(i,j,t)},v_{(i,j,t)})\end{array}.$$

(9)

The above equation can be expressed as the following map:

$$\left(\begin{array}{c}u\\ v\end{array}\right)\to \left(\begin{array}{c}u+\tau u(1-u)-\frac{\tau auv}{u+e_1}\\ v+\tau b(1-\frac{v}{u+e_2})v\end{array}\right).$$

(10)

Thus the fixed points of the map can be shown as:

$$(u_1,v_1)=(0,0),$$

(11a)

$$(u_2,v_2)=(0,e_2),$$

(11b)

$$(u_3,v_3)=(1,0),$$

(11c)

$$(u_4,v_4)=(u^{\ast },v^{\ast }).$$

(11d)

in which, (u*,v*) is the positive solution of the following equations:

$$\{\begin{array}{c}{u}^2+(a+e_1-1)u+a{e}_2-e_1=0\\ v=u+e_2\end{array}.$$

(12)

Therefore, the sufficient and necessary condition of u* > 0, v* > 0 can be obtained as:

$${(a+e_1-1)}^2-4(a{e}_2-e_1)>0,1-a-e_1+\sqrt{(a+e_1-1)-4(a{e}_2-e_1)}>0.$$

(13)

The stability of the fixed points can be studied through a Jacobi matrix with spatially homogeneous perturbations, that is shown as:

$$J=\left[\begin{array}{cc}1+\tau [1-2u-\frac{av{e}_1}{{(u+e_1)}^2}]& \frac{-\tau au}{u+e_1}\\ \frac{\tau b{v}^2}{{(u+e_2)}^2}& 1+\tau (b-\frac{2bv}{u+e_2})\end{array}\right].$$

(14)

Substituting (u₁,v₁) into J, we can get:

$${J|}_{(u_1,v_1)}=\left[\begin{array}{cc}1+\tau & 0\\ 0& 1+\tau b\end{array}\right].$$

(15)

As 1 + τ > 1, 1 + τb > 1, thus (u₁,v₁) is unstable. Substituting $(u_2,v_2)$ into $J$, we can get:

$${J|}_{(u_2,v_2)}=\left[\begin{array}{cc}1+\tau (1-\frac{a{e}_2}{e_1})& 0\\ \tau b& 1-\tau b\end{array}\right].$$

(16)

When 1 < ae₂/e₁ < 1 + 2/τ and 0 < b < 2/τ, (u₂,v₂) is stable. Substituting (u₃,v₃) into J, we can get:

$${J|}_{(u_3,v_3)}=\left[\begin{array}{cc}1-\tau & -\frac{\tau a}{1+e_1}\\ 0& 1+\tau b\end{array}\right].$$

(17)

The stability requires $\left|1-\tau \right|<1$ and $\left|1+\tau b\right|<1$, that is $0<\tau <2$ and $\tau <0$. But this is a contradiction. Thus $(u_3,v_3)$ is unstable.

Substituting $(u_4,v_4)$ into $J$, we can get:

$${J|}_{(u_4,v_4)}=\left[\begin{array}{cc}{a}_{11}& a_{12}\\ a_{21}& a_{22}\end{array}\right].$$

(18)

In which:

$$a_{11}=1+\tau [1-2u^{\ast }-\frac{a(u^{\ast }+e_2)e_1}{{(u^{\ast }+e_1)}^2}],a_{12}=\frac{-\tau a{u}^{\ast }}{u^{\ast }+e_1},a_{21}=\tau b,a_{22}=1-\tau b$$

(19)

Let:

$$t{r}_0=a_{11}+a_{22},$$

(20)

$${\Delta }_0=a_{11}{a}_{22}-a_{12}{a}_{21}.$$

(21)

The eigenvalues can be obtained as:

$${\lambda }_{1,2}=\frac{t{r}_0\pm \sqrt{t{r_0}^2-4{\Delta }_0}}{2}.$$

(22)

Therefore $(u_4,v_4)$ is stable only when:

$$\{\begin{array}{c}{\Delta }_0<1\\ -(1+{\Delta }_0)<t{r}_0<1+{\Delta }_0\end{array}.$$

(23)

2.3. Bifurcation Analysis

In this subsection, we focus on the bifurcation analysis at fixed point (u₄,v₄). Parameter conditions of three bifurcations are obtained, including Neimark-Sacker bifurcation, flip bifurcation and Turing bifurcation. The detailed progress of each bifurcation analysis is shown in Appendix A.

3. Numerical Simulations

Simulations will be carried out for each bifurcation calculated in Section 2.3. (Appendix A). Bifurcation diagrams and phase portraits will be shown to interpret the system dynamics when Neimark-Sacker bifurcation or flip bifurcation conditions occur. The emergence of Turing bifurcation is also demonstrated by simulating the variations of eigenvalues (curve of Equation (A43)). And then self-organized patterns will be shown with the corresponding parameters under each bifurcation condition.

3.1. Bifurcation Diagram and Phase Portrait

Figure 1 shows variations of u versus the parameter τ when the parameter values satisfy the flip bifurcation conditions. When τ > 2.3267, the fixed point is asymptotically stable. When τ = 2.3267, the system starts to bifurcate around a fixed point. With the increase of $\tau$, the stable states of the system go through (not only these states) period-2 (τ = 2.5 as shown in Figure 1b), period-4 (τ = 2.72 as shown in Figure 1c), period-10 (τ = 2.77 as shown in Figure 1d) and then complex periodic oscillations (τ = 2.83 as shown in Figure 1e).

Figure 2 shows the variations of u versus parameter τ when the parameter values satisfy the Neimark-Sacker bifurcation conditions. Note that this bifurcation diagram has no visible periodic windows. When τ < 2.8291, the fixed point is asymptotically stable. When τ = 2.8291, the system starts to bifurcate around a fixed point. With the increase of τ, the states of the system increase, but the system always follows the invariant circle. For example, when τ = 3, Figure 2b shows how the bifurcation drives the system from a fixed point (with 1% random perturbations) to the invariant circle in the phase plane (u,v).

Figure 3 shows another Neimark-Sacker bifurcation situation. Note that this bifurcation diagram has several periodic windows. When τ < 3.1087, the fixed point is asymptotically stable and the bifurcation point of the system is determined at τ = 2.1087.

As the value of τ increases, the stable states of the system experience several stages, such as invariant circle (τ = 3.1750 as shown in Figure 3b), period-6 (τ = 3.5250 as shown in Figure 3c), and then invariant circle again (τ = 3.69 as shown in Figure 3d).

Turing bifurcation can be shown through the variations of eigenvalues λ(k,l) as shown in Figure 4. In Figure 4a, we can see that the effects of the perturbation numbers k and l are symmetric. Thus we let k = l, and we can get the variations of eigenvalues versus l as shown in Figure 4b. When there is no perturbation, the system stays at the fixed point. When the diffusion coefficient δ = 2, the eigenvalues of the system remain at less than 1 with the increase of perturbation number l. But when the diffusion coefficient δ increases more than δ_c, the eigenvalues will exceed 1 with the increase of l and Turing bifurcation occurs.

3.2. Formation of Self-Organized Patterns

In order to investigate the formation of self-organized patterns under the above bifurcation conditions, simulations will be carried out on 100 × 100 lattices with periodic boundary conditions. Initial conditions are set as fixed points with heterogeneous random disturbance (1%). Given the parameter values under each bifurcation condition, the formation of patterns can be obtained after t = 2000. Patterns will be shown only in terms of variable u, as the patterns of variable v are similar.

Figure 5 shows the patterns induced by Neimark-Sacker bifurcation. The parameters a = 1.1, e₁ = 0.3, e₂ = 0.2, b = 0.3739 and τ = 2.84 satisfy the Neimark-Sacker bifurcation conditions. We can see that spot patterns are formed through self-organization of variable u. With the increase of τ, the patterns become more isolated.

Figure 6 shows the patterns of variable $u$ induced by Neimark-Sacker bifurcation. Parameters a = 1.1, e₁ = 0.3, e₂ = 0.2, b = 0.373 and τ satisfy the Neimark-Sacker bifurcation conditions. We can see that with the increase of τ, patterns are formed through the self-organization of u. Moreover, these patterns transit from the irregular pattern in Figure 6a to spiral patterns in Figure 6b–d, and the wavelength of the spirals decreases gradually. Spiral patterns are some of the patterns that are often recorded in the studies of predator-prey systems [41].

Figure 7 shows the patterns of variable $u$ induced by Turing bifurcation. Parameters a = 1.1, e₁ = 0.3, e₂ = 0.2, b = 0.1781 and τ = 0.1, 0.2 do not satisfy Neimark-Sacker bifurcations or flip bifurcation conditions, but the addition of δ = 30 and h = 6 make Turing bifurcation occur. We can see that a labyrinth pattern is formed through the self-organization of u. With the increase of τ, the stripes in the labyrinth pattern become broader.

Figure 8 shows the pattern of variable u induced by Turing bifurcation. Parameters a = 1.1, e₁ = 0.3, e₂ = 0.2, b = 0.3739 and τ = 2.8, 2.82 do not satisfy Neimark-Sacker bifurcations or flip bifurcation conditions, but addition of δ = 3.2 and h = 6 make Turing bifurcation occur. We can see that complex patterns are formed through the self-organization of u. The parameter values of Figure 7 and Figure 8 both satisfy the Turing bifurcation conditions, but the types are quite different from each other. The patterns in Figure 8 are like spiral patterns coupled with irregular mosaics. Note that the patterns in Figure 8 are similar as those in Figure 6. The reason may be that the value of τ is very close to the Neimark-Sacker bifurcation point, that is τ = 2.8291.

Figure 9 shows the pattern induced by flip bifurcation. Parameters a = 1.1, e₁ = 0.3, e₂ = 0.2, b = 1.2 and τ satisfy the flip bifurcation conditions. We can see that patterns are formed through the self-organization of variable u. The type of pattern is quite difficult to define, but they can be reflected by the phase portrait above, that are similar to those in Figure 1. Figure 9a–d show the patterns generated by period-2, period-4, period-6 and multi-period behaviors. Due to the effect of period-doubling cascade, patterns become more and more fractal.

4. Discussion and Conclusions

We have transformed a well-recognized continuous predator-prey model into a discrete model. During the transformation, the iteration process is considered in two stages: diffusion stage and then reaction stage. Fixed points and stability analysis are studied. Three types of bifurcation analysis are performed for the discrete model, including flip bifurcation, Neimark-Sacker bifurcation, and Turing bifurcation. With parameter values satisfying bifurcation conditions, numerical simulations are carried out to show the self-organized patterns of the predator-prey model.

From the pattern simulations, we can see that Turing instability leads to the formation of spatial patterns of labyrinth and complex patterns of spirals coupled with mosaics, Neimark-Sacker bifurcation induces the formation of spots and spiral patterns, and flip bifurcation induces the formation of irregular patterns corresponding to period-2, period-4 and many complex periodic behaviors. We can see that these patterns are quite different from each other, which shows the distinguishing effects of these three types of instabilities (based on three bifurcation situations).

These self-organized patterns can also be classified into two categories: stationary patterns and oscillatory patterns. Oscillatory patterns include the spiral patterns, and stationary patterns include the other patterns. In the stationary patterns, the spatial distribution of predator and prey remains stationary and the system dynamics will not change with time. In the oscillatory patterns, the dynamics of predator and prey is always varying spatially and temporally as time goes on. One of important characteristics of the oscillatory patterns is spatiotemporal chaos, which results in the formation of complex and diverse spiral patterns.

Compared with previous studies on a continuous predator-prey model with Lesie-Gower functional response [36], the simulations on our discrete model with same functional response show not only that all the patterns were obtained by the continuous model, but also many special patterns can be induced by flip bifurcation and a special type of pattern is induced by Turing bifurcation. This proves that discrete models can generate more self-organized patterns. The most important reason is, especially in our model, the import of time interval τ gives more possibility to the system to generate these complex dynamics. Furthermore, we think the variations of τ reflect the multiple time scale in the real ecosystems. The conclusions of this research can be summarized as follows:

• The discrete predator-prey model with Lesie-Gower functional response can generate many complex dynamics including three types of bifurcations, which are flip bifurcation, Neimark-Sacker bifurcation and Turing bifurcation.
• A variety of self-organized patterns can be formed through the discrete predator-prey model with Lesie-Gower functional response and the above three bifurcations. These patterns consist of spots, transitional patterns from spots to spirals, spirals, spirals coupled with mosaics, labyrinths, and many other complex patterns generated by flip bifurcation.
• Among the studies on predator-prey models with Lesie-Gower functional response, this research may develop a special perspective to interpret the how self-organized patterns are generated.