Interaction between water and plants: Rich dynamics in a simple model

An ordinary differential equation model describing interaction of water and plants in ecosystem is proposed. Despite its simple looking, it is shown that the model possesses surprisingly rich dynamics including multiple stable equilibria, backward bifurcation of positive equilibria, supercritical or subcritical Hopf bifurcations, bubble loop of limit cycles, homoclinic bifurcation and Bogdanov-Takens bifurcation. We classify bifurcation diagrams of the system using the rain-fall rate as bifurcation parameter. In the transition from global stability of bare-soil state for low rain-fall to the global stability of high vegetation state for high rain-fall rate, oscillatory states or multiple equilibrium states can occur, which can be viewed as a new indicator of catastrophic environmental shift.


1.
Introduction. In the past four decades, the disappearance of vegetation which is referred as desertification has become a serious environmental problem, and the global ecosystems are threatened today more than ever before [7]. Biologically it indicates the ecological unbalance between the limited water resource and ecosystem engineers such as animals, plants or microorganisms [12,13]. How ecosystem engineers affect ecosystems is one of the main frontiers in ecology [25] and it has fascinated many ecologists [7,12,13]. The disappearance of vegetation may be a slow and gradual process in many forms, such as spotted vegetation patches, fairy rings, tiger bush bands, and bands, which have been reported and studied in different environment [15,14,22,23,25,28]. Mathematically it means formation of non-uniform spatiotemporal patterns, such as symmetry breaking, instabilities, and coexistence of stable states [7]. Pattern formation in water-plant systems in the semi-arid climatic zone has been investigated by many authors [17,31] with the applicability of Turing's idea [30]. The dynamic of this class of water-plant systems is governed by the competition of plants (shrubs, trees) for common resource (water, light) which admits spatial dynamics. The corresponding reaction-diffusion systems undergo a Turing-like stability change, which means a uniform steady state is stable with respect to the ordinary differential equation system but unstable with respect to the corresponding reaction-diffusion system, that may generate spotted vegetation patches, fairy rings, tiger bush bands, and bands in some parameter range. Understanding the underlying mechanism for generation of vegetation patterns and their observed resilience is considered as an important step toward a comprehension of the desertification process, where environmental effects such as climate changing and grazing destroy the natural balance toward stable aridity [28]. Vegetation pattern formation is a complex process involving not only the spatial physical movement of substance described by diffusion, advection and chemotaxis, but also the chemical interactions between resource and biomass. In this article, we propose a model of ordinary differential equations of water and biomass, based on a previous model proposed in Shnerb et al. [28] describing interaction of one plant species (shrubs or trees) and one resource. We show that the proposed two-dimensional simple ODE model possesses surprisingly rich dynamics including multiple stable equilibria, backward bifurcation of positive equilibria, supercritical or subcritical Hopf bifurcations, bubble loop of limit cycles, homoclinic bifurcations and saddlenode bifurcation. This model will serve as the kinetic model for our subsequent studies of spatiotemporal pattern formation for the corresponding partial differential equation model with effect of diffusion, advection and cross-diffusion.
Based on the model in [28], we propose a nondimensionalized water-biomass model in the following form Here w(t) is the amount of available water and b(t) is the amount of shrubs biomass; R is the rain-fall rate; the term −w represents water losses by percolation and evaporation; the term wb represents growth of shrubs as they consume water, while −λwb is the corresponding consumption of water by shrubs and λ is the water consumption rate in the presence of biomass; µ(b)b is the biomass death rate. A flow diagram for (1) is shown in Figure 1. In [28], a reaction-diffusion type system was proposed with the kinetic system (1) with µ(b) ≡ µ (a constant), and they showed that no Turing-like instability can occur for that system. In this study, we consider a biomass per capita death rate in the following form [20] (see Figure 2) where µ 0 and µ 1 are positive constants. Here µ(b) is assumed to be a decreasing function of biomass b in order to capture the infiltration feedback [6,7]. Biologically it can be explained as follows. The biomass can improve the soil environment locally which in turn increases infiltration at vegetation patch. Thus the larger biomass density results in the higher infiltration and the more soil water available and the smaller death rate of the biomass, whereas upon a constant death rate µ 0 . When µ 1 = 0 (that is, the biomass per capita death rate is a constant), system (1) is reduced to a well-known epidemic model (with some proper rescaling): SIR model with birth/death. In that context, w(t) is the susceptible population, b(t) is the infective population, and the removed population can be determined by the conservation of total population. The dynamics in that situation is well-known [3]. Indeed in this special case (1) has a bare-soil equilibrium (w 0 , b 0 ) = (R, 0) for all parameter values and a unique positive equilibrium (w + , b + ) = µ 0 , R − µ 0 λµ 0 if R > µ 0 . By using Lyapunov theory, we have the following results (detailed proof is given in subsection 2.1).
The bifurcation diagram of (1) with µ 1 = 0 looks like the one in Figure 3 (a) when using the rain-fall rate R as the bifurcation parameter. Hence one can conclude that the dynamics of (1) with µ 1 = 0 is a relatively simple one: when increasing 2974 XIAOLI WANG, JUNPING SHI AND GUOHONG ZHANG the rain-fall, there is a direct transition from the bare-soil state to the vegetation Figure 3. Possible bifurcation diagrams for (1)-(2) for different (λ, µ 0 , µ 1 ). In all diagrams, the horizontal axis is the rain-fall rate R, and the vertical axis is the biomass b. (a) forward bifurcation and no cycle; (b) forward bifurcation and bubble loop of cycles; (c) backward bifurcation and no cycle; (d) backward bifurcation and bubble loop of cycles; (e) backward bifurcation, and bubble loop of cycles with one subcritical and one supercritical Hopf bifurcations and a saddle-node bifurcation of cycles; (f) backward bifurcation, subcritical Hopf bifurcation, and branch of cycles into a homoclinic bifurcation; (g) backward bifurcation, supercritical Hopf bifurcation, and branch of cycles into a homoclinic bifurcation; (h) backward bifurcation, subcritical Hopf bifurcation, saddle-node bifurcation of cycles and branch of cycles into a homoclinic bifurcation.
In this paper, we consider (1) with a decreasing biomass per capita death rate as in (2) which captures the infiltration feedback, and we show that when µ 1 > 0, system (1) can exhibit surprisingly rich dynamics. While mathematical details of our result will appear in later sections, we summarize our main findings using the bifurcation diagrams in Figure 3. In all bifurcation diagrams shown in Figure 3, the horizontal axis is the rain-fall rate R, and the vertical axis is the biomass b. We use solid curve to represent stable equilibria, and dashed curve for unstable one; and we use an upper-lower pair of solid curve to represent stable periodic orbits, and dashed pair for unstable one. When µ 1 = 0, a simple diagram as in Figure 3 (a) 1. Backward bifurcation of positive equilibria. The positive equilibria emerges from the bare-soil ones through a transcritical bifurcation when increasing R. The ones in Figure 3 (a) and (b) are forward bifurcation, and the ones in (c)-(h) are backward ones. The phenomenon of backward bifurcations has been found in many epidemic models [2,3,9,34], and here it occurs as a result of higher infiltration for higher biomass density. 2. Saddle-node bifurcation of positive equilibria. If there is a backward bifurcation of positive equilibria, note that there is no positive equilibria for small positive R, then the branch of positive equilibria necessarily turns back at a saddlenode bifurcation point for equilibria as seen in Figure 3 (c)-(h). The two equilibria for R near the bifurcation point can be a saddle and a stable node (as in (c)-(e)), but they can also be a pair of a saddle and an unstable node (as in (f)-(h)). We remark that this was also found in some other water-plant models [14]. 3. Closed loop ("bubble" or "heart") of periodic orbits. For certain parameter range, oscillatory patterns emerge at the intermediate rain-fall rates. Indeed, we show that periodic orbits do not exist for small or large R. When there are two Hopf bifurcation points, the branch of periodic orbits form a closed loop starting and ending on the high vegetation equilibrium (Figure 3 (b), (d) and (e)). We call the ones in (b) and (d) "bubble", and the one in (e) "heart". The difference is that in (b) and (d), both Hopf bifurcations are supercritical and bifurcating limit cycles are stable, while in (e), one of Hopf bifurcations is subcritical. In the latter case, a saddle-node bifurcation of periodic orbits always occurs so the shape of the closed loop is like a heart. The bubble loop of limit cycles are more frequently found in delay differential equation models with delay as bifurcation parameter, for example [19,29], but not so often in purely ODE models. 4. Homoclinic orbit, homoclinic bifurcation and open ended branch ("lotus" or "pepper") of periodic orbits. In Figure 3 (f)-(h), an unstable node merges from the saddle-node bifurcation of equilibria, and there is only one Hopf bifurcation point when increasing R. Hence the branch of periodic orbits cannot be a closed loop with two ends. Instead the other end of the branch is a homoclinic orbit in this case. This is known as a homoclinic bifurcation for which the period of limit cycles tends to infinity when approaching to the homoclinic bifurcation point. When there is no more bifurcation on the branch of cycles, we call the open ended branch of periodic orbits to be a "lotus" ((f) and (g)), and we call it a "pepper" if there is another saddle-node bifurcation of cycles on the branch (h). Note that the "lotus" can either be forward (f) or backward (g). 5. Bogdanov-Takens bifurcation. All bifurcation phenomena described above are codimension one with parameter R. By using R and µ 1 as bifurcation unfolding parameters, we show that a codimension two Bogdanov-Takens bifurcation (see [16]) occurs at the intersection of the Hopf bifurcation curve and the equilibrium saddle-node bifurcation curve. In this bifurcation scenario, the saddle-node bifurcation produces a saddle and unstable node, the unstable undergoes a Hopf bifurcation generating a limit cycle, and the cycle continues in the "lotus" to a homoclinic orbit based on the saddle equilibrium. Compared with the simple monostable dynamics when µ 1 = 0, the system (1) with (2) shows bistability of (i) two equilibria (Figure 3 (c), (d)-(f) for some R); (ii) one equilibrium and one limit cycle (Figure 3 (d)-(h) for some R), or even possibly tristability of (iii) two equilibria and one limit cycle (Figure 3 (h) for some R). Specific conditions for achieving the bifurcation diagrams in Figure 3 and bistability/tristability are given in Sections 2-4. The system also have multiple locally stable equilibria, multiple periodic orbits and homoclinic orbit. Note that various bistability also appear in predator-prey type models with Allee effect in prey growth [33], while in (1) the "prey" (water) does not have an Allee effect type growth.
The bifurcation phenomena revealed above have been documented in textbooks of dynamical systems [16], but it is rare that all of them appear in the same simple model like (1)- (2). This is perhaps a minimal model to have all these dynamical behavior of backward equilibrium transcritical bifurcation, saddle-node bifurcations for equilibria and limit cycles, Hopf bifurcations, limit cycle bubble/heart, homoclinic bifurcation, and Bogdanov-Takens bifurcation.
Our mathematical results here have interesting biological implications. In all bifurcation scenarios in Figure 3, the bare-soil equilibrium is globally asymptotically stable for low rain-fall rate R, and the higher positive vegetation equilibrium is globally asymptotically stable for high rain-fall rate R (see subsection 2.3), which is not surprising. However for intermediate rain-fall rate, many different dynamical behavior are shown in Figure 3. Note that each of Figure 3 (c)-(h) has a backward bifurcation and a hysteresis loop when R decreases. The hysteresis loop predicts a catastrophic shift from the higher vegetation state to the bare-soil state at the saddle-node bifurcation point, which has been the main feature of many studies of ecological or physical phase transitions [21,25,27]. Moreover the early warning sign of catastrophic shift has been a recent hot topic, as it is important to predict and prevent such catastrophic changes [4,26]. Our results shown in Figure 3 (d)-(e) suggest that oscillatory patterns can occur at the intermediate rain-fall rate, that can be regarded as an early warning sign of getting close to the tipping point (saddle-node bifurcation point).
While the system (1) models the interaction between water and plants, it can also be interpreted as an epidemic model as mentioned earlier for the case of µ 1 = 0.
In that context, the term can be broken into two terms: µ 0 b represents moving from infected to removed class, and indicates the treatment of infected individuals. Similar SIR models with birth/death and treatment have been considered in [32,34,38,39]. In all these work, backward bifurcation has been discovered and some also found Hopf bifurcation and existence/nonexistence of limit cycles. But none of these work considered limit cycle bubble, homoclinic bifurcation, and Bogdanov-Takens bifurcation. This paper is organized as follows. In Section 2, we analyze the existence, local and global stability of equilibria. In Section 3, we show the occurrence of Hopf bifurcations and homoclinic bifurcations. The Bogdanov-Takens bifurcation analysis is given in Section 4. We end with some more discussions in Section 5.

Equilibria and stability.
2.1. Preliminaries. The system corresponding to (1) with death rate given in (2) is in form It is clear that F (w, b), G(w, b) are analytic functions in the first quadrant R 2 We also assume all the parameters R, λ, µ 0 , µ 1 are positive. We first show that system (3) is as "well behaved" as one intuits from the biological problem.
Lemma 2.1. Any solution of (3) is positive for t > 0, exists for t ∈ (0, ∞) and is eventually bounded, i.e., let then for any ε > 0, there exists T > 0, So the solution is positive. Next we show the boundedness of the solution. Firstly, from the water equation we see that lim sup t→∞ w(t) ≤ R. Then for any ε > 0, there exists a T 1 > 0 such that w(t) ≤ R + ε for t > T 1 . Let which yields that Then there exists a T 2 > 0 such that b(t) ≤ R + ε λk 1 for t > T 2 . Thus, we have This shows that there exists a T 3 > 0 such that w(t) ≥ Rk 1 R + k 1 + ε for t > T 3 .
Let T = max{T 1 , T 2 , T 3 }, we obtain (5) which also implies the global existence in time.

XIAOLI WANG, JUNPING SHI AND GUOHONG ZHANG
The system (3) has a trivial (bare-soil) equilibrium (w 0 , b 0 ) = (R, 0) for all parameters. A positive equilibrium (w, b) satisfies So b is a positive root of the quadratic equation Define First we have the following result about the positive equilibria of (3) from elementary calculation.
Theorem 2.2. Assume that the parameters R, λ, µ 0 , µ 1 are all positive and let ∆ 0 , R 0 , R 1 , b ± be defined as in (8) and (9). Then , b ± when R 1 < R < R 0 , and has no positive equilibrium when In the following (w + , b + ) is referred as the higher vegetation state, while (w − , b − ) is the lower vegetation state. It is useful to note that for any b > 0, there is at most (7), it is easy to solve that Note that in the (R, b) bifurcation diagram, R = R 1 is a saddle-node bifurcation point and R = R 0 is a transcritical bifurcation point where positive equilibria bifurcate from the bare-soil equilibrium. To determine whether the bifurcation at R = R 0 is supercritical or subcritical for the positive equilibria, from (10) we have the slope of the positive equilibrium branch at zero biomass density: Then the bifurcation at R = R 0 is subcritical (backward) if and it is supercritical (forward) if See Figure 3 (a)-(b) for illustration of forward bifurcations, and Figure 3 (c)-(h) for backward ones. We end this subsection with a proof of Proposition 1.1.
Proof of Proposition 1.1. For part 1, we use a Lyapunov function Its derivative along the solutions of (1) iṡ For part 2, we use a Lyapunov function Using the fact that R = λw + b + + w + and Then from LaSalle's invariance principle, the positive equilibrium (w + , b + ) is globally asymptotically stable for all nonnegative initial values except (R, 0) whenever it exists.
Remark 2.3. When µ 1 = 0, the model considered in Proposition 1.1 is the same as (after a rescaling) the classical SIR epidemic model with birth/death: (here a, α, β > 0) with basic reproduction number R 0 = β/α. In (15), S(t), I(t) and R(t) denote the number of individuals that are susceptible to infection, that are infectious, and that have been removed with immunity, respectively. The parameter a is the birth (death) rate, α is the incidence rate and β is the natural recovery rate of the infective individuals. Clearly (1) with µ 1 = 0 is equivalent to the first two equations of (15) with a "basic reproduction number" R 0 = R/µ 0 . We shall return to the discussion of connection between (1) (with µ 1 > 0) and epidemic models in Section 5.

2.2.
Local stability of equilibria. From Theorem 2.2, the system (3) may have three non-negative equilibria: (R, 0), (w − , b − ) and (w + , b + ). In this subsection we investigate the local stability of these equilibria. First we have the following results for the equilibria (R, 0) and (w − , b − ).
Proposition 2.4. Assume that the parameters R, λ, µ 0 , µ 1 are all positive. Then 1. The bare-soil equilibrium (R, 0) is locally asymptotically stable when R < R 0 , and it is unstable when Proof. The associated Jacobian matrix at an equilibrium (w, b) is given by Then at the bare-soil equilibrium (w 0 , b 0 ) = (R, 0), which shows that (R, 0) is locally asymptotically stable when R < R 0 , and it is unstable when R > R 0 . At the equilibrium (w − , b − ), the associated Jacobian matrix is On the other hand, by using (10), a direct calculation shows that is always a saddle equilibrium.
Next we discuss the local stability of (w + , b + ). We consider the two cases that λ satisfies 0 < λ < 1 and λ ≥ 1 separately. For the case of 0 < λ < 1, we define and Define the subregions I, II, III, IV and V in the µ 0 − µ 1 plane as follows (see
On the other hand, for the case that λ ≥ 1, we define the subregions VI, VII in µ 0 − µ 1 plane (see Figure 4 right panel): Then the local stability of (w + , b + ) is as follows: Theorem 2.6. Assume that the parameters R, λ, µ 0 , µ 1 are all positive and λ ≥ 1.
Let µ * 1 be defined as (19) and VI, VII be defined as in (23). Then we have The proofs of Theorems 2.5 and 2.6 are technical, and we postpone them to subsection 2.4.

2.3.
Global stability and nonexistence of periodic orbits. In this subsection, we discuss the global dynamical behavior of the system (3). First of all, the baresoil equilibrium is globally asymptotically stable if the rain-fall rate R is low. More precisely, we have Proposition 2.7. Assume that the parameters R, λ, µ 0 , µ 1 are all positive and let k 1 be defined as (4). If then (R, 0) is globally asymptotically stable for system (3).
Proof. According to Lemma 2.1, for any ε > 0, there exists a T 1 > 0 such that Thus there exists a T 2 such that b(t) < ε for t > T 2 , and This clearly shows that lim t→∞ w(t) = R.
Our next result is that when the rain-fall rate R is either low or high, the system (3) cannot have a periodic orbit. Define Then by the Dulac criterion we can eliminate the possibility of periodic orbits in the system (3) in the following result for low or high rain-fall rate R.
Proposition 2.8. Assume that the parameters R, λ, µ 0 , µ 1 are all positive and let R * and R * be defined as in (27). Then for 0 < R ≤ R * or R ≥ R * , system (3) has no periodic orbits.
Proof. We prove by using the Dulac Criterion. Define a Dulac function in the form (see [10,11]) where s(w) and r(b) are to be determined later. Then along solution orbits of system (3), we have Let s(w) = e cw where c is a constant to be determined. Then s (w)/s(w) = c and (28) is reduced to In (29), we choose r(b) = e λcb b c−1 so that r(b) satisfies From (29) and (30), it follows that Now we choose c = 1. Then which is negative if R ≤ R * = 1 + µ 0 . Thus from Dulac Criterion, there are no nontrivial periodic solutions if 0 < R ≤ R * and this completes the proof for the case of 0 < R ≤ R * .
Theorem 2.9 show that when µ 1 is small, then the dynamics of system (3) is dominated by equilibria. That is, any solution orbit tends to a non-negative equilibrium. When µ 1 is near 0 (satisfying (34) or (35)), the dynamics of (3) has a transition from the unique vegetation equilibrium to the bare-soil state as the rainfall rate R decreases across R 0 . The transition is a gradual one as the bifurcation is forward at R = R 0 . When µ 1 is in an intermediate range (satisfying (36)), there is a bistable parameter regime R 1 < R < R 0 , and there is a sudden transition from high vegetation equilibrium (w + , b + ) to the bar-soil one at R = R 1 . Hence a hysteresis effect occurs in this case. In all these parameter ranges for µ 1 , there is no limit cycle for (3), and the biomass always tends to an equilibrium level asymptotically. However as shown in Theorems 2.5 and 2.6, when µ 1 is larger, the system (3) can have oscillatory patterns for intermediate rain-fall rates. This will be shown in the Section 3.

2.4.
Proof of Theorems 2.5 and 2.6. In this subsection we prove Theorems 2.5 and 2.6. Recall the function in (6):

With a change of variable
The properties of functionR(B) can be summarized as follows: Let B + = b + + 1 then B + > 1 and we always haveR (B + ) > 0.
Proof of Theorem 2.5. Note that Det(J(w + , b + )) > 0. Then to consider the stability of (w + , b + ) we just need to consider the sign of the trace of J(w + , b + ), which has the same sign withf (B + ).
3. Hopf bifurcations and homoclinic bifurcations. From Theorem 2.5 and Theorem 2.6, there exist values of R, which denoted by R i , i = 2, 3 such that the corresponding characteristic matrix has a pair of complex roots, denoted by σ(R), From last section, we have wheref (B) is defined as (37). Here, we use the fact thatf (B 2 ) > 0 andf (B 3 ) < 0 which result in α (R 2 ) > 0 and α (R 3 ) < 0. Then we conclude that a Hopf bifurcation occurs at R = R i , i = 2, 3.
Theorem 3.1. Assume that the parameters R, λ, µ 0 , µ 1 are all positive and R = R i , i = 2, 3 are as defined in Theorems 2.5 and 2.6. Then system (3) undergo a Hopf bifurcation at R 2 and R 3 if exists.
We note that there are well-established methods (see [8,24]) to determine the direction of the Hopf bifurcations at R = R 2 and R = R 3 . Here we omit the detailed procedure as the computation is tedious and inconclusive. In the following we show by cases and examples that both supercritical and subcritical Hopf bifurcations can occur for different parameter ranges.
First when (µ 0 , µ 1 ) is in the parameter regions III, V, or VII, there are two Hopf bifurcation points R = R 2 and R = R 3 . In this case, a bounded branch of limit cycles emerges from the equilibrium (w + , b + ) at R = R 2 and returns to the same equilibrium at R = R 3 , so the branch is a loop or a "bubble". In Figure  5, the bifurcation diagrams and phase portraits of limit cycles for three sets of parameters in III, V, and VII respectively are shown. In each case, the bifurcation diagram (produced with MatCont) in the left panel shows two Hopf bifurcation points (labelled by "H"), and the numerically calculated first Lyapunov coefficients at both bifurcation points are negative in all three cases so the bifurcating periodic orbits are stable ones. Note that in Figure 5 (a)-(b), the dynamics is bistable with a stable limit cycle and a stable equilibrium (R, 0) when R 2 < R < R 3 , while in (c)-(d) and (e)-(f), the limit cycle is globally asymptotically stable when R 2 < R < R 3 . In all cases shown in Figure 5, the limit cycle appears to be unique, though we do 11.12 11.13 11.14 11.15 11.16 0 0.5   not have an analytic proof of that. For the parameter values in Figure 5 (a)-(b), the period of the (unique) limit cycle is decreasing in R as shown in 6 (c). In Figure 6 (a)-(b), the bifurcation diagram and phase portrait of limit cycles for a different value (µ 0 , µ 1 ) ∈ III is shown. Here the Hopf bifurcation at R = R 2 is supercritical but the one at R = R 3 is subcritical. The branch of limit cycles is still a closed loop but we call it a "heart" to indicate the existence of multiple periodic orbits for some values of R. Indeed here a saddle-node bifurcation for the periodic orbits occurs at some R = R 4 > R 3 , and there exist two periodic orbits for R 3 < R < R 4 (see Figure 6 (b) for phase portrait and Figure 6 (e)-(f) for the time series of the two periodic orbits). The dependence of period of periodic orbits on R is shown in 6 (d): for R close to R 2 , the period appears to be not monotone, and for R close to R 3 , the period is multi-valued as there are multiple periodic orbits.
Secondly, if (µ 0 , µ 1 ) is in the parameter region II, then there is only one Hopf bifurcation point R = R 3 . Here the branch of limit cycles is not a closed loop as in the "bubble" or the "heart" case above, but is an "open-ended" one. However the results in Proposition 2.8 imply that there is no periodic orbits for large R or small R, hence the branch of periodic orbits is bounded in R direction. On the other hand, Lemma 2.1 implies that all periodic orbits are bounded on the phase plane for bounded R values. Then from the global Hopf bifurcation theorem [1,5], the period of the periodic orbits on the bifurcating branch must be unbounded. We call such a bifurcation diagram a "lotus" with the Hopf bifurcation point as the base and the homoclinic bifurcation point the top of the lotus. In Figure 7 (a)-(b) and (c)-(d), the bifurcation diagrams and phase portraits of limit cycles for two different (µ 0 , µ 1 ) ∈ II are shown. For (µ 0 , µ 1 ) ∈ II, there is always a saddle-node bifurcation point at R = R 1 , which is labeled by "LP" in Figure 7 (a) and (c). The Hopf bifurcation at R = R 3 can be supercritical as in (a), or subcritical as in (c). In the former case, the bifurcating periodic orbits are stable while in the latter case they are unstable. In Figure 7 (e) and (f), one can observe that the period of cycles is increasing in R for the bifurcation diagram corresponding to (a), and the period is decreasing in R for the one corresponding to (c). In both cases, the period approach to ∞ when R tends to some R = R HL , which is a homoclinic bifurcation point where the limit cycles converge to a homoclinic orbit based on the saddle equilibrium (w − , b − ). This can be observed in Figure 7 (a) and (c) as the limit of minimum of the cycles is the lower positive equilibrium. Note that the sequence of Hopf bifurcation and homoclinic bifurcation is also the signature of the codim-2 Bogdanov-Takens bifurcation, which will be discussed later. Figure 8 shows the evolution of phase portraits of system (3) with λ = 0.2 and (µ 0 , µ 1 ) = (5, 5.2) ∈ II (which is the one in Figure 7 (a)) when R increases. One can easily calculate that the bifurcation values R 0 = 10.2 (backward transcritical bifurcation point), R 1 = 9.119216 (saddle-node bifurcation point), and R 3 = 9.161059 (Hopf bifurcation point). Here the Hopf bifurcation is subcritical so the bifurcating periodic orbits are unstable. Moreover the homoclinic bifurcation point is R HL = 9.1830408. Apparently when 0 < R < R 1 , the bare-soil equilibrium (R, 0) is the attractor; when R 1 < R < R 3 , both positive equilibria (w + , b + ) (unstable node) and (w − , b − ) (saddle) are unstable, and all positive orbits except the positive equilibria and the stable manifold (the green orbits) of (w − , b − ) tend to the trivial bare soil equilibrium (w 0 , b 0 ) = (R, 0) (see Figure 8 (a)); when R increases to R 3 < R < R HL , the equilibrium (w + , b + ) becomes stable and an unstable limit cycle emerges from the subcritical Hopf bifurcation (see Figure 8 (b)). In this case, 9.12 9.14 9.16 9. the unstable cycle serves as the separatrix between the basins of attraction of the equilibria (w + , b + ) and (R, 0). At R = R HL , a homoclinic orbit exists and it is the separatrix between the basins of attraction of the equilibria (w + , b + ) and (R, 0) (see Figure 8 (c)). When R increases to R HL < R < R 0 , the basins of attraction for the two locally stable equilibria are separated by the stable manifold of the saddle equilibrium (w − , b − ) (see Figure 8 (d)). In each of Figure 8 (b), (c) and (d), system (3) has a bistability of the bare-soil equilibrium (R, 0) and the positive equilibrium (w + , b + ), but the separatrix is different. When R > R 0 , the system is monostable as all solutions tend to the positive equilibrium (w + , b + ).  On the other hand, Figure 9 shows the evolution of phase portraits of system (3) with λ = 0.2 and (µ 0 , µ 1 ) = (7, 4.9) ∈ II (which is the one in Figure 7 (c)) when R increases. Now the bifurcation values are R 0 = 11.9, R 1 = 11.265296, and R 3 = 11.391151 and the Hopf bifurcation at R = R 3 is supercritical. Furthermore the homoclinic bifurcation point is R HL = 11.363757. When R > R 3 , the equilibrium (w + , b + ) is stable in Figure 9 (d). Then system (3) admits the bistability of the bare-soil equilibrium and a positive equilibrium when R 3 < R < R 0 and the region of attraction for these equilibria are separated by the separatrix (the green orbits) of saddle (w − , b − ). As R decreases to R HL < R < R 3 , the equilibrium (w + , b + ) becomes unstable and a stable limit cycle emerges from a supercritical Hopf bifurcation in Figure 9 (c). In this case, system (3) admits the bistability of the bare-soil equilibrium and a stable limit cycle with the stable manifolds (the green orbits) of the saddle (w − , b − ) as the separatrix of attraction. When R decreases to R HL , the period of the limit cycle tends to infinity and a homoclinic orbit emerges in Figure  9 (b). As R decreases to R 1 < R < R HL , the homoclinic orbit is broken and the equilibrium (w + , b + ) is still unstable. All of the positive orbits except for positive equilibria (w ± , b ± ) and stable manifold (the green orbits) of (w − , b − ) tend to the trivial bare-soil equilibrium (w 0 , b 0 ) = (R, 0) in Figure 9 (a).  Finally we mention that homoclinic bifurcation and multiple periodic orbits can occur for the same (µ 0 , µ 1 ) value. In Figure 10, the bifurcation diagram and phase portraits of limit cycles are shown for λ = 0.2 and (µ 0 , µ 1 ) = (7, 5) ∈ II. Here a subcritical Hopf bifurcation occurs at R = R 3 . In some range of R only one stable limit cycle exists (Figure 10 (b)), and the multi-valued period of the limit cycles in Figure 10 (c) implies that two periodic orbits can coexist (Figure 10 (d)). The period of the limit cycles tend to infinity at a critical value R HL = 11.47377 which is the homoclinic bifurcation point. Time series of the two limit cycles for the same parameter values are shown in Figure 10 (e) and (f). We call this bifurcation diagram a "pepper" as it resembles the shape of bell pepper.  Figure 11 shows the evolution of phase portraits when λ = 0.2 and (µ 0 , µ 1 ) = (7, 5) ∈ II (the one in Figure 10). The bifurcation points are R 0 = 12, R 1 = 11.332864, and R 3 = 11.497537 and the Hopf bifurcation is subcritical. Two additional bifurcation points which cannot be analytically solved are the limit cycle saddle-node bifurcation point R LP C = 11.49838 and the homoclinic bifurcation point R HL = 11.47377. Some parts of evolution are same as the ones in Figure  9, but in Figure 11 (d) there are indeed three locally stable states: the positive equilibrium (w + , b + ), the large amplitude limit cycle, and the bare-soil equilibrium (R, 0). This is the only "tristable" dynamics for (3), while most other cases for (3) are bistable or monostable.
Let u = w, v = b − c 3 wb and rewrite u, v as w, b respectively. Then system (43) becomes where P 3 (w, b), Q 3 (w, b) are smooth functions in (w, b) at least of order three. Fi- where P 4 (w, b), Q 4 (w, b) are smooth functions in (w, b) at least of order three. Note that and if µ 0 = µ * 0 . Thus, it follows from [16,Lemma 8.6] that the equilibrium (w * , b * ) of (2.1) is a cusp of codimension 2, i.e. it is a Bogdanov-Takens singularity if µ 0 = µ * 0 . Theorem 4.1 implies that system (3) may admit a Bogdanov-Takens bifurcation. From [16,Theorem 8.4] and symbolic calculation in Maple, we can find the versal unfolding and approximating bifurcation curves depending on the original parameters for (2.1) near the Bogdanov-Takens bifurcation point (see Appendix for details). Figure 12 illustrates the Bogdanov-Takens bifurcation in the R − µ 1 plane for several (λ, µ 0 ) values. In each panel of Figure 12, the cyan curve Γ 1 represents the Hopf bifurcation curve, the blue curve Γ 2 represents the saddle-node bifurcation curve (which only exists when 0 < λ < 1) and the black line Γ 3 is the transcritical bifurcation curve R = R 0 = µ 1 +µ 0 . When 0 < λ < 1, the saddle-node curve Γ 1 and the Hopf curve Γ 2 divide the R−µ 1 plane into three parts, marked by "nonexistence" (no positive equilibrium), "stable" ((w + , b + ) is stable) and "unstable" ((w + , b + ) is unstable), respectively (see Figure 12  parameter regions with 1 or 2 positive equilibria. The intersection point of the saddle-node curve Γ 1 and the Hopf curve Γ 2 is the Bogdanov-Takens bifurcation point marked by "BT", and the intersection of Γ 2 and Γ 3 is the cusp bifurcation point marked by "CP" (again by using package MatCont). Another degenerate bifurcation point marked by "GH" is the generalized Hopf bifurcation point where the first Lyapunov coefficient vanishes while the second Lyapunov coefficient does not vanish, and that is where the Hopf bifurcation changes from subcritical to supercritical. In Figure 12 (a), the Hopf bifurcation curve is a monotone curve, which implies that there is only one Hopf bifurcation point R = R 3 when using R as bifurcation parameter. On the other hand, in Figure 12 (c) the Hopf bifurcation curve is not monotone which means that there may exist two Hopf bifurcation points R = R 2 and R = R 3 . In Figure 12 (c), the Hopf bifurcation curve Γ 1 also intersects with the transcritical bifurcation curve Γ 3 , which means that the Hopf bifurcation points R 2 and R 3 may be larger than R 0 . Figure 12 (d) shows the case of λ ≥ 1. Here the R − µ 1 plane is divided by the Hopf bifurcation curve Γ 1 and the transcritical bifurcation curve Γ 3 . There is no Bogdanov-Takens bifurcation point here as there is no saddle-node bifurcation point, and we suspect all Hopf bifurcations are supercritical so there is no generalized Hopf bifurcation point. Note that Figure 12 does not exhaust all possible situations, but in other cases (as shown in Theorems 2.5 and 2.6), there is no Hopf bifurcation curve then the dynamics is simpler.

5.
Conclusions. We propose an ordinary differential equation model (3) for the interaction between a plant (shrubs or trees) and a resource (water). The biomass death rate is assumed to be decreasing respect to the biomass in order to capture the infiltration feedback [6,7]. We find rich bifurcation structure in (3) including (forward/backward) transcritical bifurcations, saddle-node bifurcations, (subcritical/supercritical) Hopf bifurcations, homoclinic bifurcation, and Bogdanov-Takens bifurcations. The type of bifurcations are listed in Table 1 according to the values of (µ 0 , µ 1 ) in different regions. From analysis and simulations, system (3) admits monostable, bistable or tristable dynamics. There exist bistable regions in which the bare-soil equilibrium coexists with a positive equilibrium, i.e. the extinction of the biomass coexists with the "steady state persistence" of the biomass, or the bare-soil equilibrium coexists with a stable limit cycle, i.e. the extinction of the biomass coexists with the "oscillatory persistence" of the biomass, or the positive equilibrium coexists with a stable limit cycle, and tristable regions in which two equilibria coexist with one limit cycle. In the bistable and tristable parameter regions, the asymptotic biomass is determined by the initial water and biomass density.
A hysteresis loop exists in the model (3) if the rain-fall rate varies. When the rainfall rate is high, the system stays at a high vegetation equilibrium; when the rain-fall rate drops, the vegetation equilibrium decreases and it may switch to a oscillatory state; when the rain-fall rate further decreases, the oscillatory state may shrink back to an equilibrium (Figure 3 (d)-(e)), or may suddenly go to global extinction (Figure 3 (f)-(h)). In the former case, the high vegetation state will also disappear in the saddle-node bifurcation. So either case lead to a catastrophic shift. But in both scenarios, time-periodic pattern occurs, which may be an early indicator of the catastrophic shift. Previous work on catastrophic shifts focus on mostly equilibrium behavior changes in the models (without spatial effect) [14,28,31], and our work indicates that time-periodic pattern is also an important state for such systems. System (3) provides a fundamental theoretical foundation for alternative stable states and oscillations in species-resource systems, and it will also serve as a basic kinetic model for spatiotemporal pattern formation with effect of diffusion, cross-diffusion and nonlocal dispersal.
As mentioned earlier, our model (3) can also be considered as an epidemic model. Indeed one can consider an SIR epidemic model with birth/death and treatment: where α, β, a, b, c > 0. Then methods used in this paper can also be applied to (46). Similar SIR epidemic models with treatment have been considered in [32,34,38,39]. Backward bifurcation has been observed in all these models, and some also studied Hopf bifurcation and existence/nonexistence of limit cycles. But none of these work considered limit cycle bubble, homoclinic bifurcation, and Bogdanov-Takens bifurcation. By using our methods here, one can also obtain a detailed partition of the parameter space and classify the dynamics of (46) into various bifurcation diagrams including "bubble", "heart", "lotus" and "pepper".