# Analysis of Spatiotemporal Dynamic and Bifurcation in a Wetland Ecosystem

Authors: Yi Wang; Min Zhao; Hengguo Yu; Chuanjun Dai; Debing Mei; Qi Wang; Zengling Ma
Journal: Discrete Dynamics in Nature and Society
 (2015)

---

## Abstract

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 the wetland ecosystem can 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.


---

## Body

## 1. 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–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–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–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 ofCarex 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] studiedCarex 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 inCarex 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 ofCarex stricta. These models, which describe the interaction between plants and wrack, have the following general structure [8]:(1)∂P∂t=P1-PFP-sP-GP,W+d1ΔP,∂W∂t=sP-bW+d2ΔW,where P is plant biomass, W is wrack biomass, F(P) is a function describing the positive effect of plant biomass on its own growth, s is the specific rate of plant senescence, G(P,W) is a function describing the inhibiting effect of wrack on plant growth as a function of plant and wrack biomass, b is the decay rate of wrack, d1 and d2 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 scale-dependent feedback in the formation of regular spatial patterns in ecosystems. Their work is very important in the study of the spatial distribution ofCarex stricta. However, if d1 and d2 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 [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:(2)∂P∂t=P1-P-sP-aPWkP+k+d1ΔP,∂W∂t=sP-bW+d2ΔW,where k/(P+k) is added to the inhibition term to lower inhibition as P increases, k is the level of plant biomass where inhibition is reduced by half, and a is an inhibition coefficient [8].Model (2) was analyzed under the following nonzero initial conditions:(3)Px,0=P0x>0,Wx,0=W0x>0,x∈Ω=R+2,and the zero-flux boundary conditions:(4)∂P∂n=∂W∂n=0,x∈∂Ω,where n is the outward unit normal vector of the boundary ∂Ω.The rest of the paper is organized as follows. In Section2, 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.


## 2. Constant Steady States

Obviously,E1=0,0 is the trivial constant steady state of system (2). If s<1, there exists a unique positive constant steady state E∗=P∗,W∗, and(5)P∗=b-ask-kb-sb+Δ2b,W∗=sbP∗,where Δ=-b+ask+kb+sb2-4(s-1)kb2. Furthermore, if s>1, there are no positive constant steady states of system (2).


### 2.1. Dissipation

Theorem 1.
Ifs<1, any solution P(x,t),W(x,t) of system (2) satisfies (6)limt→∞sup⁡maxΩ-⁡Px,t≤1-s,limt→∞sup⁡maxΩ-⁡Wx,t≤s1-sb.Therefore, for any ε>0, the rectangle 0,1-s+ε×0,s(1-s)/b+ε is a global attractor of (2) in R+2.Proof.
From the first equation of system (2),(7)∂P∂t-d1ΔP≤P1-s-P,and then from the comparison principle, limt→∞⁡sup⁡maxΩ-⁡Px,t≤1-s. This implies that, for any ε>0, there exists t1>0 such that P(x,t)≤c-1 for all x∈Ω- and t>t1, where c-1=1-s+ε>0. Then it follows that W satisfies:(8)∂W∂t-d2ΔW≤sc-1-bW,x∈Ω-,t>t1.Let z(t) be a solution of the following ordinary differential equation:(9)z′t=sc-1-bzt,t>t1,zt1=maxΩ-⁡zx,t1.Then limt→∞⁡z(t)=sc-1/b. It follows from the comparison principle that W(x,t)≤z(t) and limt→∞⁡sup⁡maxΩ-⁡W(x,t)≤s(1-s)/b. This completes the proof.


### 2.2. Persistence

Theorem 2.
If1-s-as(1-s)/b>0 and s<1, then system (2) is persistent.Proof.
From the proof of Theorem1, it can be concluded that there exists t2>t1 such that for any t>t2, W(x,t)≤c-2, where c-2=sc-1/b+ε>0. From the first equation of system (2), for all x∈Ω- and t>t2,(10)∂P∂t-d1ΔP=P1-P-sP-aPWkP+k≥P1-s-P-ac-2.Because 1-s-as(1-s)/b>0, then from the comparison principle, for any ε>0, there exists t3>t2 such that, for any t>t3,(11)Px,t≥c_1>0,where c_1=1-s-ac-2-ε.
ThenW satisfies(12)∂W∂t-d2ΔW≥sc_1-bW,x∈Ω-,t>t3.Let y(t) be a solution of the following ordinary differential equation:(13)y′t=sc_1-byt,t>t3,yt3=minΩ-⁡yx,t3.Then limt→∞⁡y(t)=sc_1/b. It follows from the comparison principle that there exists t4>t3 such that, for any t>t4,(14)Wx,t≥yt≥c_2>0,where c_2=sc_1/b-ε.
From (11) and (14), it can be easily obtained that limt→∞⁡inf⁡minΩ-⁡P(x,t)≥1-s-as(1-s)/b and limt→∞⁡inf⁡minΩ-⁡W(x,t)≥s1-s-as(1-s)/b/b. Therefore, together with Theorem 1, we can conclude that system (2) is persistent.


### 2.3. 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 thats<1, which guarantees that system (2) has a unique positive constant steady state E∗.Recall that-Δ under the homogeneous Neumann boundary condition has eigenvalues 0=μ0<μ1<μ2<⋯<∞ and limi→∞⁡μi=∞. Let S(μi) be the eigenspace corresponding to μi with multiplicity mi≥1. Let ϕij, 1≤j≤mi, be the normalized eigenfunctions corresponding to μi. Let Xij=c·ϕij:c∈R2, where ϕij:1≤j≤dim⁡S(μi) is an orthonormal basis of Sμi. Since X:=U=(u,v)∈C1Ω-2:∂u/∂n=∂v/∂n=0, then X=⊕i=1∞Xi and Xi=⊕j=1dimSμiXij.Theorem 3.
For system (2),(a)
the trivial constant steady stateE1=0,0 is unstable;(b)
ifaksP∗2<b(P∗+b)(P∗+k)2 and d1≥d2, the positive constant steady state E∗=P∗,W∗ is locally asymptotically stable.Proof.
The linearization of (2) at a constant solution U=P,W can be expressed by Ut=(DΔ+FU(U))U, where U=(P(x,t),W(x,t))T, D=diag⁡(d1,d2),(15)FU=P1-P-sP-aPWkP+k,sP-bWT,FUU=1-s-2P-ak2WP+k2-akPP+ks-b=J11J12J21J22.
For eachi≥0, Xi is invariant under the operator DΔ+FU(U), and λ is an eigenvalue of the operator on Xi if and only if it is an eigenvalue of the matrix(16)Ji=-μiD+FUU=-d1μi+J11J12J21-d2μi+J22.Therefore, the stability is reduced to consider the characteristic equation:(17)det⁡λI-Ji=λ2-trJi·λ+det⁡Ji,where tr(Ji)=-(d1+d2)μi+J11+J22, and det(Ji)=d1d2μi2-(d2J11+d1J22)μi+J11J22-J12J21.(a)
IfU∗=(0,0), then Ji=1-s-d1μi0s-b-d2μi. Obviously, if s<1, for i=0, one of the eigenvalues is 1-s>0. Therefore, E1=0,0 is unstable.(b)
IfU∗=P∗,W∗, then Ji=P∗aksP∗/b(P∗+k)2-1-d1μi-akP∗/P∗+ks-b-d2μi. It is easy to check that tr(Ji)=P∗aksP∗/b(P∗+k)2-1-b-(d1+d2)μi<0 and det⁡(Ji)=d1d2μi2+bd1-d2P∗aksP∗/b(P∗+k)2-1μi+ak2sP∗/(P∗+k)2+bP∗>0 for all i∈N0, if aksP∗2<b(P∗+b)(P∗+k)2 and d1≥d2. Therefore, E∗=P∗,W∗ is locally asymptotically stable.
Now, we give the result of global stability ofE∗=P∗,W∗ for system (2), which implies that plants and wrack will be spatially homogeneously distributed as time tends to infinity.Theorem 4.
Assume thataksP∗<b(P∗+k)2 and b>as; then E∗ is globally asymptotically stable for (2); that is, for any initial values P(x,0)>0, W(x,0)>0,(18)limt→∞⁡Px,t=P∗,limt→∞⁡Wx,t=W∗.Proof.
IfaksP∗<b(P∗+k)2, from Theorem 3, we can conclude that E∗ is locally asymptotically stable. Since b>as and s<1, then 1-s-as(1-s)/b>0. From Theorems 1 and 2, for any t>t4, (19)c_1≤Px,t≤c-1,c_2≤Wx,t≤c-2,and c_1, c_2, c-1, and c-2 satisfy(20)0≥1-c-1-s-akc_2c-1+k,0≥sc-1-bc-2,0≤1-c_1-s-akc-2c_1+k,0≤sc_1-bc_2.
Inequalities (20) show that (c-1,c-2) and (c_1,c_2) are a pair of coupled upper and lower solutions of system (2) as in the definition in [25, 26] because the nonlinearities in (2) are mixed quasimonotone. Therefore, it is clear that there exists L>0 such that, for any (c_1,c_2)≤(P1,W1),(P2,W2)≤(c-1,c-2),(21)P11-P1-sP1-aP1W1kP1+k-P21-P2+sP2+aP2W2kP2+k≤LP1-P2+W1-W2,sP1-bW1-sP2+bW2≤LP1-P2+W1-W2.We define two iteration sequences (c-1(n),c-2(n)) and (c_1(n),c_2(n)) as follows: for n≥1,(22)c-1n=c-1n-1+1Lc-1n-11-c-1n-1-s-akc_2n-1c-1n-1+k,c-2n=c-2n-1+1Lsc-1n-1-bc-2n-1,c_1n=c_1n-1+1Lc_1n-11-c_1n-1-s-akc-2n-1c_1n-1+k,c_2n=c_2n-1+1Lsc_1n-1-bc_2n-1,where (c-1(0),c-2(0))=(c-1,c-2) and (c_1(0),c_2(0))=(c_1,c_2). Then for n≥1,(23)c_1,c_2≤c_1n,c_2n≤c_1n+1,c_2n+1≤c-1n+1,c-2n+1≤c-1n,c-2n≤c-1,c-2,and there exist (c~1,c~2) and (c^1,c^2) such that(24)c_1,c_2≤c^1,c^2≤c~1,c~2≤c-1,c-2.Therefore, limn→∞⁡c_1n=c^1, limn→∞⁡c_2n=c^2, limn→∞⁡c-1n=c~1, limn→∞⁡c-2(n)=c~2, and(25)1-c~1-s-akc^2c~1+k=0,sc~1-bc~2=0,1-c^1-s-akc~2c^1+k=0,sc^1-bc^2=0.Simplifying (25) yields(26)1-s-c~1c~1+k=aksc^1b,1-s-c^1c^1+k=aksc~1b.Subtracting the first equation of (26) from the second equation, we have(27)c~1-c^11-s-k+aksb-c~1-c^1=0.If we assume that(28)c~1≠c^1,then1-s-k+aksb=c~1+c^1.Substituting (28) into (26),(29)1-s-c~1c~1+k=aks1-s-k+aks/b-c~1b,1-s-c^1c^1+k=aks1-s-k+aks/b-c^1b.Hence, the equation(30)1-s-YY+k=aks1-s-k+aks/b-Ybhas two positive roots c~1 and c^1. Equation (30) can be written as follows:(31)Y2+s+k-1-aksbY+kasb-1aksb+1-s=0.
Sinceb>as and s<1, (31) cannot have two positive roots. Hence, c~1=c^1, and, consequently, c~2=c^2. Then, from the results in [25, 26], the solution (P(x,t),W(x,t)) of system (2) satisfies(32)limt→∞⁡Px,t=P∗,limt→∞⁡Wx,t=W∗,uniformly forx∈Ω-.
Therefore, from the above analysis, it can be concluded that the constant steady stateE∗ is globally asymptotically stable for (2) if aksP∗<b(P∗+k)2 and b>as. This completes the proof.


## 2.1. Dissipation

Theorem 1.
Ifs<1, any solution P(x,t),W(x,t) of system (2) satisfies (6)limt→∞sup⁡maxΩ-⁡Px,t≤1-s,limt→∞sup⁡maxΩ-⁡Wx,t≤s1-sb.Therefore, for any ε>0, the rectangle 0,1-s+ε×0,s(1-s)/b+ε is a global attractor of (2) in R+2.Proof.
From the first equation of system (2),(7)∂P∂t-d1ΔP≤P1-s-P,and then from the comparison principle, limt→∞⁡sup⁡maxΩ-⁡Px,t≤1-s. This implies that, for any ε>0, there exists t1>0 such that P(x,t)≤c-1 for all x∈Ω- and t>t1, where c-1=1-s+ε>0. Then it follows that W satisfies:(8)∂W∂t-d2ΔW≤sc-1-bW,x∈Ω-,t>t1.Let z(t) be a solution of the following ordinary differential equation:(9)z′t=sc-1-bzt,t>t1,zt1=maxΩ-⁡zx,t1.Then limt→∞⁡z(t)=sc-1/b. It follows from the comparison principle that W(x,t)≤z(t) and limt→∞⁡sup⁡maxΩ-⁡W(x,t)≤s(1-s)/b. This completes the proof.


## 2.2. Persistence

Theorem 2.
If1-s-as(1-s)/b>0 and s<1, then system (2) is persistent.Proof.
From the proof of Theorem1, it can be concluded that there exists t2>t1 such that for any t>t2, W(x,t)≤c-2, where c-2=sc-1/b+ε>0. From the first equation of system (2), for all x∈Ω- and t>t2,(10)∂P∂t-d1ΔP=P1-P-sP-aPWkP+k≥P1-s-P-ac-2.Because 1-s-as(1-s)/b>0, then from the comparison principle, for any ε>0, there exists t3>t2 such that, for any t>t3,(11)Px,t≥c_1>0,where c_1=1-s-ac-2-ε.
ThenW satisfies(12)∂W∂t-d2ΔW≥sc_1-bW,x∈Ω-,t>t3.Let y(t) be a solution of the following ordinary differential equation:(13)y′t=sc_1-byt,t>t3,yt3=minΩ-⁡yx,t3.Then limt→∞⁡y(t)=sc_1/b. It follows from the comparison principle that there exists t4>t3 such that, for any t>t4,(14)Wx,t≥yt≥c_2>0,where c_2=sc_1/b-ε.
From (11) and (14), it can be easily obtained that limt→∞⁡inf⁡minΩ-⁡P(x,t)≥1-s-as(1-s)/b and limt→∞⁡inf⁡minΩ-⁡W(x,t)≥s1-s-as(1-s)/b/b. Therefore, together with Theorem 1, we can conclude that system (2) is persistent.


## 2.3. 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 thats<1, which guarantees that system (2) has a unique positive constant steady state E∗.Recall that-Δ under the homogeneous Neumann boundary condition has eigenvalues 0=μ0<μ1<μ2<⋯<∞ and limi→∞⁡μi=∞. Let S(μi) be the eigenspace corresponding to μi with multiplicity mi≥1. Let ϕij, 1≤j≤mi, be the normalized eigenfunctions corresponding to μi. Let Xij=c·ϕij:c∈R2, where ϕij:1≤j≤dim⁡S(μi) is an orthonormal basis of Sμi. Since X:=U=(u,v)∈C1Ω-2:∂u/∂n=∂v/∂n=0, then X=⊕i=1∞Xi and Xi=⊕j=1dimSμiXij.Theorem 3.
For system (2),(a)
the trivial constant steady stateE1=0,0 is unstable;(b)
ifaksP∗2<b(P∗+b)(P∗+k)2 and d1≥d2, the positive constant steady state E∗=P∗,W∗ is locally asymptotically stable.Proof.
The linearization of (2) at a constant solution U=P,W can be expressed by Ut=(DΔ+FU(U))U, where U=(P(x,t),W(x,t))T, D=diag⁡(d1,d2),(15)FU=P1-P-sP-aPWkP+k,sP-bWT,FUU=1-s-2P-ak2WP+k2-akPP+ks-b=J11J12J21J22.
For eachi≥0, Xi is invariant under the operator DΔ+FU(U), and λ is an eigenvalue of the operator on Xi if and only if it is an eigenvalue of the matrix(16)Ji=-μiD+FUU=-d1μi+J11J12J21-d2μi+J22.Therefore, the stability is reduced to consider the characteristic equation:(17)det⁡λI-Ji=λ2-trJi·λ+det⁡Ji,where tr(Ji)=-(d1+d2)μi+J11+J22, and det(Ji)=d1d2μi2-(d2J11+d1J22)μi+J11J22-J12J21.(a)
IfU∗=(0,0), then Ji=1-s-d1μi0s-b-d2μi. Obviously, if s<1, for i=0, one of the eigenvalues is 1-s>0. Therefore, E1=0,0 is unstable.(b)
IfU∗=P∗,W∗, then Ji=P∗aksP∗/b(P∗+k)2-1-d1μi-akP∗/P∗+ks-b-d2μi. It is easy to check that tr(Ji)=P∗aksP∗/b(P∗+k)2-1-b-(d1+d2)μi<0 and det⁡(Ji)=d1d2μi2+bd1-d2P∗aksP∗/b(P∗+k)2-1μi+ak2sP∗/(P∗+k)2+bP∗>0 for all i∈N0, if aksP∗2<b(P∗+b)(P∗+k)2 and d1≥d2. Therefore, E∗=P∗,W∗ is locally asymptotically stable.
Now, we give the result of global stability ofE∗=P∗,W∗ for system (2), which implies that plants and wrack will be spatially homogeneously distributed as time tends to infinity.Theorem 4.
Assume thataksP∗<b(P∗+k)2 and b>as; then E∗ is globally asymptotically stable for (2); that is, for any initial values P(x,0)>0, W(x,0)>0,(18)limt→∞⁡Px,t=P∗,limt→∞⁡Wx,t=W∗.Proof.
IfaksP∗<b(P∗+k)2, from Theorem 3, we can conclude that E∗ is locally asymptotically stable. Since b>as and s<1, then 1-s-as(1-s)/b>0. From Theorems 1 and 2, for any t>t4, (19)c_1≤Px,t≤c-1,c_2≤Wx,t≤c-2,and c_1, c_2, c-1, and c-2 satisfy(20)0≥1-c-1-s-akc_2c-1+k,0≥sc-1-bc-2,0≤1-c_1-s-akc-2c_1+k,0≤sc_1-bc_2.
Inequalities (20) show that (c-1,c-2) and (c_1,c_2) are a pair of coupled upper and lower solutions of system (2) as in the definition in [25, 26] because the nonlinearities in (2) are mixed quasimonotone. Therefore, it is clear that there exists L>0 such that, for any (c_1,c_2)≤(P1,W1),(P2,W2)≤(c-1,c-2),(21)P11-P1-sP1-aP1W1kP1+k-P21-P2+sP2+aP2W2kP2+k≤LP1-P2+W1-W2,sP1-bW1-sP2+bW2≤LP1-P2+W1-W2.We define two iteration sequences (c-1(n),c-2(n)) and (c_1(n),c_2(n)) as follows: for n≥1,(22)c-1n=c-1n-1+1Lc-1n-11-c-1n-1-s-akc_2n-1c-1n-1+k,c-2n=c-2n-1+1Lsc-1n-1-bc-2n-1,c_1n=c_1n-1+1Lc_1n-11-c_1n-1-s-akc-2n-1c_1n-1+k,c_2n=c_2n-1+1Lsc_1n-1-bc_2n-1,where (c-1(0),c-2(0))=(c-1,c-2) and (c_1(0),c_2(0))=(c_1,c_2). Then for n≥1,(23)c_1,c_2≤c_1n,c_2n≤c_1n+1,c_2n+1≤c-1n+1,c-2n+1≤c-1n,c-2n≤c-1,c-2,and there exist (c~1,c~2) and (c^1,c^2) such that(24)c_1,c_2≤c^1,c^2≤c~1,c~2≤c-1,c-2.Therefore, limn→∞⁡c_1n=c^1, limn→∞⁡c_2n=c^2, limn→∞⁡c-1n=c~1, limn→∞⁡c-2(n)=c~2, and(25)1-c~1-s-akc^2c~1+k=0,sc~1-bc~2=0,1-c^1-s-akc~2c^1+k=0,sc^1-bc^2=0.Simplifying (25) yields(26)1-s-c~1c~1+k=aksc^1b,1-s-c^1c^1+k=aksc~1b.Subtracting the first equation of (26) from the second equation, we have(27)c~1-c^11-s-k+aksb-c~1-c^1=0.If we assume that(28)c~1≠c^1,then1-s-k+aksb=c~1+c^1.Substituting (28) into (26),(29)1-s-c~1c~1+k=aks1-s-k+aks/b-c~1b,1-s-c^1c^1+k=aks1-s-k+aks/b-c^1b.Hence, the equation(30)1-s-YY+k=aks1-s-k+aks/b-Ybhas two positive roots c~1 and c^1. Equation (30) can be written as follows:(31)Y2+s+k-1-aksbY+kasb-1aksb+1-s=0.
Sinceb>as and s<1, (31) cannot have two positive roots. Hence, c~1=c^1, and, consequently, c~2=c^2. Then, from the results in [25, 26], the solution (P(x,t),W(x,t)) of system (2) satisfies(32)limt→∞⁡Px,t=P∗,limt→∞⁡Wx,t=W∗,uniformly forx∈Ω-.
Therefore, from the above analysis, it can be concluded that the constant steady stateE∗ is globally asymptotically stable for (2) if aksP∗<b(P∗+k)2 and b>as. This completes the proof.


## 3. Nonconstant Steady States

In this section, we discuss the existence and nonexistence of nonconstant steady states of (2). The corresponding steady-state problem of (2) is the elliptic system:(33)-d1ΔP=P1-P-sP-aPWkP+k,x∈Ω,t>0,-d2ΔW=sP-bW,x∈Ω,t>0,∂P∂n=∂W∂n=0,x∈∂Ω,t>0.


### 3.1. A Priori Estimates

Now, we first cite two important lemmas. For notational convenience, we denoteΛ=Λa,b,s,k.Lemma 5 (maximum principle [27]).
Suppose thatg∈CΩ-×R.(a)
Assume thatω∈C2(Ω)∩C1(Ω-) satisfies(34)Δωx+gx,ωx≥0,x∈Ω,∂ω∂n≤0,x∈∂Ω.If ωx0=maxΩ-⁡ω, then g(x0,ω(x0))≥0.(b)
Assume thatω∈C2(Ω)∩C1(Ω-) satisfies(35)Δωx+gx,ωx≤0,x∈Ω,∂ω∂n≥0,x∈∂Ω.If ωx0=minΩ-⁡ω, then g(x0,ω(x0))≤0.Lemma 6 (Harnack inequality [28]).
Letω∈C2Ω∩C1Ω- be a positive solution to Δω(x)+c(x)ω(x)=0, where c∈C(Ω-), satisfying the homogeneous Neumann boundary condition. Then there exists a positive constant C∗=C∗n,Ω,c∞, such that(36)maxΩ-⁡ω≤C∗minΩ-⁡ω.Theorem 7 (upper bounds).
For any positive solution(P(x),W(x)) of (33),(37)maxΩ-⁡Px≤1-s,maxΩ-⁡Wx≤s1-sb.Proof.
Let(38)Px1=maxΩ-⁡Px,Px2=minΩ-⁡Px,Wy1=maxΩ-⁡Wx,Wy2=minΩ-⁡Wx.Applying Lemma 5, we can get(39)Px11-Px1-sPx1-akPx1Wx1Px1+k≥0,sPy1-bWy1≥0.Then(40)Px1≤1-s,Wy1≤sPy1b.Therefore,(41)maxΩ-⁡Px≤1-s,maxΩ-⁡Wx≤s1-sb.Theorem 8 (lower bounds).
For any positive constantd, there exists a positive constant C_(Λ,Ω,d) such that every positive solution (P(x),W(x)) of (33) satisfies P(x)≥C_ and W(x)≥sC_/b if d1≥d.Proof.
Applying Lemma5 again, we can get(42)Px21-Px2-sPx2-akPx2Wx2Px2+k≤0,sPy2-bWy2≤0.Then(43)Px22+s+k-1Px2+akWx2-k1-s≥0,(44)Wy2≥sPy2b.
Note thatP(x2)≤P(x1) and W(x2)≤W(y1)≤sP(y1)/b≤sP(x1)/b from (40). Equation (43) implies that maxΩ-⁡P(x)=P(x1)≥C1 for some positive constant C1=C1(Λ).
Letc(x)=d1-1(1-s-P-akW/(P+k)). Then c(x)≤d-1(2+as(1-s)/b) if d1≥d. The Harnack inequality shows that there exists a positive constant C∗=C∗(Λ,Ω,d) such that(45)maxΩ-⁡Px≤C∗minΩ-⁡Px.Combining (45) with maxΩ-⁡P(x)≥C1, we can get minΩ-⁡P(x)≥C_ for some positive constant C_=C_(Λ,Ω,d). It follows from (44) that minΩ-⁡W(x)=W(y2)≥sP(y2)/b≥sP(x2)/b≥sC_/b. This completes the proof.


### 3.2. Nonexistence of Nonconstant Steady States

In this subsection, we show the nonexistence of nonconstant positive solutions of (33) by the effect of large diffusivity. For convenience, we denote(46)fP,W=P1-P-sP-aPWkP+k,gP,W=sP-bW.Theorem 9.
LetD1>(1-s)/μ1 be a fixed positive constant. Then there exists a positive constant D2(a,k,s,D1) such that (33) has no nonconstant positive solutions if d1≥D1 and d2≥D2.Proof.
Assume that(P(x),W(x)) is a positive solution of (33). Let (47)P-=Ω-1∫ΩPxdx,W-=Ω-1∫ΩWxdx;then(48)∫ΩP-P-dx=∫ΩW-W-dx=0.
Multiplying the first equation in (33) by P-P- and integrating over Ω, (49)d1∫Ω∇P-P-2dx=∫ΩfP,WP-P-dx=∫ΩfP,W-fP-,W-P-P-dx=∫Ω1-s-P+P-P-P-2dx-∫ΩakPWP+k-P-W-P-+kP-P-dx.Similarly, we multiply the second equation in (33) by W-W- to obtain(50)d2∫Ω∇W-W-2dx=∫ΩgP,WW-W-dx=∫ΩgP,W-gP-,W-W-W-dx=∫ΩsP-P-W-W-dx-∫ΩbW-W-2dx.Adding the above two equations, we can get(51)d1∫Ω∇P-P-2dx+d2∫Ω∇W-W-2dx=∫Ω1-s-P+P-+akPWP+kP-+k-akWP-+kP-P-2dx+∫Ωs-akP-P-+kP-P-·W-W-dx-b∫ΩW-W-2dx≤∫Ω1-s-P+P-P-P-2dx+∫Ωs+akP-P-W-W-dx-b∫ΩW-W-2dx≤∫Ω1-sP-P-2dx+s+ak2ε∫ΩP-P-2dx+s+ak2ε·∫ΩW-W-2dx-b∫ΩW-W-2dx=∫Ω1+s+ak2ε-sP-P-2dx+s+ak2ε-b∫ΩW-W-2dx.
From the Poincaré inequality, we can obtain(52)∫Ωd1μ1P-P-2+d2μ1W-W-2dx≤∫Ω1+s+ak2ε-sP-P-2+s+ak2ε-bW-W-2dx.
Sinced1μ1>1-s, there exists a sufficiently small ε0>0 such that d1μ1>1+s+ak/2ε0-s. Let D2=1/μ1(s+ak/2ε0-b); then we can conclude that P=P- and W=W- if d1≥D1 and d2≥D2. This completes the proof.


### 3.3. Existence of Nonconstant Steady States

In this subsection, we discuss the existence of nonconstant positive steady-state solutions of (33), which guarantees the existence of stationary patterns [29–32], when the diffusion coefficients d1, d2 vary and the other parameters are kept fixed. From Theorem 3, we can easily conclude that E∗ is locally asymptotically stable if aksP∗/b(P∗+k)2-1≤0. In this case, there might be no nonconstant positive steady-state solutions. Therefore, we will restrict our discussion to aksP∗/b(P∗+k)2-1>0.For simplicity, we denote(53)a1=P∗aksP∗bP∗+k2-1,a2=akP∗P∗+k,B=a1-a2s-b=FUU∗.Then, a1,a2>0 and (33) can be written as follows:(54)-ΔU=D-1FU,x∈Ω,∂U∂n=0,x∈∂Ω.Therefore, U∗ is the solution of (33) if and only if it satisfies(55)fd1,d2;U=U-I-Δ-1D-1FU+U=0,where (I-Δ)-1 is the inverse of I-Δ with the homogeneous Neumann boundary condition. By simple computation, we can get(56)DUfd1,d2;U=I-I-Δ-1D-1B+I.It is known thatλ is an eigenvalue of DUf(d1,d2;U) on Xi if and only if λ(1+μi) is an eigenvalue of the following matrix:(57)Mi=μiI-D-1B=μi-d1-1a1d1-1a2-d2-1sμi+d2-1b.Obviously,(58)det⁡Mi=d1-1d2-1d1d2μi2+d1b-d2a1μi+a2s-a1b,trMi=2μi+d2-1b-d1-1a1.Let(59)Hd1,d2;μ=d1d2μ2+d1b-d2a1μ+a2s-a1b.Then Hd1,d2;μi=d1d2det(Mi). If(60)d1b-d2a12>4d1d2a2s-a1b,then Hd1,d2;μ=0 has two real roots defined by(61)μ+d1,d2=d2a1-d1b+d1b-d2a12-4d1d2a2s-a1b2d1d2,μ-d1,d2=d2a1-d1b-d1b-d2a12-4d1d2a2s-a1b2d1d2.Denote Γ=Γd1,d2=μ:μ≥0,μ-(d1,d2)<μ<μ+(d1,d2), SP=μ0,μ1,μ2,μ3,…, and let m(μi) be the multiplicity of μi.To calculate the index off(d1,d2;·) at U∗, we introduce the following lemma.Lemma 10 (see [33]).
SupposeHd1,d2;μi≠0 for all μi∈SP, then (62)index⁡fd1,d2;·,U∗=-1σ,where(63)σ=∑μi∈Γ∩SPmμi,Γ∩SP≠Φ,0,Γ∩SP=Φ.In particular, if Hd1,d2;μ>0 for all μ≥0, then σ=0.
From Lemma10, in order to calculate the index of f(d1,d2;·) at U∗, we need to determine the range of μ for which Hd1,d2;μ<0.Theorem 11.
IfaksP∗/b(P∗+k)2-1>0, a1/d1∈(μq,μq+1) for some q≥1, and σq=∑i=1qm(μi) is odd, then there exists a positive constant d∗ such that (33) has at least one nonconstant positive solution for all d2≥d∗.Proof.
Sincea1>0, it follows that if d2 is large enough, then (60) holds and μ+d1,d2>μ-d1,d2>0. Furthermore, (64)limd2→+∞⁡μ+d1,d2=a1d1,limd2→+∞⁡μ-d1,d2=0.If a1/d1∈(μq,μq+1), there exists d0≫1 such that(65)μ+d1,d2∈μq,μq+1,0<μ-d1,d2<μ1,∀d2≥d0.From Theorem 9, we know that there exists d>d0 such that (33) with d1=d and d2≥d has no nonconstant positive solution. Let d>0 be large enough such that a1/d1<μ1, then there exists d∗>d such that(66)0<μ-d1,d2<μ+d1,d2<μ1,∀d2≥d∗.
In the following, we will prove that, for anyd2≥d∗, (33) has at least one nonconstant positive solution. By way of contradiction, assume that the assertion is not true for some d2∗≥d∗. By using the homotopy invariance of the topological degree, we can derive a contradiction.
Fixingd2=d2∗, for t∈0,1, we define(67)Dt=td1+1-td00td2+1-td∗,and consider the following problem:(68)-ΔU=D-1tFU,x∈Ω,∂U∂n=0,x∈∂Ω.Therefore, U is a nonconstant positive solution of (33) if and only if it is a solution of (68) with t=1. It is obvious that U∗=(P∗,W∗) is the unique positive constant solution of (68). For any t∈0,1, U is a nonconstant positive solution of (68) if and only if it is a solution of the following problem:(69)hU;t=U-I-Δ-1D-1tFU+U=0.
Our earlier analysis has shown that (69) has no nonconstant positive solution for t=0, and we have assumed that there is no such solution for t=1 at d2=d2∗. It is obvious that(70)hU;1=fd1,d2;U,hU;0=fd,d∗;U,DUfd1,d2;U∗=I-I-Δ-1D-1B+I,DUfd,d∗;U∗=I-I-Δ-1D~-1B+I,where D~=diag⁡d,d∗. From (65) and (66), we have (71)Γd1,d2∩SP=μ1,μ2,μ3,…,μq,Γd,d∗∩SP=Φ.Because σq is odd, Lemma 10 yields(72)index⁡h·;1,U∗=index⁡fd1,d2;·,U∗=-1σq=-1,index⁡h·;0,U∗=index⁡fd,d∗;·,U∗=-10=1.
From Theorems7 and 8, there exist positive constants C_ and C- such that the positive solutions of (69) satisfy C_<P(x),W(x)<C- on Ω- for all t∈0,1.
Set(73)Σ=U∈X:C_<ux,vx<C-.Then hU;t≠0 for all U∈∂Σ and t∈0,1. By the homotopy invariance of the Leray-Schauder degree [34],(74)deg⁡h·;0,Σ,0=deg⁡h·;1,Σ,0.Since both equations hU;0=0 and hU;1=0 have the unique positive solution U∗=(P∗,W∗) in Σ, then(75)deg⁡h·;0,Σ,0=index⁡h·;0,U∗=1,deg⁡h·;1,Σ,0=index⁡h·;1,U∗=-1.This contradicts (74), and the proof is complete.


## 3.1. A Priori Estimates

Now, we first cite two important lemmas. For notational convenience, we denoteΛ=Λa,b,s,k.Lemma 5 (maximum principle [27]).
Suppose thatg∈CΩ-×R.(a)
Assume thatω∈C2(Ω)∩C1(Ω-) satisfies(34)Δωx+gx,ωx≥0,x∈Ω,∂ω∂n≤0,x∈∂Ω.If ωx0=maxΩ-⁡ω, then g(x0,ω(x0))≥0.(b)
Assume thatω∈C2(Ω)∩C1(Ω-) satisfies(35)Δωx+gx,ωx≤0,x∈Ω,∂ω∂n≥0,x∈∂Ω.If ωx0=minΩ-⁡ω, then g(x0,ω(x0))≤0.Lemma 6 (Harnack inequality [28]).
Letω∈C2Ω∩C1Ω- be a positive solution to Δω(x)+c(x)ω(x)=0, where c∈C(Ω-), satisfying the homogeneous Neumann boundary condition. Then there exists a positive constant C∗=C∗n,Ω,c∞, such that(36)maxΩ-⁡ω≤C∗minΩ-⁡ω.Theorem 7 (upper bounds).
For any positive solution(P(x),W(x)) of (33),(37)maxΩ-⁡Px≤1-s,maxΩ-⁡Wx≤s1-sb.Proof.
Let(38)Px1=maxΩ-⁡Px,Px2=minΩ-⁡Px,Wy1=maxΩ-⁡Wx,Wy2=minΩ-⁡Wx.Applying Lemma 5, we can get(39)Px11-Px1-sPx1-akPx1Wx1Px1+k≥0,sPy1-bWy1≥0.Then(40)Px1≤1-s,Wy1≤sPy1b.Therefore,(41)maxΩ-⁡Px≤1-s,maxΩ-⁡Wx≤s1-sb.Theorem 8 (lower bounds).
For any positive constantd, there exists a positive constant C_(Λ,Ω,d) such that every positive solution (P(x),W(x)) of (33) satisfies P(x)≥C_ and W(x)≥sC_/b if d1≥d.Proof.
Applying Lemma5 again, we can get(42)Px21-Px2-sPx2-akPx2Wx2Px2+k≤0,sPy2-bWy2≤0.Then(43)Px22+s+k-1Px2+akWx2-k1-s≥0,(44)Wy2≥sPy2b.
Note thatP(x2)≤P(x1) and W(x2)≤W(y1)≤sP(y1)/b≤sP(x1)/b from (40). Equation (43) implies that maxΩ-⁡P(x)=P(x1)≥C1 for some positive constant C1=C1(Λ).
Letc(x)=d1-1(1-s-P-akW/(P+k)). Then c(x)≤d-1(2+as(1-s)/b) if d1≥d. The Harnack inequality shows that there exists a positive constant C∗=C∗(Λ,Ω,d) such that(45)maxΩ-⁡Px≤C∗minΩ-⁡Px.Combining (45) with maxΩ-⁡P(x)≥C1, we can get minΩ-⁡P(x)≥C_ for some positive constant C_=C_(Λ,Ω,d). It follows from (44) that minΩ-⁡W(x)=W(y2)≥sP(y2)/b≥sP(x2)/b≥sC_/b. This completes the proof.


## 3.2. Nonexistence of Nonconstant Steady States

In this subsection, we show the nonexistence of nonconstant positive solutions of (33) by the effect of large diffusivity. For convenience, we denote(46)fP,W=P1-P-sP-aPWkP+k,gP,W=sP-bW.Theorem 9.
LetD1>(1-s)/μ1 be a fixed positive constant. Then there exists a positive constant D2(a,k,s,D1) such that (33) has no nonconstant positive solutions if d1≥D1 and d2≥D2.Proof.
Assume that(P(x),W(x)) is a positive solution of (33). Let (47)P-=Ω-1∫ΩPxdx,W-=Ω-1∫ΩWxdx;then(48)∫ΩP-P-dx=∫ΩW-W-dx=0.
Multiplying the first equation in (33) by P-P- and integrating over Ω, (49)d1∫Ω∇P-P-2dx=∫ΩfP,WP-P-dx=∫ΩfP,W-fP-,W-P-P-dx=∫Ω1-s-P+P-P-P-2dx-∫ΩakPWP+k-P-W-P-+kP-P-dx.Similarly, we multiply the second equation in (33) by W-W- to obtain(50)d2∫Ω∇W-W-2dx=∫ΩgP,WW-W-dx=∫ΩgP,W-gP-,W-W-W-dx=∫ΩsP-P-W-W-dx-∫ΩbW-W-2dx.Adding the above two equations, we can get(51)d1∫Ω∇P-P-2dx+d2∫Ω∇W-W-2dx=∫Ω1-s-P+P-+akPWP+kP-+k-akWP-+kP-P-2dx+∫Ωs-akP-P-+kP-P-·W-W-dx-b∫ΩW-W-2dx≤∫Ω1-s-P+P-P-P-2dx+∫Ωs+akP-P-W-W-dx-b∫ΩW-W-2dx≤∫Ω1-sP-P-2dx+s+ak2ε∫ΩP-P-2dx+s+ak2ε·∫ΩW-W-2dx-b∫ΩW-W-2dx=∫Ω1+s+ak2ε-sP-P-2dx+s+ak2ε-b∫ΩW-W-2dx.
From the Poincaré inequality, we can obtain(52)∫Ωd1μ1P-P-2+d2μ1W-W-2dx≤∫Ω1+s+ak2ε-sP-P-2+s+ak2ε-bW-W-2dx.
Sinced1μ1>1-s, there exists a sufficiently small ε0>0 such that d1μ1>1+s+ak/2ε0-s. Let D2=1/μ1(s+ak/2ε0-b); then we can conclude that P=P- and W=W- if d1≥D1 and d2≥D2. This completes the proof.


## 3.3. Existence of Nonconstant Steady States

In this subsection, we discuss the existence of nonconstant positive steady-state solutions of (33), which guarantees the existence of stationary patterns [29–32], when the diffusion coefficients d1, d2 vary and the other parameters are kept fixed. From Theorem 3, we can easily conclude that E∗ is locally asymptotically stable if aksP∗/b(P∗+k)2-1≤0. In this case, there might be no nonconstant positive steady-state solutions. Therefore, we will restrict our discussion to aksP∗/b(P∗+k)2-1>0.For simplicity, we denote(53)a1=P∗aksP∗bP∗+k2-1,a2=akP∗P∗+k,B=a1-a2s-b=FUU∗.Then, a1,a2>0 and (33) can be written as follows:(54)-ΔU=D-1FU,x∈Ω,∂U∂n=0,x∈∂Ω.Therefore, U∗ is the solution of (33) if and only if it satisfies(55)fd1,d2;U=U-I-Δ-1D-1FU+U=0,where (I-Δ)-1 is the inverse of I-Δ with the homogeneous Neumann boundary condition. By simple computation, we can get(56)DUfd1,d2;U=I-I-Δ-1D-1B+I.It is known thatλ is an eigenvalue of DUf(d1,d2;U) on Xi if and only if λ(1+μi) is an eigenvalue of the following matrix:(57)Mi=μiI-D-1B=μi-d1-1a1d1-1a2-d2-1sμi+d2-1b.Obviously,(58)det⁡Mi=d1-1d2-1d1d2μi2+d1b-d2a1μi+a2s-a1b,trMi=2μi+d2-1b-d1-1a1.Let(59)Hd1,d2;μ=d1d2μ2+d1b-d2a1μ+a2s-a1b.Then Hd1,d2;μi=d1d2det(Mi). If(60)d1b-d2a12>4d1d2a2s-a1b,then Hd1,d2;μ=0 has two real roots defined by(61)μ+d1,d2=d2a1-d1b+d1b-d2a12-4d1d2a2s-a1b2d1d2,μ-d1,d2=d2a1-d1b-d1b-d2a12-4d1d2a2s-a1b2d1d2.Denote Γ=Γd1,d2=μ:μ≥0,μ-(d1,d2)<μ<μ+(d1,d2), SP=μ0,μ1,μ2,μ3,…, and let m(μi) be the multiplicity of μi.To calculate the index off(d1,d2;·) at U∗, we introduce the following lemma.Lemma 10 (see [33]).
SupposeHd1,d2;μi≠0 for all μi∈SP, then (62)index⁡fd1,d2;·,U∗=-1σ,where(63)σ=∑μi∈Γ∩SPmμi,Γ∩SP≠Φ,0,Γ∩SP=Φ.In particular, if Hd1,d2;μ>0 for all μ≥0, then σ=0.
From Lemma10, in order to calculate the index of f(d1,d2;·) at U∗, we need to determine the range of μ for which Hd1,d2;μ<0.Theorem 11.
IfaksP∗/b(P∗+k)2-1>0, a1/d1∈(μq,μq+1) for some q≥1, and σq=∑i=1qm(μi) is odd, then there exists a positive constant d∗ such that (33) has at least one nonconstant positive solution for all d2≥d∗.Proof.
Sincea1>0, it follows that if d2 is large enough, then (60) holds and μ+d1,d2>μ-d1,d2>0. Furthermore, (64)limd2→+∞⁡μ+d1,d2=a1d1,limd2→+∞⁡μ-d1,d2=0.If a1/d1∈(μq,μq+1), there exists d0≫1 such that(65)μ+d1,d2∈μq,μq+1,0<μ-d1,d2<μ1,∀d2≥d0.From Theorem 9, we know that there exists d>d0 such that (33) with d1=d and d2≥d has no nonconstant positive solution. Let d>0 be large enough such that a1/d1<μ1, then there exists d∗>d such that(66)0<μ-d1,d2<μ+d1,d2<μ1,∀d2≥d∗.
In the following, we will prove that, for anyd2≥d∗, (33) has at least one nonconstant positive solution. By way of contradiction, assume that the assertion is not true for some d2∗≥d∗. By using the homotopy invariance of the topological degree, we can derive a contradiction.
Fixingd2=d2∗, for t∈0,1, we define(67)Dt=td1+1-td00td2+1-td∗,and consider the following problem:(68)-ΔU=D-1tFU,x∈Ω,∂U∂n=0,x∈∂Ω.Therefore, U is a nonconstant positive solution of (33) if and only if it is a solution of (68) with t=1. It is obvious that U∗=(P∗,W∗) is the unique positive constant solution of (68). For any t∈0,1, U is a nonconstant positive solution of (68) if and only if it is a solution of the following problem:(69)hU;t=U-I-Δ-1D-1tFU+U=0.
Our earlier analysis has shown that (69) has no nonconstant positive solution for t=0, and we have assumed that there is no such solution for t=1 at d2=d2∗. It is obvious that(70)hU;1=fd1,d2;U,hU;0=fd,d∗;U,DUfd1,d2;U∗=I-I-Δ-1D-1B+I,DUfd,d∗;U∗=I-I-Δ-1D~-1B+I,where D~=diag⁡d,d∗. From (65) and (66), we have (71)Γd1,d2∩SP=μ1,μ2,μ3,…,μq,Γd,d∗∩SP=Φ.Because σq is odd, Lemma 10 yields(72)index⁡h·;1,U∗=index⁡fd1,d2;·,U∗=-1σq=-1,index⁡h·;0,U∗=index⁡fd,d∗;·,U∗=-10=1.
From Theorems7 and 8, there exist positive constants C_ and C- such that the positive solutions of (69) satisfy C_<P(x),W(x)<C- on Ω- for all t∈0,1.
Set(73)Σ=U∈X:C_<ux,vx<C-.Then hU;t≠0 for all U∈∂Σ and t∈0,1. By the homotopy invariance of the Leray-Schauder degree [34],(74)deg⁡h·;0,Σ,0=deg⁡h·;1,Σ,0.Since both equations hU;0=0 and hU;1=0 have the unique positive solution U∗=(P∗,W∗) in Σ, then(75)deg⁡h·;0,Σ,0=index⁡h·;0,U∗=1,deg⁡h·;1,Σ,0=index⁡h·;1,U∗=-1.This contradicts (74), and the proof is complete.


## 4. Bifurcation

### 4.1. 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 E∗=P∗,W∗, and we will show the existence of spatially homogeneous and nonhomogeneous periodic orbits. For convenience, in this section, we assume that all eigenvalues μi(i∈N0) are simple and that the space Ω is one-dimensional spatial domain (0,lπ), where l is a positive constant. In fact, this discussion also holds for the generic class of domains in higher dimensions. Let m=P∗. In the following, we use m as the main bifurcation parameter. To look for a possible Hopf bifurcation value, we rewrite the following necessary and sufficient condition from [35, 36].(A1)
There existsi∈N0 such that (76)TimH=0,DimH>0,TjmH≠0,DjmH≠0,∀j∈N0∖i.And for the unique pair of complex eigenvalues near the imaginary axis α(m)±iω(m), α′(mH)≠0.For convenience, (2) can be translated into the following system by the translation of P^=P-P∗ and W^=W-W∗, while still letting P and W denote P^ and W^ respectively. This yields(77)∂P∂t=P+P∗1-P∗-P-sP+P∗-akP+P∗W+W∗P+P∗+k+d1ΔP,x∈Ω,t>0,∂W∂t=sP+P∗-bW+W∗+d2ΔW,x∈Ω,t>0,∂P∂n=∂W∂n=0,x∈∂Ω,t>0.Now, we consider the linearization near (0,0) of the above system:(78)Lm=d1Δ+aksm2bm+k2-m-akmm+ksd2Δ-b,(79)Lim=-d1μi+aksm2bm+k2-m-akmm+ks-d2μi-b.The characteristic equation of Li(m) is (80)λ2-Timλ+Dim=0,i=0,1,2,…,where μi=i2/l2 and (81)Tim=-m+aksm2bm+k2-b-d1+d2μi,Dim=d1d2μi2+bd1-d2maksmbm+k2-1μi+mak2sm+k2+b.We will identify the Hopf bifurcation valuesm which satisfy condition (A1).(a) When(82)i=0,ifm0+bm0+k2=aksm02b,then(83)T0m0=0,D0m0=m0ak2sm0+k2+b>0,α′m0H=T0′m02=-m03-3km02+k22as/b-3m0-k32m0+k3≠0.Obviously, if (82) and (83) hold, then T0(m0H)=0 and Tj(m0H)<0 for any j≥1. And Dj(m0H)=d1d2μi2+b(d1-d2)μi+mak2s/(m+k)2+b>0 for any j≥1, if d1≥d2. This corresponds to the Hopf bifurcation of spatially homogeneous periodic solution. Furthermore, m0H satisfying (82) is the unique case in which a Hopf bifurcation of spatially homogeneous periodic solution emerges for any l>0. From the above discussion, we can easily get the following theorem.Theorem 12.
Assume thatd1≥d2. If (82) and (83) hold, then system (2) undergoes a spatially homogeneous Hopf bifurcation.In the following, we look for spatially nonhomogeneous Hopf bifurcation fori≥1. From Theorem 3, if aksP∗2<b(P∗+b)(P∗+k)2 and d1≥d2, then E∗=P∗,W∗ is locally asymptotically stable. In this case, there cannot be any Hopf bifurcation from E∗. Hence, any potential spatially nonhomogeneous Hopf bifurcation point mjH can make aksm2>b(m+b)(m+k)2 hold. In the following, we assume that aksm2>b(m+b)(m+k)2.(b) Wheni≥1, let f(m)=-m+aksm2/b(m+k)2-b. Then(84)f0=-b<0,limm→+∞⁡fm=-∞.Sinceaksm2>b(m+b)(m+k)2, there exist m0 and m^(m^>m0>0), such that f(m)>0(m∈(m0,m^)), where m0 is a spatially homogeneous Hopf bifurcation point.Obviously,(85)f′m=-m3-3km2+k22as/b-3m-k3m+k3,limm→+∞⁡f′m=-1,f′0<0.If(86)-1-s3-3k1-s2+k22asb-31-s-k3>0,then f′(1-s)>0, and there exist m_ and m-(0<m_<1-s<m-), such that f′(m_)=f′(m-)=0. Note that m=P∗∈0,1-s, and therefore f(m) attains its maximum M- at 1-s, and M-=f(1-s)>0. Define(87)ln=nd1+d2M-,n∈N,whereM-=f1-s>0.Then, for ln<l≤ln+1 and 1≤j≤n, we define mj,-H to be the roots of f(m)=(d1+d2)j2/l2. These points satisfy (88)0<m_<m0H<m1,-H<m2,-H<⋯<mn,-H<1-s<m-<m^.Clearly,Tj(mj,-H)=0 and Ti(mj,-H)≠0 for i≠j. It is obvious that α′(mj,-H)=Ti′(mj,-H)/2=f′(mj,-H)/2≠0. Now it is only necessary to consider whether Di(mj,-H)≠0 for all i≠j and Dj(mj,-H)>0. From (81), if d1/d2≥m/baksm/b(m+k)2-1, then Di(m)>0 for all i∈N0. Summarizing the above analysis, we can obtain the main results in this section.Theorem 13.
Assume thataksm2>b(m+b)(m+k)2, ln are defined as in (87). If (86) and d1/d2≥m/baksm/b(m+k)2-1 hold, then for any l∈ln,ln+1, there exist n points mj,-H(l)(1≤j≤n) satisfying(89)0<m_<m0H<m1,-H<m2,-H<⋯<mn,-H<1-s<m-<m^,such that system (2) undergoes a Hopf bifurcation at m=m0H or m=mj,-H, and the bifurcation periodic orbit near (m,P,W)=(mjH,P∗,W∗) can be parameterized in the form of m(θ)=mjH+o(θ) for θ∈(0,δ) and δ a small positive constant. Then (90)Pθx,t=P∗+θljcos⁡ωmjHtϕjx+oθ,Wθx,t=W∗mjH+θnjcosωmjHtϕjx+oθ,where ω(mjH)=Dj(mjH) is the corresponding time frequency, ϕj(x) is the corresponding spatial eigenfunction, and (lj,nj) is the corresponding eigenvector:(91)LmjH-iωmjHIlj,njTϕjx=0,0T,where L(m) is defined as in (78) and I is the identity map. Moreover,(1)
the bifurcation periodic orbits fromm=m0H are spatially homogeneous, which coincide with the periodic orbits of the corresponding ODE system;(2)
the bifurcation periodic orbits fromm=mj,-H are spatially nonhomogeneous.The stability and bifurcation direction of the spatially homogeneous periodic orbits bifurcating fromm=m0H 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 E∗=P∗,W∗ is unstable, L(mj,-H) has at least one pair of eigenvalues with positive real part. Therefore, the spatially nonhomogeneous periodic orbits bifurcating from m=mj,-H(1≤j≤n) are all unstable.


### 4.2. Turing Bifurcation

In this section, we will study Turing bifurcation of the positive constant steady stateE∗=P∗,W∗ 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 E∗ of system (2):(92)P=P∗+εexp⁡Kzi+λt,W=W∗+τexp⁡Kzi+λt,where ε, τ are small enough and K is the wave number. Substituting (92) into (2), we can obtain its characteristic matrix:(93)A=P∗akW∗P∗+k2-1-d1K2-akP∗P∗+ks-b-d2K2=A11-d1K2A12s-b-d2K2.Correspondingly, λ1 and λ2 are the roots of the following characteristic equation:(94)λ2-trKλ+ΔK=0,where(95)trK=A11-b-d1+d2K2,ΔK=d1d2K4+bd1-d2A11K2-bA11-sA12.The roots of (94) yield the dispersion relation:(96)λ1,2K=trK±trK2-4ΔK2.It is obvious that the stability conditions forE∗ of the corresponding ODE system are given by tr0=A11-b<0 and -bA11-sA12>0. It is clear that trK<tr0<0. Therefore, the stability of E∗ in system (2) changes with the sign of ΔK. It is easy to establish that ΔK<0 for K12<K2<K22, where (97)K1,22=d2A11-bd1∓bd1-d2A112+4d1d2bA11+sA122d1d2.IfK1,22 has positive values, then the range of instability for E∗ can be obtained, which is called Turing space. In this case, the positive constant steady state E∗ of system (2) is unstable. Consequently, Turing patterns emerge.


## 4.1. 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 E∗=P∗,W∗, and we will show the existence of spatially homogeneous and nonhomogeneous periodic orbits. For convenience, in this section, we assume that all eigenvalues μi(i∈N0) are simple and that the space Ω is one-dimensional spatial domain (0,lπ), where l is a positive constant. In fact, this discussion also holds for the generic class of domains in higher dimensions. Let m=P∗. In the following, we use m as the main bifurcation parameter. To look for a possible Hopf bifurcation value, we rewrite the following necessary and sufficient condition from [35, 36].(A1)
There existsi∈N0 such that (76)TimH=0,DimH>0,TjmH≠0,DjmH≠0,∀j∈N0∖i.And for the unique pair of complex eigenvalues near the imaginary axis α(m)±iω(m), α′(mH)≠0.For convenience, (2) can be translated into the following system by the translation of P^=P-P∗ and W^=W-W∗, while still letting P and W denote P^ and W^ respectively. This yields(77)∂P∂t=P+P∗1-P∗-P-sP+P∗-akP+P∗W+W∗P+P∗+k+d1ΔP,x∈Ω,t>0,∂W∂t=sP+P∗-bW+W∗+d2ΔW,x∈Ω,t>0,∂P∂n=∂W∂n=0,x∈∂Ω,t>0.Now, we consider the linearization near (0,0) of the above system:(78)Lm=d1Δ+aksm2bm+k2-m-akmm+ksd2Δ-b,(79)Lim=-d1μi+aksm2bm+k2-m-akmm+ks-d2μi-b.The characteristic equation of Li(m) is (80)λ2-Timλ+Dim=0,i=0,1,2,…,where μi=i2/l2 and (81)Tim=-m+aksm2bm+k2-b-d1+d2μi,Dim=d1d2μi2+bd1-d2maksmbm+k2-1μi+mak2sm+k2+b.We will identify the Hopf bifurcation valuesm which satisfy condition (A1).(a) When(82)i=0,ifm0+bm0+k2=aksm02b,then(83)T0m0=0,D0m0=m0ak2sm0+k2+b>0,α′m0H=T0′m02=-m03-3km02+k22as/b-3m0-k32m0+k3≠0.Obviously, if (82) and (83) hold, then T0(m0H)=0 and Tj(m0H)<0 for any j≥1. And Dj(m0H)=d1d2μi2+b(d1-d2)μi+mak2s/(m+k)2+b>0 for any j≥1, if d1≥d2. This corresponds to the Hopf bifurcation of spatially homogeneous periodic solution. Furthermore, m0H satisfying (82) is the unique case in which a Hopf bifurcation of spatially homogeneous periodic solution emerges for any l>0. From the above discussion, we can easily get the following theorem.Theorem 12.
Assume thatd1≥d2. If (82) and (83) hold, then system (2) undergoes a spatially homogeneous Hopf bifurcation.In the following, we look for spatially nonhomogeneous Hopf bifurcation fori≥1. From Theorem 3, if aksP∗2<b(P∗+b)(P∗+k)2 and d1≥d2, then E∗=P∗,W∗ is locally asymptotically stable. In this case, there cannot be any Hopf bifurcation from E∗. Hence, any potential spatially nonhomogeneous Hopf bifurcation point mjH can make aksm2>b(m+b)(m+k)2 hold. In the following, we assume that aksm2>b(m+b)(m+k)2.(b) Wheni≥1, let f(m)=-m+aksm2/b(m+k)2-b. Then(84)f0=-b<0,limm→+∞⁡fm=-∞.Sinceaksm2>b(m+b)(m+k)2, there exist m0 and m^(m^>m0>0), such that f(m)>0(m∈(m0,m^)), where m0 is a spatially homogeneous Hopf bifurcation point.Obviously,(85)f′m=-m3-3km2+k22as/b-3m-k3m+k3,limm→+∞⁡f′m=-1,f′0<0.If(86)-1-s3-3k1-s2+k22asb-31-s-k3>0,then f′(1-s)>0, and there exist m_ and m-(0<m_<1-s<m-), such that f′(m_)=f′(m-)=0. Note that m=P∗∈0,1-s, and therefore f(m) attains its maximum M- at 1-s, and M-=f(1-s)>0. Define(87)ln=nd1+d2M-,n∈N,whereM-=f1-s>0.Then, for ln<l≤ln+1 and 1≤j≤n, we define mj,-H to be the roots of f(m)=(d1+d2)j2/l2. These points satisfy (88)0<m_<m0H<m1,-H<m2,-H<⋯<mn,-H<1-s<m-<m^.Clearly,Tj(mj,-H)=0 and Ti(mj,-H)≠0 for i≠j. It is obvious that α′(mj,-H)=Ti′(mj,-H)/2=f′(mj,-H)/2≠0. Now it is only necessary to consider whether Di(mj,-H)≠0 for all i≠j and Dj(mj,-H)>0. From (81), if d1/d2≥m/baksm/b(m+k)2-1, then Di(m)>0 for all i∈N0. Summarizing the above analysis, we can obtain the main results in this section.Theorem 13.
Assume thataksm2>b(m+b)(m+k)2, ln are defined as in (87). If (86) and d1/d2≥m/baksm/b(m+k)2-1 hold, then for any l∈ln,ln+1, there exist n points mj,-H(l)(1≤j≤n) satisfying(89)0<m_<m0H<m1,-H<m2,-H<⋯<mn,-H<1-s<m-<m^,such that system (2) undergoes a Hopf bifurcation at m=m0H or m=mj,-H, and the bifurcation periodic orbit near (m,P,W)=(mjH,P∗,W∗) can be parameterized in the form of m(θ)=mjH+o(θ) for θ∈(0,δ) and δ a small positive constant. Then (90)Pθx,t=P∗+θljcos⁡ωmjHtϕjx+oθ,Wθx,t=W∗mjH+θnjcosωmjHtϕjx+oθ,where ω(mjH)=Dj(mjH) is the corresponding time frequency, ϕj(x) is the corresponding spatial eigenfunction, and (lj,nj) is the corresponding eigenvector:(91)LmjH-iωmjHIlj,njTϕjx=0,0T,where L(m) is defined as in (78) and I is the identity map. Moreover,(1)
the bifurcation periodic orbits fromm=m0H are spatially homogeneous, which coincide with the periodic orbits of the corresponding ODE system;(2)
the bifurcation periodic orbits fromm=mj,-H are spatially nonhomogeneous.The stability and bifurcation direction of the spatially homogeneous periodic orbits bifurcating fromm=m0H 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 E∗=P∗,W∗ is unstable, L(mj,-H) has at least one pair of eigenvalues with positive real part. Therefore, the spatially nonhomogeneous periodic orbits bifurcating from m=mj,-H(1≤j≤n) are all unstable.


## 4.2. Turing Bifurcation

In this section, we will study Turing bifurcation of the positive constant steady stateE∗=P∗,W∗ 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 E∗ of system (2):(92)P=P∗+εexp⁡Kzi+λt,W=W∗+τexp⁡Kzi+λt,where ε, τ are small enough and K is the wave number. Substituting (92) into (2), we can obtain its characteristic matrix:(93)A=P∗akW∗P∗+k2-1-d1K2-akP∗P∗+ks-b-d2K2=A11-d1K2A12s-b-d2K2.Correspondingly, λ1 and λ2 are the roots of the following characteristic equation:(94)λ2-trKλ+ΔK=0,where(95)trK=A11-b-d1+d2K2,ΔK=d1d2K4+bd1-d2A11K2-bA11-sA12.The roots of (94) yield the dispersion relation:(96)λ1,2K=trK±trK2-4ΔK2.It is obvious that the stability conditions forE∗ of the corresponding ODE system are given by tr0=A11-b<0 and -bA11-sA12>0. It is clear that trK<tr0<0. Therefore, the stability of E∗ in system (2) changes with the sign of ΔK. It is easy to establish that ΔK<0 for K12<K2<K22, where (97)K1,22=d2A11-bd1∓bd1-d2A112+4d1d2bA11+sA122d1d2.IfK1,22 has positive values, then the range of instability for E∗ can be obtained, which is called Turing space. In this case, the positive constant steady state E∗ of system (2) is unstable. Consequently, Turing patterns emerge.


## 5. 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: a, b, s, k, d1, and d2. If we choose them satisfying the conditions of Theorem 3 (b), then E∗ 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 b, d1, and d2 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 a, b, s, k, d1, and d2 are chosen suitably, then system (2) can generate nonconstant positive steady-state 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 s=0.15, b=0.1, a=10, k=0.2, d1=0.0001, and d2=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].Figure 1
Stable positive constant steady stateE∗=(0.03345633121,0.08364082803) of system (2) for s=0.25, b=0.1, a=10, k=0.2, d1=1, and d2=1 and initial conditions (P0,W0)=(0.03346+0.03sin⁡x,0.08364+0.03sin⁡x).
(a)
(b)Figure 2
Illustration of the existence of spatially homogeneous periodic solutions of system (2) for s=0.15, b=0.15373, a=10, k=0.2, d1=3, and d2=1 and initial conditions (P0,W0)=(0.12,0.1166).
(a)
(b)Figure 3
Illustration of the existence of spatially nonhomogeneous periodic solutions of system (2) for s=0.15, b=0.08, a=10, k=0.2, d1=1.05, and d2=0.7 and initial conditions (P0,W0)=(0.0539+0.01sin⁡x,0.101+0.01sin⁡x).
(a)
(b)Figure 4
Illustration of the existence of nonconstant positive steady-state solutions of system (2) for s=0.25, b=0.1, a=10, k=0.2, d1=0.06, and d2=10 and initial conditions (P0,W0)=(0.033456+0.03sin⁡5x,0.08364+0.06sin⁡5x).
(a)
(b)In Figure5, we use b and d2 as parameters to plot the bifurcation diagram of model (2). From Theorem 3, it is known that E1=(0,0) is always unstable. In region I, E∗=(P∗,W∗) 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, E∗ 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.Figure 5
Bifurcation diagram of system (2) for s=0.15, a=4, k=0.2, and d1=0.1. The green line is the Hopf bifurcation line, and the red curve is the Turing bifurcation curve.


## 6. 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 steady-state 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 d1, d2 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 b and the diffusion constants d1, d2 have an important effect on the spatial distribution of plants and wrack. Different values of b, d1, and d2 can lead to different spatial distributions of vegetation. Fixing s, a, and k as in Figures 2 and 3, if b<0.15373 and d1, d2 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 ofCarex 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 b 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 b, d1, and d2 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 P-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 Figures1 and 4, it is apparent that different diffusion coefficients d1 and d2 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 scale-dependent 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–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.



---

*Source: 185432-2015-10-04.xml*