Spatiotemporal patterns induced by four mechanisms in a tussock sedge model with discrete time and space variables

In this paper, we investigate the spatiotemporal patterns of a freshwater tussock sedge model with discrete time and space variables. We first analyze the kinetic system and show the parametric conditions for flip and Neimark–Sacker bifurcations respectively. With spatial diffusion, we then show that the obtained stable homogeneous solutions can experience Turing instability under certain conditions. Through numerical simulations, we find periodic doubling cascade, periodic window, invariant cycles, chaotic behaviors, and some interesting spatial patterns, which are induced by four mechanisms: pure-Turing instability, flip-Turing instability, Neimark–Sacker–Turing instability, and chaos.


Introduction
Tussock sedge is a common North American carex stricata. They grow almost exclusively on the top of tussocks in dense tufted clumps and have wiry stems and leaves. The spaces between the clumps can provide the habitat for small animals. Tussock sedge plays an important role in wetland ecosystem, such as dike reinforcement, slope protection, and soil erosion prevention. In general, the dead wracks lie at the foot of the tussock, gather in the low inter-tussock areas, and inhibit the growth of plant. On the other hand, they activate the wracks. The interaction between the plant and the wrack as well as the unequal diffusions of the two are essential for the formation of spatiotemporal patterns [1]. For the theoretical significance and application value in the real world, it is very interesting to investigate the formation mechanism of spatiotemporal patterns from mathematical point of view. In [2], the authors studied the spatial patterns of the tussock sedge by setting an experimental site and found that the regular patterns have a close connection with the space scale. The models they studied are characterized by reaction diffusion equations. In recent years, many biological and ecological phenomena, especially for spatial patterns, have been characterized by reaction diffusion equations. Different diffusion forms, crossdiffusion [3], self-diffusion, sub-diffusion [4], super-diffusion [5,6], and so on, for example, are widely used in reaction diffusion equations, which are the key factors that lead the system to undergo Turing instability [1] and form abundant patterns: spots, stripes, labyrinth, gaps, spirals, circles, and so forth [7][8][9][10].
In recent decades, lots of works about Turing instability concentrating on continuous or semi-discrete reaction diffusion systems have been carried out, please see [11,12] and the references therein. For one of the tussock sedge models in [2] with continuous time and space variables, the authors in [13] found layered Turing patterns which are generated by an equilibrium and a limit cycle. In [14], the authors showed that the formation of patterns is related to the domain size, the growth rate of the tussock sedge grasses, and the carrying capacity of land. Moreover, they derived the conditions for the existence of the patterns by computing the Leray-Schauder degree. In [15], the author considered another plant-wrack model in a two-dimensional space, but with half inhibition in the plant biomass equation, and found spotted, labyrinth, and coexistence of stripe-like and spotted patterns. In many situations, the data acquisition is not time continuous, and the distribution of biological population is also spatially discontinuous. Furthermore, the algorithm of simulating spatial patterns is on the basis of the corresponding discrete form of the continuous system. Based on these considerations, a discrete model or an effective discretization method is also very important in connecting the actual model with the simulations. Therefore, in this article, we mainly concentrate on the pattern formation mechanisms of a tussock sedge model in [2] but with the discrete time and space variables. We will apply the coupled map lattices (CMLs) method to discretize the continuous tussock sedge model to get the corresponding system with discrete time and space variables. The system obtained by the CMLs method is usually called CMLs model. Surprisingly, dynamics and patterns of the CMLs model are more abundant compared with a continuous system. For the CMLs model, some effective methods, such as the Jacobian matrix for stability analysis, center manifold reduction [16], and calculation of spatial discrete operator [17], and theoretical results about Turing instability criterion for discrete systems [18] have been obtained and made accessible to researchers [19][20][21][22][23]. An advantage of the CMLs model is that it retains the inherent properties of the original system; in the meantime, it represents a numerical simulation algorithm [24,25]. The CMLs method has been widely used in many aspects, such as chemical oscillator [26], ecological system [23], neural dynamic system [27,28]. For more related literature works, please see [7,[29][30][31][32] and the references therein.
There is not much research on the dynamic behavior and pattern formation for the tussock sedge model of CMLs type. Through a series of theoretical analyses, we find that the model shows abundant dynamics. If given suitable parameters, the CMLs model can undergo spatiotemporal bifurcation simultaneously, which leads to four kinds of pattern formation mechanisms: flip-Turing bifurcation, Neimark-Sacker-Turing bifurcation, chaotic oscillation, and pure-Turing bifurcation. For this reason, the patterns are more abundant compared with the continuous form.
We organize the paper as follows. In Sect. 2, we first develop the corresponding CMLs model in Sect. 2.1, and then implement the theoretical analysis on the dynamics of the homogeneous stationary state for the kinetic system in the reaction stage in Sect. 2.2. In Sect. 3, we mainly consider the bifurcation behaviors of the homogeneous stationary state, including the flip bifurcation in Sect. 3.1, Neimark-Sacker bifurcation in Sect. 3.2, and Turing bifurcation in Sect. 3.3. In Sect. 4, simulations are carried out for the purpose to explain the theoretical conclusions and illustrate the dynamics and spatial pattern. Finally, we further show our conclusions and discussion.

The CMLs type model and stability analysis
We first construct the CMLs type model corresponding to the continuous tussock sedge model, and then analyze its stability.

The development of CMLs tussock sedge model
One of the continuous tussock sedge models in [2] is where P(x, y, t) and W (x, y, t) are the plant and wrack biomass at time t and space (x, y) respectively, and = ∂ 2 ∂x 2 + ∂ 2 ∂y 2 . The term s represents a specific rate of leaf senescence. The decomposition rate of wrack is denoted by b. The diffusion rates of plant and wrack are d p and d w respectively, and d w > d p . In this model, the wrack has an inhibition on the growth of plant, namely the linear term -θ PW , where θ represents an inhibition rate. For more information about model (1), please see [2].
In order to keep the notation concise and intuitive, we drop -, replace T with t, and rewrite system (2) into the following dimensionless form: (3) Based on system (3), we build up our CMLs model. On a two-dimensional rectangular area, we define some n × n lattice sites. We specify two numbers, the plant biomass P (i,j,t) and the wrack biomass W (i,j,t) , at time t ∈ Z + on each site (i, j), (i, j ∈ {1, 2, 3, . . . , n}). Between different sites, we presume that there are local reactions and spatial diffusions [23,25], and the biomass of plant and wrack at each site follow the system dynamics with time.
From t to t + 1, the dynamical behaviors of plant and wrack for the CMLs model contain two stages: "reaction" stage and "diffusion" stage [7,[23][24][25]. The diffusion behavior occurs before the reaction behavior. Introducing the time step τ and space step δ and discretizing the spatial form of system (3), we get the equations governing the dispersal process: where P (i,j,t) and W (i,j,t) are the biomass of plant and wrack that take part in the next reaction stage. The Laplacian operator in discrete form is described by Via discretizing the non-spatial form of (3), we achieve the equations that control the reaction process: , where Equations (4)-(7) are the CMLs model of system (3). All the parameters are positive, and P (i,j,t) ≥ 0, W (i,j,t) ≥ 0. In the next, we investigate the CMLs model's dynamics under the following periodic boundary conditions: P (0,j,t) = P (n,j,t) , P (1,j,t) = P (n+1,j,t) , W (0,j,t) = W (n,j,t) , W (1,j,t) = W (n+1,j,t) .
The discrete time and space tussock sedge model have spatially homogeneous and heterogeneous behaviors. For all i, j, and t, the homogeneous behavior satisfies Together with equations (4)-(7), the homogeneous dynamics, ignoring the spatial sites index, are governed by which can be written into the following maps: The homogeneous dynamics of CMLs model (4)-(7) then can be obtained via analyzing map (10). For heterogeneous behavior, it requires at least one group of i, j, and t to make ∇ 2 d P (i,j,t) = 0 and ∇ 2 d W (i,j,t) = 0.

Dynamics of the homogenous stationary state
In this subsection, we try to get the parametric conditions that make the homogeneous stationary state stable. Since the fixed point of map (10) is exactly equivalent to the homogeneous stationary state of CMLs model (4)-(7), we can analyze the dynamics of the fixed point of map (10) instead of analyzing the homogeneous stationary state of the CMLs model. By solving the following equations: we get two fixed points of map (10): (P 1 , W 1 ) = (0, 0) and (P 2 , W 2 ) = ( s 1+ξ , s 1+ξ ). For the stability of the two fixed points, the following theorem gives a series of precise parametric conditions.

Theorem 1
The fixed point (P 1 , W 1 ) of map (10) is a saddle. For the fixed point (P 2 , W 2 ): (1) If one of (H 1 ) and (H 2 ) holds, it is a saddle, where (2) If one of (SN 1 ) and (SN 2 ) holds, it is a stable node, and if one of (UN 1 ) and (UN 2 ) holds, it is an unstable node, where Furthermore, if one of (H 3 ) and (H 4 ) holds, it is a stable degenerate node, and if one of (H 5 ) and (H 6 ) holds, it is unstable, where (3) If (SF) holds, it is a stable focus, and if (UF) holds, it is an unstable focus, where Proof By computing the eigenvalues of the Jacobian matrix of map (10) at the associated fixed points respectively, we can determine the stability conditions of the fixed points. For the fixed point (P 1 , W 1 ) = (0, 0), the Jacobian matrix is Obviously, the eigenvalues of J(P 1 , W 1 ) are 1 + τ s > 1 and 1τ < 1. According to [33], the fixed point (P 1 , W 1 ) is a saddle. For the fixed point (P 2 , W 2 ), the corresponding Jacobian matrix is The two eigenvalues are According to [33], we can obtain the results by direct calculation.
In the sequel, the bifurcation behaviors and pattern formation mechanisms around the fixed point (P 2 , W 2 ) are investigated, since the fixed point (0, 0) is always unstable. For description convenience, we call the fixed point of map (10) the homogeneous stationary state of CMLs model (4)-(7).

Bifurcation analysis of the homogeneous stationary state
In this section, we treat τ as the main bifurcation parameter and explore the bifurcation behavior of the homogeneous stationary state, such as flip bifurcation, Neimark-Sacker bifurcation, and Turing bifurcation respectively. Based on this theoretical analysis, we get the parametric conditions to support the formation of spatial patterns.

Analysis of flip bifurcation
The flip bifurcation can lead the stable node (P 2 , W 2 ) to be an unstable one. In the meantime, new period-2 points appear. The first requirement for the emergence of flip bifurcation is that λ 1 or λ 2 for the Jacobian matrix J(P 2 , W 2 ) at the critical value of flip bifurcation is -1, and the other is neither -1 nor 1. This requirement can be satisfied by the following conditions: Solving the above conditions, we can get the critical value of flip bifurcation and arrive at the following conclusion.
If the re-scaled conversion rate s satisfies 0 < s < s 1 or s > s 2 , as τ increases from less than τ 1 to more than τ 1 , the positive homogeneous stationary state (P 2 , W 2 ) changes from stable node to unstable node. Considering η 1 = 0 and η 2 > 0, CMLs model (4)-(7) demonstrates flip bifurcation at (P 2 , W 2 ), and a stable period-2 orbit is bifurcated on the right-hand side of the critical value of the flip bifurcation τ 1 . This implies that the plant and wrack biomass may coexist in an oscillatory state between the period-2 points.
Under condition (17), we make the transformation w = P -P 2 , z = W -W 2 to translate the fixed point (P 2 , W 2 ) to the origin (0, 0), which can simplify the latter description and analysis. In this coordinate transformation, map (10) can be written as follows: The coefficients a 1 , a 2 , a 3 , a 4 , b 1 , b 2 are as previously defined in (12), but with τ 1 replaced with τ * . Because the inverse translation will not change the qualitative behavior of the fixed point, the eigenvalues of the Jacobian matrix of map (18) at (0, 0) are also conjugate complex numbers and the modulus is 1. For notational convenience, we denote the two eigenvalues by λ(τ * ) andλ(τ * ), and where J(τ * ) = J(P 2 , W 2 )| τ =τ * , i 2 = -1 and |λ(τ * )| = |λ(τ * )| = 1. Except for condition (17), the Neimark-Sacker bifurcation also needs the eigenvalues to satisfy the nonzero transversality condition d|λ(τ * )| dτ = 0. Direct computations show that Moreover, the Neimark-Sacker bifurcation requires that the two eigenvalues are neither real nor imaginary. This requirement can be guaranteed by the following condition: namely, For the sake of the Neimark-Sacker bifurcation, we also need to obtain the normal form of map (18) by carrying out the center manifold reduction to get the last criterion. Let w = a 2w , z = (αa 1 )wβz, this invertible transformation leads map (18) tõ where To ensure the emergence of Neimark-Sacker bifurcation for map (22), we demand the determinative quantity σ satisfying where a 2a 1 ) , , Based on the above analysis and computation, we have the following.
If the re-scaled conversion rate s satisfies s 1 < s < s 2 , as 0 < τ < τ * , the homogeneous stationary state (P 2 , W 2 ) is a stable focus, when τ > τ * , it is an unstable focus. The dynamic transition implies the possibility of the occurrence of Neimark-Sacker bifurcation. Considering σ = 0, this indicates that the Neimark-Sacker bifurcation emerges when τ = τ * . If σ < 0, the bifurcated cycle is stable. This implies that the plant and wrack coexist in the form of uniform quasi-periodic oscillation.

Analysis of Turing bifurcation
The breaking of spatial symmetry is the main reason for the occurrence of Turing bifurcation. If Turing instability occurs, the stable homogeneous stationary state of the CMLs model is driven to be unstable due to the uneven spatial diffusions. This induces the formation of spatial patterns. For the occurrence of Turing bifurcation, two conditions are essential [7,8,34]. First we need the nontrivial homogeneous stationary state to be stable about time. Second we need the stable nontrivial homogeneous stationary state to be unstable under one or more kinds of spatially heterogeneous perturbations. From Theorem 1, we know that if one of conditions (SN 1 ), (SN 2 ), (SF), (H 3 ), and (H 4 ) holds, then (P 2 , W 2 ) is stable about time.
Throughout this subsection, we suppose that one of (SN 1 ), (SN 2 ), (SF), (H 3 ), and (H 4 ) holds, and explore the Turing bifurcation of the homogeneous stationary state. To get the conditions that support the Turing instability, we first discuss the eigenvalue problems of the discrete Laplacian operator ∇ 2 d . Taking into account the following equation: adapted to the periodic boundary conditions X i,0 = X i,n , X i,1 = X i,n+1 , X 0,j = X n,j , X 1,j = X n+1,j .
As in [17], we have that the eigenvalue λ kl of ∇ 2 d satisfies λ kl = 4 sin 2 (k -1)π n + sin 2 (l -1)π n := 4 sin 2 φ k + sin 2 φ l , where k, l ∈ {1, 2, 3, . . . , n}. Corresponding to the eigenvalue λ kl , the eigenfunction is denoted by X ij kl , namely, ∇ 2 d X ij kl + λ kl X ij kl = 0. In the sequel, we investigate the stability of (P 2 , W 2 ) under small spatial heterogenous perturbation. LetP (i,j,t) = P (i,j,t) -P 2 ,W (i,j,t) = W (i,j,t) -W 2 , and take the perturbation into CMLs model (4) where a 1 , a 2 , b 1 , and b 2 are given in (12). When the perturbation is small, the linear terms dominate the dynamics of system (26). Multiplying X ij kl at both ends of equation (26), we obtain X ij klP (i,j,t+1) = X ij kl (a 1P(i,j,t) + a 2W(i,j,t) ) For equation (27), we sum all i and j together to get the following system: LetP t+1 = n i,j=1 X ij klP (i,j,t) ,W t+1 = n i,j=1 X ij klW (i,j,t) , system (28) can be written as follows: At all the lattices, the dynamics of the spatially heterogeneous perturbations are governed by system (29). If the fixed point (0, 0) of system (29) is stable, then the spatially homogenous stationary state (P 2 , W 2 ) of CMLs model (4)- (7) is stable. Otherwise, the state (P 2 , W 2 ) is unstable. In the latter situation, Turing patterns will come into formation. In the following, we calculate the eigenvalues of the Jacobian matrix of system (29) at (0, 0). We are very interested in finding the parametric conditions to make at least one eigenvalue with module greater than 1. The Jacobian of system (29) at (0, 0) is .
The above computations signify the following theorem.
Concluding the whole analysis and computations, we conclude that if the homogeneous stationary state (P 2 , W 2 ) only experiences Turing instability, then the patterns are induced by pure-Turing bifurcation. If the flip bifurcation and Turing instability occur simultaneously, the patterns are induced by flip-Turing instability. If Neimark-Sacker bifurcation and Turing instability happen at the same time, the spatial patterns are caused by Neimark-Sacker-Turing instability.

Numerical simulation
We carry out numerical simulations to illustrate the dynamic evolvement of flip, Neimark-Sacker bifurcations, and Turing instability as well as the related spatiotemporal patterns in this section.

Flip bifurcation and the related Turing patterns
We first explain the theoretical results obtained in Proposition 2 and Theorem 3. Then we combine the flip bifurcation with the Turing bifurcation and present some spatial patterns induced by pure-Turing instability, flip-Turing instability, and chaos mechanisms respectively. In this subsection, we set ξ = 5.58, s = 5, and the computational rectangular grid is n = 150.
The location of the fixed point (P 2 , W 2 ) = ( s 1+ξ , s 1+ξ ) is independent of τ , but its stability is closely related with τ . The fixed point is (P 2 , W 2 ) = (0.7599, 0.7599). Further direct calculation shows that s 1 = 0.51769, s 2 = 2.68605. Obviously, the condition (C 2 ) : s > s 2 in Proposition 2 is satisfied. The critical value of flip bifurcation is τ 1 = 0.501806. Set τ = τ 1 , then the eigenvalues are -1 and 0.370478, and η 1 = -3.98561 < 0, η 2 = 0.989406 > 0. As stated by Theorem 3, the bifurcated period-2 orbit is stable when τ is in the right neighborhood of τ 1 . We plot the corresponding bifurcation diagram with τ ∈ [0.4, 0.7], please see Fig. 1(a). The initial value of flip bifurcation simulation is (P 2 + 0.0001, W 2 + 0.0015). From Fig. 1(a) we can clearly see the period-doubling cascade of the wrack biomass W . In order to get some detailed information from the period window, we plot a local amplification diagram in Fig. 1(b) with τ ∈ [0.656, 0.66]. In this period window, we see a period-6 orbit.
Corresponding to the flip bifurcation diagram, we draw the maximum Lyapunov exponents in Fig. 2, which can help us determine the chaotic and non-chaotic behavior quantitatively. From Fig. 2(b), we can see that the maximum Lyapunov exponent is above zero when τ is around 0.659. This means that the chaotic behavior may occur.
According to the flip bifurcation diagram and the maximum Lyapunov exponents, we show the dynamic evolution of map (10). We first illustrate the dynamics of the plant and wrack biomass changing over time for different given τ values. Then we plot the phase portraits to illustrate the dynamic transition of map (10).
Set τ = 0.45 < τ 1 , then the eigenvalues are -0.793523 and 0.435469. We choose an initial state (P 2 + 0.0001, W 2 -0.0015) and demonstrate the biomass of plant P and wrack W  Fig. 2(a) changing over time in Fig. 3(a). We can see that the fixed point (P 2 , W 2 ) = (0.7599, 0.7599) is a stable node. Set τ = 0.502, then the eigenvalues are -1.00078 and 0.370234. The fixed point (P 2 , W 2 ) is an unstable node. The loss of the stability of the node results in the emergence of stable period-2 points (0.7491, 0.7633) and (0.7704, 0.7562). Starting from the same initial state (P 2 + 0.0001, W 2 -0.0015), the biomass of the plant P and wrack W converge to the period-2 points. The biomass of the plant P finally oscillates between 0.7491 and 0.7704. Please see the red dots in Fig. 3(b). The biomass of the plant W finally oscillates between 0.7633 and 0.7562, please see the blue dots in Fig. 3(b). This implies that the bifurcated period-2 orbit is stable, and the biomass of the plant and wrack finally oscillates between the period-2 points: (0.7491, 0.7633) and (0.7704, 0.7562). Continue increasing the value of τ , we observe the period doubling behaviors. For example, let τ = 0.61, we can see period-4 points (0.4451, 0.8203), (0.9506, 0.5941), (0.5495, 0.7831), and (0.9352, 0.6406) as shown in Fig. 3(c). Correspondingly, the dynamics changing over time are shown in Fig. 3(d). We can see that the plant biomass P (denoted by red dots) Further suitably increasing the value of τ , the period will be further doubled. For example, set τ = 0.6252, we see a stable period-8 orbit, please see Fig. 4(a). The biomass of the plant and wrack (P, W ) will converge to the period-8 points, and finally oscillate between eight alternate states. Since the dynamics are similar to those shown in Fig. 3(b) and Fig. 3(d), we do not show the dynamics changing over time.
In the following, we concentrate on the dynamic transition of map (10) as the value of τ increases. In Fig. 4(a), we see a period-8 orbit. For τ = 0.6302 and τ = 0.632, we get another two different periodic orbits. Please see Fig. 4(b)-(c). Next, we set τ = 0.6575, then we arrive at a period-6 orbit corresponding to the period window in Fig. 1(b), please  Fig. 4(d). At last, we exhibit a chaotic attractor with τ = 0.659 in Fig. 4(e). The local amplification of Fig. 4(e) is shown in Fig. 4(f ) to present some detailed information about the chaotic attractor.
From Figs. 3 and 4, we can clearly see the dynamic variation on the route from a stable fixed point to chaos as shown in the flip bifurcation diagram with various τ values.
In the sequel, we first demonstrate the Turing bifurcation diagram to determine the corresponding critical value τ . And then we draw the pattern formation regions. According to different regions, we exhibit the transition of patterns.
Set d P = 0.3, d W = 0.6, δ = 10, we plot the Z mτ diagram, please see Fig. 5(a). With these given parametric values, we get the critical value for Turing bifurcation τ ≈ 0.5018062. Combining the flip bifurcation curve τ = τ 1 with Turing bifurcation curve τ = τ , we show the pattern formation regions in Fig. 5(b) as d W varying from 0 to 8. Three regions are obtained: homogeneous stationary state region, pure-Turing instability region, and flip-Turing instability region.
In the following, we simulate the spatial patterns of the re-scaled wrack biomass W , which are induced by the pure-Turing, flip-Turing instability and chaos mechanisms in Fig. 6 corresponding to the cases in Fig. 3 and Fig. 4. All the figures shown about the spatial patterns are taken as the spatial distribution of the CMLs model at t = 100,000. The initial state is a random perturbation to the homogeneous stationary state (P 2 , W 2 ).
Set d W = 0.6, for τ = 0.45, then neither Turing bifurcation nor flip bifurcation will occur, the stable homogeneous stationary state is locally uniformly stable. Therefore, no spatial patterns will be formed. The biomass W is projected onto the 150 × 150 rectangular grid, please see Fig. 6(a). For τ = 0.502, CMLs model (4)-(7) undergoes Turing bifurcation and flip bifurcation simultaneously. At this time the CMLs model will form spatial heterogeneous patterns, which are induced by the flip-Turing instability mechanism, please see Fig. 6(b). The spatial distribution of the wrack are mainly concentrated on two alternative states (represented by two different colors: yellow (the wrack biomass W = 0.7633) and blue (the wrack biomass W = 0.7562)), namely, the period-2 points. Let τ = 0.61, then we  In the meantime, the related patterns also show the characteristics of chaos. The pattern is mosaic type. We cannot tell how many colors are in the pattern. From the above simulations, we can see that the patterns gradually transit to fragment type, and finally chaos as τ varies from 0.45 to 0.659. In the transition process, the patterns also exhibit period-doubling phenomenon in space.

Neimark-Sacker bifurcation and Turing patterns
In this subsection, we first mainly give numerical illustrations about the results of Neimark-Sacker bifurcation obtained in Theorem 4. And then combined with the Turing bifurcation, we show the related patterns. Throughout this subsection, the parameters ξ and s are fixed as ξ = 5.58, s = 0.9. The computational grid is n × n = 200 × 200.
From the local amplification graphs of Fig. 7(b), as shown in Fig. 8, we find that when suitably choosing the τ value, such as τ = 2.151 and 2.21, the maximum Lyapunov exponent can be above zero, this implies that chaotic behavior may occur.
In the following, we first simulate the dynamic behavior of map (10) for the stable focus. Then we use the phase portraits to exhibit the dynamic transition from fixed point to chaotic behaviors associated with Neimark-Sacker bifurcation diagram.
Let τ = 1.95, then τ < τ * . We calculate the eigenvalues of the Jacobian matrix J(P 2 , W 2 ) and obtain λ 1,2 = -0.7191 ± 0.6832i and |λ 1,2 | = 0.9919 < 1. Therefore, the homogeneous stationary state (P 2 , W 2 ) = (0.1368, 0.1368) of map (10) is a stable focus. Choosing an initial state (P 2 + 0.02, W 2 + 0.002) (the green dot in Fig. 9), we obtain a trajectory with a discrete sequence of points (the blue dots in Fig. 9). For the sake of explaining the dynamic evolve- ments from the initial state to the focus, we connect these points in order with line segments and use arrows to indicate the direction from one state to the next state. Please see Fig. 9. Except for the initial state, we index the first 15 points with ST 1 ∼ ST 15 and assign two arrows for the first 15 points. One arrow tells the previous state, the other indicates the next state. The trajectory starting from the initial state proceeds to ST 1 , then to ST 2 , to ST 3 , to ST 4 , and so on. From the evolution of the first 15 points, we can conclude that the trajectory converges to the stable focus (P 2 , W 2 ) anticlockwise.
In the sequel, associated with the Neimark-Sacker bifurcation diagram Fig. 7(a), we use phase portraits to explicitly demonstrate the dynamic transition from the stable focus to chaotic behavior as τ increases. Please see Fig. 10. In Fig. 10(a), we show the stable focus with τ = 1.95. Increasing the value of τ to 1.9592, then τ > τ * , Theorem 4 and the Neimark-Sacker bifurcation diagram Fig. 7(a) imply that (P 2 , W 2 ) loses its stability to be an unstable focus. An attracting invariant circle is bifurcated, please see Fig. 10(b). Moreover, let τ = 2.055, 2.1, 2.12, we obtain some other invariant cycles with larger amplitudes as shown in Fig. 10(c), (e), and (f ). Continuing to increase the value of τ , such as τ = 2.07, τ = 2.14, and τ = 2.19, we arrive at the periodic windows corresponding to Fig. 7(a) respectively, where periodic-8, periodic-11, and period-6 orbits appear, please see Fig. 10(d) and  Fig. 11(a) and (c). According to the maximum Lyapunov exponents in Fig. 8, let τ = 2.151 and 2.21, we get two different chaotic attractors, please see Fig. 11(b) and (d). From Fig. 10 and Fig. 11, as τ increases, map (10) demonstrates a dynamic transition from the focus to Figure 11 Continuation of Fig. 10 invariant cycles, experiencing periodic window intermediately, and finally to chaotic attractors.
In the following, we illustrate the formation of spatial patterns on n × n = 200 × 200 lattices with space step δ = 10. The initial state is a random perturbation of the homogeneous stationary state (P 2 , W 2 ). All the figures shown about the spatial patterns are taken as the spatial distribution of the CMLs model at t = 150,000 -250,000.
In Fig. 12(a), we plot the Z m versus τ graph and get the critical value for Turing instability, namely τ = 1.95912137. In Fig. 12(b), we draw the Neimark-Sacker bifurcation curve and Turing bifurcation curve respectively. Then we obtain two pattern formation regions: pure-Turing instability region and Neimark-Sacker-Turing instability region. The right two regions are both none pattern formation regions, one is a homogeneous stationary state region, the other is a pure-Neimark-Sacker bifurcation region.
Set d P = 0.02, d W = 1.0, we demonstrate the spatial patterns self-organized under pure-Turing and Neimark-Sacker-Turing bifurcation and chaos.
When τ = 1.95, corresponding to Fig. 10(a), then τ < τ * and τ < τ . At this moment, neither Neimark-Saker bifurcation nor Turing bifurcation will occur. From Fig. 13(a), we see a uniform distribution at the homogeneous level of the wrack biomass. At τ = 1.959125, then τ < τ < τ * , the CMLs model undergoes pure-Turing instability. From Fig. 13(b), we can see spot patterns. Let τ = 1.9592, then τ > τ and τ > τ * , the CMLs model undergoes Neimark-Sacker-Turing instability. From Fig. 13(c), we can see that the spot pattern is deformed. The radius of the spot increases. Continue changing the value of τ , we can see other types of patterns, such as circles and spirals, as shown in Fig. 13(d)-(f ) and Fig. 14(a)-(b). Finally, when τ arrives at around 2.151, the patterns begin to be fragmented. The mosaic type patterns are presented. As τ increases, the patterns finally present complete disorder and chaos, please see Fig. 14(c)-(e). In the transition of spatial patterns, we find that the spot patterns break into stripe, then into spiral, and finally into mosaic type. The irregularity increases. Notice that at the periodic window the spatial patterns are different with τ = 2.14 (spiral) and τ = 2.19 (mosaic). For the two different chaotic attractors, the spatial patterns are both of mosaic type, but with larger value of τ , the spatial patterns are more fragmented. Moreover, during the simulation, we found that the patterns also have a period. The period may have something to do with the period of the invariant cycle bifurcated from the fixed point due to the Neimark-Sacker bifurcation.

Conclusions and discussions
In this article, we investigate the homogeneous and heterogeneous dynamics of a tussock sedge model with discrete time and space variables. Firstly, we develop our system by means of the coupled map lattice method. Secondly, by a series of theoretical analyses, we respectively give some precise conditions for the parameters to support the emergence of stable homogeneous stationary state, flip bifurcation, Neimark-Sacker bifurcation, and Turing instability. At last, we show some numerical examples to illustrate the theoretical results and present the abundant spatial patterns.
It is important to note that the dynamics of the CMLs model are more abundant. For the homogeneous behavior, the CMLs tussock sedge model can experience flip bifurcation and chaos, which cannot happen in the corresponding continuous system. Therefore, the homogeneous dynamics of the discrete form of the tussock sedge model are more abundant, especially for the chaotic behavior. It is worth noting that there are other bifurcation types for the fixed point, such as saddle-node bifurcation. This is our next research work. For the heterogeneous behavior, the Turing patterns can be generated by four mechanisms: pure Turing instability, flip-Turing instability, Neimark-Sacker-Turing instability, and chaos. In general, the patterns induced by flip-Turing instability present fragment doubling in space. When the fixed point is a focus, the patterns induced by pure-Turing instability can be spot type. Furthermore, the spot patterns show periodic deformation about time. If the CMLs model experiences Neimark-Sacker-Turing instability, the in- duced patterns have the characteristic of a spiral. In the transition of the spatial patterns, we observe that the change of τ value influences the sizes of spots in the patterns, please see Fig. 13(b)-(c). The patterns induced by chaos present disorder behavior both in time and space. We choose two different values of τ to see whether the patterns are the same. In fact, both patterns present irregular patterns. For larger τ , the patterns are more complex and fragmental.
The richness of the homogeneous dynamics of the discrete system, for example, the multiple period orbits, invariant cycles, and chaotic attractors, generated by flip bifurcation and Neimark-Sacker bifurcation, enhances the complexity and diversity of the patterns generated under the combination mechanisms. The various patterns imply that the discrete system can present more abundant spatiotemporal self-organization structures.