Analysis of Spatiotemporal Dynamic and Bifurcation in a Wetland Ecosystem

A wetland ecosystem is studied theoretically and numerically to reveal the rules of dynamics which can be quite accurate to better describe the observed spatial regularity of tussock vegetation. Mathematical theoretical works mainly investigate the stability of constant steady states, the existence of nonconstant steady states, and bifurcation, which can deduce a standard parameter control relation and in return can provide a theoretical basis for the numerical simulation. Numerical analysis indicates that the theoretical works are correct and thewetland ecosystemcan show rich dynamical behaviors not only regular spatial patterns.Our results further deepen and expand the study of dynamics in the wetland ecosystem. In addition, it is successful to display tussock formation in the wetland ecosystem may have important consequences for aquatic community structure, especially for species interactions and biodiversity. All these results are expected to be useful in the study of the dynamic complexity of wetland ecosystems.


Introduction
Ecosystems consist of organisms, the environment, and the interactions with each other [1].Many ecosystems exhibit remarkable complexity in both structure and dynamics, which has been used to predict the response of ecosystems to human interference [2,3].Complexity theory has successfully predicted that localized ecological interactions may strongly affect the organization of ecosystems on larger spatial scales, which may induce spatial self-organization.There is a growing body of literature emphasizing the importance of self-organization in determining the spatial complexity of ecosystems.Localized ecological interactions can generate striking large-scale spatial patterns in ecosystems through self-organization [4][5][6][7].In recent decades, regular spatial patterning, which is common in ecosystems, has been a dominant topic in the study of spatial self-organization [8,9].In addition, a number of theoretical ecologists have proposed various self-organization mechanisms to explain regular spatial patterns in real ecosystems [10][11][12][13][14].Many regular spatial patterns can be observed in various ecosystems, including arid ecosystems [12], savanna ecosystems [15], mussel beds [9], coral reefs [16], ribbon forests [17], intertidal mudflats [18], and marsh tussocks [8,19].In the study of regular spatial patterns, reaction-diffusion equations [20][21][22] have proved to be an effective tool and are extensively used.
Tussock formation is a global phenomenon that enhances microtopography and increases biodiversity by adding structure to ecological communities [19].Recently, many ecologists have paid increasing attention to experimental investigation of Carex stricta [8,19,23,24].Carex stricta, the tussock sedge, is a species with widespread distribution in fresh water marshes of North America.Lawrence and Zedler [19] studied Carex stricta tussock size in relation to elevation at a range of sites in southern Wisconsin, USA, and tested the effect of five hydroperiods and N+P addition on tussock formation during a three-year mesocosm experiment.Their work is significant to the study of tussock development in relation to environmental factors.van de Koppel and Crain [8] presented an experimental investigation of regular spatial patterning in Carex stricta.When van de Koppel and Crain studied a marsh tussocks ecosystem, they proposed a model to investigate alternative theoretical mechanisms driving regular spatial patterning in general.To obtain the best theoretical explanation of the observed spatial regularity of tussock vegetation, van de Koppel and Crain analyzed and compared these alternative mechanisms, using three simple, spatially explicit models to describe the growth of Carex stricta.These models, which describe the interaction between plants and wrack, have the following general structure [8]: where  is plant biomass,  is wrack biomass, () is a function describing the positive effect of plant biomass on its own growth,  is the specific rate of plant senescence, (, ) is a function describing the inhibiting effect of wrack on plant growth as a function of plant and wrack biomass,  is the decay rate of wrack,  1 and  2 are diffusion constants describing lateral movement of plants and wrack, and Δ is the usual Laplacian operator in two-dimensional space.van de Koppel and Crain [8] compared three alternative models to examine potential mechanisms driving regularity and found that scale-dependent inhibition provides the best explanation for the observed regular tussock spacing.van de Koppel and Crain emphasized the importance of scaledependent feedback in the formation of regular spatial patterns in ecosystems.Their work is very important in the study of the spatial distribution of Carex stricta.However, if   [8] was investigated in detail using analysis theory.Although the model is relatively simple, it is important and meaningful because it can reveal the theoretical mechanism that explains tussock formation in wetland ecosystems.The model can be expressed as follows: where /( + ) is added to the inhibition term to lower inhibition as  increases,  is the level of plant biomass where inhibition is reduced by half, and  is an inhibition coefficient [8].
Model (2) was analyzed under the following nonzero initial conditions: and the zero-flux boundary conditions: where  is the outward unit normal vector of the boundary Ω.
The rest of the paper is organized as follows.In Section 2, the properties of constant steady-state solutions of (2) are analyzed, including dissipation, persistence properties, and local and global stability.In Section 3, we give a priori estimates of the positive solutions and establish the existence and nonexistence of nonconstant positive steady-state solutions of this reaction-diffusion system.In Section 4, we perform bifurcation analysis, including Hopf bifurcation and Turing bifurcation.In Section 5, some numerical simulations are presented to verify the feasibility of the theoretical results.The paper ends with some conclusions and a brief discussion.

Local and Global Stability of Constant Steady States.
In this section, we will discuss the local and global stability of constant steady states.In the remaining part of this paper, we always assume that  < 1, which guarantees that system (2) has a unique positive constant steady state  * .

Nonconstant Steady States
In this section, we discuss the existence and nonexistence of nonconstant steady states of (2).The corresponding steadystate problem of (2) is the elliptic system: (33)

Existence of Nonconstant Steady States.
In this subsection, we discuss the existence of nonconstant positive steadystate solutions of (33), which guarantees the existence of stationary patterns [29][30][31][32], when the diffusion coefficients  1 ,  2 vary and the other parameters are kept fixed.From Theorem 3, we can easily conclude that  * is locally asymptotically stable if  * /( * + ) 2 − 1 ≤ 0. In this case, there might be no nonconstant positive steady-state solutions.Therefore, we will restrict our discussion to  * /( * + ) 2 − 1 > 0.
For simplicity, we denote Then,  1 ,  2 > 0 and (33) can be written as follows: Therefore,  * is the solution of (33) if and only if it satisfies where ( − Δ) −1 is the inverse of  − Δ with the homogeneous Neumann boundary condition.By simple computation, we can get It is known that  is an eigenvalue of   ( 1 ,  2 ; ) on   if and only if (1 +   ) is an eigenvalue of the following matrix: Obviously, tr Let Then ( 1 ,  2 ;   ) =  1  2 det(  ).If then ( 1 ,  2 ; ) = 0 has two real roots defined by , . ..}, and let (  ) be the multiplicity of   .
To calculate the index of ( 1 ,  2 ; ⋅) at  * , we introduce the following lemma.
for some  ≥ 1, and   = ∑  =1 (  ) is odd, then there exists a positive constant  * such that (33) has at least one nonconstant positive solution for all  2 ≥  * .Proof.Since  1 > 0, it follows that if  2 is large enough, then (60) holds and  + ( 1 ,  2 ) >  − ( 1 ,  2 ) > 0. Furthermore, lim If  1 / 1 ∈ (  ,  +1 ), there exists  0 ≫ 1 such that From Theorem 9, we know that there exists  >  0 such that (33) with  1 =  and  2 ≥  has no nonconstant positive solution.Let  > 0 be large enough such that  1 / 1 <  1 , then there exists  * >  such that In the following, we will prove that, for any  2 ≥  * , (33) has at least one nonconstant positive solution.By way of contradiction, assume that the assertion is not true for some  * 2 ≥  * .By using the homotopy invariance of the topological degree, we can derive a contradiction.

Hopf Bifurcation.
In this section, we study the existence of periodic solutions of system (2) by analyzing Hopf bifurcation from the positive constant steady state  * = ( * ,  * ), and we will show the existence of spatially homogeneous and nonhomogeneous periodic orbits.For convenience, in this section, we assume that all eigenvalues   ( ∈  0 ) are simple and that the space Ω is one-dimensional spatial domain (0, ), where  is a positive constant.In fact, this discussion also holds for the generic class of domains in higher dimensions.Let  =  * .In the following, we use  as the main bifurcation parameter.To look for a possible Hopf bifurcation value, we rewrite the following necessary and sufficient condition from [35,36].
For convenience, (2) can be translated into the following system by the translation of P =  −  * and Ŵ =  −  * , while still letting  and  denote P and Ŵ respectively.This yields Now, we consider the linearization near (0, 0) of the above system: The characteristic equation of   () is where   =  2 / 2 and (81) We will identify the Hopf bifurcation values  which satisfy condition (A1).
(a) When then Obviously, if (82) and (83) hold, then  0 (  0 ) = 0 and   (  0 ) < 0 for any  ≥ 1.And This corresponds to the Hopf bifurcation of spatially homogeneous periodic solution.Furthermore,   0 satisfying (82) is the unique case in which a Hopf bifurcation of spatially homogeneous periodic solution emerges for any  > 0. From the above discussion, we can easily get the following theorem.
The stability and bifurcation direction of the spatially homogeneous periodic orbits bifurcating from  =   0 can be determined by calculating the normal form according to [35].Because the calculation process is lengthy and complex, it will not be discussed here, and the interested reader may refer to [35] for details.Because the positive constant steady state  * = ( * ,  * ) is unstable, (  ,− ) has at least one pair of eigenvalues with positive real part.Therefore, the spatially nonhomogeneous periodic orbits bifurcating from  =   ,− (1 ≤  ≤ ) are all unstable.

Turing Bifurcation.
In this section, we will study Turing bifurcation of the positive constant steady state  * = ( * ,  * ) of system (2).Mathematically speaking, Turing instability occurs when a stable equilibrium point is driven unstable by the local dynamics and diffusion of species.In the following, we consider small space-and time-dependent perturbations for  * of system (2): where ,  are small enough and  is the wave number.Substituting (92) into (2), we can obtain its characteristic matrix: Correspondingly,  1 and  2 are the roots of the following characteristic equation: where The roots of (94) yield the dispersion relation: It is obvious that the stability conditions for  * of the corresponding ODE system are given by tr 0 =  11 −  < 0 and − 11 −  12 > 0. It is clear that tr  < tr 0 < 0. Therefore,  the stability of  * in system (2) changes with the sign of Δ  .It is easy to establish that Δ  < 0 for  2  1 <  2 <  2 2 , where If  2 1,2 has positive values, then the range of instability for  * can be obtained, which is called Turing space.In this case, the positive constant steady state  * of system (2) is unstable.Consequently, Turing patterns emerge.

Numerical Simulations
In this section, we present some numerical simulations to verify the feasibility of the theoretical results obtained in previous sections.In these numerical simulations, we consider one-dimensional spatial domain Ω = (0, 10).The reaction-diffusion system (2) contains six parameters: , , , ,  1 , and  2 .If we choose them satisfying the conditions of Theorem 3 (b), then  * is locally asymptotically stable, as shown in Figure 1.Therefore, the model develops to a homogeneous vegetation state despite strong heterogeneity introduced by the initial conditions (Figure 1).Moreover, if we choose them satisfying the conditions of Theorem 12, then system (2) undergoes a spatially homogeneous Hopf bifurcation and generates spatially homogeneous periodic solutions, as shown in Figure 2. In this case, the plants and wrack are distributed uniformly over the whole space while changing periodically over time (Figure 2).However, if the values of parameters ,  1 , and  2 are varied and the other parameters are fixed, from Theorem 13, system (2) undergoes a spatially nonhomogeneous Hopf bifurcation and generates spatially nonhomogeneous periodic solutions, as shown in Figure 3.Under the circumstances, the plants and wrack are heterogeneously distributed in space with high density on the edge and cyclical variations over time (Figure 3).In both cases, they change periodically with time.However, the former is spatially homogeneous regardless of location, while the latter is spatially nonhomogeneous and its distribution is strongly dependent on its spatial location.Furthermore, if the six parameters , , , ,  1 , and  2 are chosen suitably, then system (2) can generate nonconstant positive steadystate solutions, as shown in Figure 4.In this case, the plants and wrack will eventually become distributed periodically in space but time-independent (Figure 4).According to the above Turing bifurcation analysis, if  = 0.15,  = 0.1,  = 10,  = 0.2,  1 = 0.0001, and  2 = 0.01 as in [8], we can obtain Turing patterns [37].The numerical results about Turing patterns for model (2) can be found in [8].
In Figure 5, we use  and  2 as parameters to plot the bifurcation diagram of model (2).From Theorem 3, it is known that  1 = (0, 0) is always unstable.In region I,  * = ( * ,  * ) is globally asymptotically stable.All positive orbits converge to this constant steady state no matter where the initial point is.That is to say, the model will develop to a homogeneous state to any perturbations as time tends to infinity.In region II,  * is locally asymptotically stable.In this case, both the plants and wrack will be spatially homogeneously distributed to small spatiotemporal perturbations.Region III is Turing space described above, in which many Turing patterns exist.In this region, the plants and wrack have various spatial distributions.Regions IV and V are unstable regions.In region V, spatially nonhomogeneous Hopf bifurcation may occur.Under certain conditions, nonconstant positive steady-state solutions of system (2) may exist in regions III, IV, and V.

Conclusions and Discussion
In this paper, we have investigated a diffusive plant-wrack model with homogeneous Neumann boundary conditions theoretically and numerically.Mathematical analysis is used to study dissipation, persistence properties, and local and global stability of constant steady states.To prove global stability, a new method based on the upper and lower solution method is applied.Furthermore, the conditions of existence and nonexistence of nonconstant positive steadystate solutions of the reaction-diffusion system have been derived.Moreover, it has been found that system (2) has no nonconstant positive steady-state solutions if the diffusion coefficients  1 ,  2 are large enough.According to Hopf bifurcation analysis, the existence of periodic solutions of system (2), including spatially homogeneous and nonhomogeneous periodic solutions, has been investigated.Moreover, the conditions for Turing instability have also been obtained by Turing bifurcation analysis.
The results presented here indicate that system (2) can exhibit rich dynamics, not only regular spatial pattern formation as described in [8].From Figures 2 and 3, we can see clearly that the decay rate of wrack  and the diffusion constants  1 ,  2 have an important effect on the spatial distribution of plants and wrack.Different values of ,  1 , and  2 can lead to different spatial distributions of vegetation.Fixing , , and  as in Figures 2 and 3, if  < 0.15373 and  1 ,  2 meet appropriate conditions, instability and spatial oscillations may occur, which causes a nonhomogeneous spatial vegetation distribution.Many field studies have shown that oscillations of population density occur in reality [38,39], and Figures 2 and 3 indicate that the decay rate of wrack and diffusion can indeed result in oscillations with different amplitudes of biomass over time.From [8], it is known that a scale-dependent inhibitory effect of Carex stricta on its own growth, which is strongest at some distance from the tussock center, is key to explaining the formation of regular tussocks in space.In wetland ecosystems, the decay rate of wrack  reflects variations in mortality.Moreover, differences in diffusion can lead to changes of biomass in space, which further enhance the scale-dependent inhibitory effect.Therefore, variations of ,  1 , and  2 can create different spatial vegetation structures.Obviously, the amplitude of the plant biomass is larger than that of wrack in Figures 2 and 3 (note the difference in the scale of the -axis).Combining differences in amplitude with scale-dependent inhibition, it can be concluded that scale-dependent inhibition has more influence on plants than wrack.
By comparing Figures 1 and 4, it is apparent that different diffusion coefficients  1 and  2 can also lead to different distributions in space.In both cases, although the model ultimately develops to a homogeneous steady state, the spatial distributions of vegetation are not the same.The distribution is spatially homogeneous in Figure 1, but nonhomogeneous in Figure 4.That is to say, scale-dependent inhibition cannot affect the spatial distribution of vegetation in Figure 1 but has an important influence in Figure 4. From Figure 4, the same conclusion can be drawn as from Figures 2 and 3 that scaledependent inhibition has a greater impact on plants than wrack.
Using the model proposed by van de Koppel and Crain [8], we have presented the first systematic theoretical analysis of the wetland ecosystem, which further deepen and expand the study of dynamics in the wetland ecosystem.Many studies [8,[40][41][42][43] have suggested that tussock formation and the development of spatial self-organization structure within vegetation may have important consequences for community structure, especially for species interactions and biodiversity.Our results indeed coincide with these studies.All these results may provide a better understanding of the dynamics of wetland ecosystems.
1 and  2 vary, what are the dynamics?And what are the distributions of plants and wrack?Furthermore, if other parameters vary, what are the resulting changes in the dynamics and spatial distribution of vegetation?These problems are what we are interested in.One model proposed by van de Koppel and Crain 1 /.It follows from the comparison principle that (, ) ≤ () and lim  → ∞ sup max Ω (, ) ≤ (1 − )/.This completes the proof.