Original article Biomath 3 (2014), 1404212, 1–18 B f Volume ░, Number ░, 20░░ BIOMATH ISSN 1314-684X Editor–in–Chief: Roumen Anguelov B f BIOMATH h t t p : / / w w w . b i o m a t h f o r u m . o r g / b i o m a t h / i n d e x . p h p / b i o m a t h / Biomath Forum Mathematical Analysis of a Size Structured Tree-Grass Competition Model for Savanna Ecosystems Valaire Yatat1, Yves Dumont2, Jean Jules Tewa1, Pierre Couteron3 and Samuel Bowong1 1UMMISCO, LIRIMA project team GRIMCAPE, Yaounde, Cameroon Email: yatatvalaire@yahoo.fr; tewajules@gmail.com; sbowong@gmail.com 2CIRAD, Umr AMAP, Montpellier, France 3IRD, Umr AMAP, Montpellier, France Email: yves.dumont@cirad.fr; pierre.couteron@ird.fr Received: 23 October 2013, accepted: 21 April 2014, published: 29 May 2014 Abstract—Several continuous-time tree-grass competition models have been developed to study conditions of long-lasting coexistence of trees and grass in savanna ecosystems according to environmental parameters such as climate or fire regime. In those models, fire intensity is a fixed parameter while the relationship between woody plant size and fire-sensitivity is not systematically considered. In this paper, we propose a mathematical model for the tree-grass interaction that takes into account both fire intensity and size-dependent sensitivity. The fire intensity is modeled by an increasing function of grass biomass and fire return time is a function of climate. We carry out a qualitative analysis that highlights ecological thresholds that summarize the dynamics of the system. Finally, we develop a non-standard numerical scheme and present some simulations to illustrate our analytical results. Keywords-Asymmetric competition, Savanna, fire, continuous-time modelling, qualitative analysis, Non-standard numerical scheme. I. Introduction Savannas are tropical ecosystems characterized by the durable co-occurrence of trees and grasses (Scholes 2003, Sankaran et al. 2005) that have been the focus of researches since many years. Savanna-like vegetations cover extensive areas, especially in Africa and understanding savannas history and dynamics is important both to under- stand the contribution of those areas to biosphere- climate interactions and to sustainably manage the natural resources provided by savanna ecosystems. At biome scale, vegetation cover is known to dis- play complex interactions with climate that often feature delays and feed-backs. For instance any shift from savanna to forest vegetation not only means increase in vegetation biomass and carbon sequestration but also may translate into changes in the regional patterns of rainfall (Scheffer et al. 2003, Bond et al. 2005). In the face of the ongoing global change, it is therefore important to understand how climate along with local factors Citation: Valaire Yatat, Yves Dumont, Jean Jules Tewa, Pierre Couteron, Samuel Bowong, Mathematical Analysis of a Size Structured Tree-Grass Competition Model for Savanna Ecosystems, Biomath 3 (2014), 1404212, http://dx.doi.org/10.11145/j.biomath.2014.04.212 Page 1 of 18 http://www.biomathforum.org/biomath/index.php/biomath http://dx.doi.org/10.11145/j.biomath.2014.04.212 V Yatat et al., Mathematical Analysis of a Size Structured Tree-Grass... drive the dynamics of savannas ecosystems. In many temperate and humid tropical biomes, forest vegetation in known to recover quickly from dis- turbances and woody species are expected to take over herbaceous species. Yet in the dry tropics, it is well-known that grassy and woody species may coexist over decades although their relative proportion may show strong variations (Scholes 2003, Sankaran et al. 2005, 2008). Savanna-like ecosystems are diverse and expla- nations found in the literature about the long- lasting coexistence of woody and grassy vegeta- tion components therefore relate to diverse factors and processes depending on the location and the ecological context. Several studies have pointed towards the role of stable ecological factors in shaping the tree to grass ratio along large-scale gradients of rainfall or soil fertility (Sankaran et al. 2005, 2008). Other studies have rather em- phasized the reaction of vegetation to recurrent disturbances such as herbivory or fire (Langevelde et al. 2003, D’Odorico et al. 2006, Sankaran et al. 2008, Smit et al. 2010, Favier et al. 2012 and references therein). Those two points of view are not mutually-exclusive since both environ- mental control and disturbances may co-occur in a given area, although their relative importance generally varies among ecosystems. Bond et al. (2003) proposed the name of climate-dependent for ecosystems that are highly dependent on cli- matic conditions (rainfall, soil moisture) and fire- dependent or herbivore-dependent for ecosystems which evolution are strongly dependent on fires or herbivores. In a synthesis gathering data from 854 sites across Africa, Sankaran et al. (2005) showed that the maximal observed woody cover appears as water-controlled in arid to semi-arid sites since it directly increase with mean annual precipitation (MAP) while it shows no obvious dependence on rainfall in wetter locations, say above c. 650 mm MAP where it is probably controlled by distur- bance regimes. Above this threshold, fire, grazing and browsing are therefore required to prevent tree canopy closure and allow the coexistence of trees and grasses. Several models using a system of ordinary dif- ferential equations (ODES) have been proposed to depict and understand the dynamics of woody and herbaceous components in savanna-like veg- etation. A first attempt (Walker et al. 1981) was orientated towards semiarid savannas and analyzed the effect of herbivory and drought on the bal- ance between woody and herbaceous biomass. This model refers to ecosystems immune to fire due to insufficient annual rainfall. Indeed, fires in savanna-like ecosystems mostly rely on herba- ceous biomass that has dried up during the dry season. As long as rainfall is sufficient, fire can thus indirectly increase the inhibition of grass on tree establishment in a way far more pervasive than the direct competition between grass tufts and woody seedlings. More recently, several attempts have been made (see Accatino et al. 2010, De Michele et al. 2011 and references therein) to model the dynamics of fire-prone savannas on the basis of the initial framework of Tilman (1994) that used coupled ODES to model the competitive interactions be- tween two kinds of plants. On analogous grounds, Langevelde et al. (2003) have developed a model taking into account fires, browsers, grazers and Walter’s (1971) hypothesis of niche separation by rooting zone depth. Models relying on stochas- tic differential equations have also been used (Baudena et al. 2010). Notably, Accatino et al. (2010) and De Michele et al. (2011) focused on the domain of stability of tree-grass coexistence with respect to influencing ”biophysical” variables (climate, herbivory). However, fire was consid- ered as a forcing factor independent of climate and vegetation, while woody cover was treated as a single variable with no distinction between seedling/saplings which are highly fire sensitive and mature trees which are largely immune to fire damages. The way in which the fundamental, indirect retroaction of grass onto tree dynamics is modeled is therefore to be questioned. In the present paper we therefore a model that differs in this respect. Thus, to take into account the role of fire in Biomath 3 (2014), 1404212, http://dx.doi.org/10.11145/j.biomath.2014.04.212 Page 2 of 18 http://dx.doi.org/10.11145/j.biomath.2014.04.212 V Yatat et al., Mathematical Analysis of a Size Structured Tree-Grass... savanna dynamics, we consider a tree-grass com- partmental model with one compartment for grass and two for trees, namely fire-sensitive individuals (like seedlings, saplings, shrubs) and non-sensitive mature trees. Based on field observations and ex- periments reported by Scholes and Archer (1997) and by Scholes (2003), we develop a system of three coupled non-linear ordinary differential equations (ODES), one equation per vegetation compartment that describes savanna dynamics. In addition, we model fire intensity (i.e. impact on sensitive woody plants) as an increasing function of grass biomass. Compared to existing models, our model aims to properly acknowledge two ma- jor phenomena, namely the fire-mediated negative feedback of grasses onto sensitive trees and the negative feed-back of grown-up, fire insensitive trees on grasses. We therefore explicitly model the asymmetric nature of tree-grass competitive interactions in fire-prone savannas. After some theoretical results of the continuous fire model, though which we highlighted some ecological thresholds that summarize savanna dy- namics and some interesting bistability, we present an appropriate non-standard numerical scheme (see Anguelov et al. 2012, 2013, 2014 and Dumont et al. 2010, 2012) for the model considered and we end with numerical simulations. We show that the fire frequency and the competition parameters are bifurcation parameters which allow the continuous fire model of asymmetric tree-grass competition to converge to different steady states. II. THE CONTINUOUS FIRE MODEL OF ASYMMETRIC TREE-GRASS COMPETITION (COFAC) As we have mentioned before, we consider the class of sensitive tree biomass (TS ), the class of non-sensitive tree biomass (TNS ) and the class of grass biomass (G). We model the fire intensity by an increasing function of grass biomass w(G). To built up our model, we consider the following assumptions. 1) The grass vs. sensitive-tree competition has a negative feedback on sensitive tree dynamics. 2) The grass vs. non sensitive-tree competition has a negative feedback on grass dynamics. 3) After an average time expressed in years, the sensitive tree biomass becomes non sensitive to fire. 4) Fire only impacts grass and sensitive Tree. We also consider the following parameters. • There exists a carrying capacity KT for tree biomass (in tons per hectare, t.ha−1). • There exists a carrying capacity KG for grass biomass (in tons per hectare, t.ha−1). • Sensitive tree biomass is made up from non sensitive tree biomass with the rate γNS (in yr−1) and from existing sensitive tree biomass with the rate γS (in yr−1). • Sensitive tree biomass has a natural death rate µS (in yr−1). • Non sensitive tree biomass has a natural death rate µNS (in yr−1). • f is the fire frequency (in yr−1). • Grass biomass has a natural death rate µG (in yr−1). • 1 ωS is the average time, expressed in year, that a sensitive tree takes to become non sensitive to fire. • 1 ωS + µS is the average time that a tree spends in the sensitive tree class without competition and fires. • σG is the competition rate, for light or/and nutrients, between sensitive tree and grass (in ha.t−1.yr−1). • σNS is the competition rate, for light or/and nutrients, between non sensitive tree and grass (in ha.t−1.yr−1). • ηS is the proportion of sensitive tree biomass that is consumed by fire. • ηG is the proportion of grass biomass that is consumed by fire. Remark 1. Competition parameters σG and σNS are asymmetric, indeed σG inhibits sensitive tree (TS ) growth and there is no reciprocal inhibition; likewise, σNS inhibits grass (G) growth. Based on these ecological premises, and taking Biomath 3 (2014), 1404212, http://dx.doi.org/10.11145/j.biomath.2014.04.212 Page 3 of 18 http://dx.doi.org/10.11145/j.biomath.2014.04.212 V Yatat et al., Mathematical Analysis of a Size Structured Tree-Grass... into account the effect of fire as a forcing con- tinuous in time, which is the classical approach, we propose a model for the savanna vegetation dynamics through a system of three interrelate non-linear equations. The COFAC is given by  dTS dt = (γS TS + γNS TNS ) ( 1 − TS + TNS KT ) −TS (µS + ωS + σGG + f ηS w(G)), dTNS dt = ωS TS −µNS TNS , dG dt = γG ( 1 − G KG ) G − (σNS TNS + f ηG + µG)G, (1) with TS (0) = TS 0 > 0, TNS (0) = TNS 0 ≥ 0 , G(0) = G0 > 0. (2) For this continuous fire model, the fire intensity function w is chosen as a sigmoidal function of grass biomass because we want first to inves- tigate the ecological consequences of the non- linear response of fire intensity to grass biomass, while nearly all published models using differen- tial equations so far assumed a linear response. Non linearity is justified since whenever grass biomass is low fires are virtually absent while fire impact increases rapidly with grass biomass before reaching saturation. Thus, w(G) = G2 G2 + g20 , (3) where G0 = g20 is the value of grass biomass at which fire intensity reaches its half saturation (g0 in tons per hectare, t.ha−1). The feasible region for system (1) is the set Ω defined by Ω = {(TS , TNS , G) ∈ R3+ | 0 ≤ TS + TNS ≤ KT , 0 ≤ G ≤ KG}. III. MATHEMATICAL ANALYSIS A. Existence of equilibria, ecological thresholds and stability analysis We set R01 = γS µNS + γNS ωS µNS (µS + ωS ) and R02 = γG f ηG + µG . 1) Existence of equilibria: Setting the right hand-side of system (1) to zero, straightforward computations lead to the following proposition Proposition 1. System (1) has four kinds of equi- libria • The desert equilibrium point E0 = (0, 0, 0) which always exists. • The forest equilibrium point ET = (T S ; T NS ; 0), with T S = KT µNS ωS + µNS 1 − 1 R01  and T NS = KT ωS ωS + µNS 1 − 1 R01  which is ecologically meaningful whenever R01 > 1. • The point EG = (0, 0, G), with G = KG 1 − 1 R02  , is ecologically meaningful when R02 > 1 The point EG when it exists is the grassland equilibrium. • The savanna equilibrium point ET G = (T∗S , T ∗ NS , G ∗), with T∗S , T ∗ NS and G ∗ given in Appendix A, has an ecological significance whenever R 0 1 > 1, R 0 2 > 1 and 0 < G ∗ < KG 1 − 1 R02  . Remark 2. The number of savanna equilibria depends on the form of the function w. • If w(G) = G, then the COFAC has at most one savanna equilibrium. • If w(G) = G G + G0 (the Holling type II function), then the COFAC has at most two savanna equilibria. Biomath 3 (2014), 1404212, http://dx.doi.org/10.11145/j.biomath.2014.04.212 Page 4 of 18 http://dx.doi.org/10.11145/j.biomath.2014.04.212 V Yatat et al., Mathematical Analysis of a Size Structured Tree-Grass... • If w(G) = G2 G2 + G0 (the Holling type III function), then the COFAC has at most three savanna equilibria. 2) Ecological thresholds interpretation: The qualitative behaviors of the COFAC depend on the following thresholds R01, R 0 2, RG1 = γS µNS + γNS ωS µNS (µS + ωS + σGG + f ηS w(G)) , R T NS 2 = γG f ηG + µG + σNS T NS , where • R01 is the sum of the average amount of biomass produced by a sensitive/young plant, without fires and competition with grass, and the average amount of biomass produced by a mature plant multiplied by the proportion of young plants which reach the mature stage. • RG1 is the sum of the average amount of biomass produced by a sensitive/young plant, in presence of fires and competition with grass, and the average amount of biomass produced by a mature plant multiplied by the proportion of young plants which reach the mature stage. • R02 is the average amount of biomass pro- duced per unit of grass biomass during its whole lifespan in presence of fires and and free from competition with non-sensitive trees. • R T NS 2 is the average biomass produced per unit of grass biomass during its whole lifespan in presence of fires and experiencing competi- tion from non-sensitive trees. Remark 3. The following relations hold R G 1 < R 0 1, R T NS 2 < R 0 2. 3) Stability analysis: Let R = R(G∗) = γG (µNS + ωS ) (γS µNS + γNS ωS ) KG KT µNS ωS σNS (σG + f ηS w′(G?)) . We have the following result: Theorem 1. If R01 < 1 and R 0 2 < 1, then the desert equilibrium E0 is globally asymptotically stable. Proof: See Appendix B. Theorem 2. If R01 > 1, then the forest equilibrium ET exists. • If RT NS2 < 1, then the forest equilibrium ET is locally asymptotically stable. • If R02 < 1, then the forest equilibrium ET is globally asymptotically stable. • If R02 > 1, R T̄NS 2 < 1, R Ḡ 1 > 1 and R < 1, then the forest equilibrium ET is globally asymptotically stable. Proof: See Appendix C. Furthermore, using the same approach as in the proof of Theorem 2, we derive the following results Theorem 3. Suppose R02 > 1 so that the grassland equilibrium EG exists. • If RG1 < 1, then the grassland equilibrium EG is locally asymptotically stable. • If R01 < 1, then the grassland equilibrium EG is globally asymptotically stable. • If R01 > 1, R T̄NS 2 > 1,R Ḡ 1 < 1 and R < 1, then the grassland equilibrium EG is globally asymptotically stable. Theorem 4. Suppose that R01 > 1, R 0 2 > 1 and R > 1. We have the following three cases: • The savanna equilibrium ET G is locally asymptotically stable (LAS) when it is unique. • When there exists two savanna equilibria, one is LAS and the other is unstable. • When there exists three savanna equilibria, two are LAS and one is unstable. Thus System (1) will converges to one of the two stable savanna equilibria depending on initial con- ditions. Proof: See Appendix D. 4) Summary table of the qualitative analysis: The qualitative behavior of system (1) is sum- marized in the following Table in which we present Biomath 3 (2014), 1404212, http://dx.doi.org/10.11145/j.biomath.2014.04.212 Page 5 of 18 http://dx.doi.org/10.11145/j.biomath.2014.04.212 V Yatat et al., Mathematical Analysis of a Size Structured Tree-Grass... only the most realistic, from an ecological point of view, case i.e R01 > 1 and R 0 2 > 1. TABLE I Summary table of the qualitative analysis of system (1) Thresholds E0 ET EG ET G R01 R 0 2 R Ḡ 1 R T̄NS 2 > 1 R > 1 U U U L ? > > > R T̄NS 2 < 1 R > 1 U L U L 1 1 1 R < 1 U G U U RḠ1 R T̄NS 2 > 1 R > 1 U U L L < R < 1 U U G U 1 RT̄NS2 < 1 R > 1 U L L L R < 1 U L L U In Table I, the notations U, L and G stand for unstable, locally asymptotically stable, globally asymptotically stable, respectively, while the notation L? means that we have the global stability if there are no periodic solutions. Remark 4. From an ecological point of view, Lignes 1 to 7 of Table I are interesting because in these cases, one ton of grass biomass will produce during it lifespan at least one ton of grass biomass (R02 > 1) and simultaneously, one ton of tree biomass (sensitive and non sensitive) will produce during it lifespan at least one ton of tree biomass (R01 > 1). Moreover, it is also in these cases that we have the most interesting situations of savanna dynamics, namely bistability cases (Lines 2, 4, 7 in Table 1) and a tristability case (Line 6 in Table 1). IV. NUMERICAL SIMULATIONS Compartmental models are usually solved using standard numerical methods, for example, Euler or Runge Kutta methods included in software package such as Scilab [18] and Matlab [19]. Un- fortunately, these methods can sometimes present spurious behaviors which are not in adequacy with the continuous system properties that they aim to approximate i.e, lead to negative solutions, exhibit numerical instabilities, or even converge to the wrong equilibrium for certain values of the time discretization or the model parameters (see Anguelov et al. 2012, Dumont et al. 2010 for further investigations). For instance, we provided in Appendix E some numerical simulations done with Runge Kutta schemes to illustrate some of its spurious behaviors. In this section, following Anguelov et al. 2012, 2013, 2014 and Dumont et al. 2010, 2012, we perform numerical simulations using an implicit nonstandard algorithm to illus- trate and validate analytical results obtained in the previous sections. A. A nonstandard scheme for the COFAC System (1) is discretized as follows: T k+1NS − T k NS φ(h) = ωS T k+1S −µNS T k+1 NS , Gk+1 − Gk φ(h) = γG ( 1 − G k KG ) Gk+1 −σNS T kNS G k+1 −(µG + f ηG)Gk+1, T k+1S − T k S φ(h) = (γS − (µS + ωS ))T k+1S + γNS T k+1 NS − γS KT T k+1S (T k S + T k NS ) − γNS KT T kNS T k+1 NS − ( γNS KT T kNS + (σGG k + f ηS w(Gk)) ) T k+1S , (4) where the denominator function φ is such that φ(h) = h + O(h2), ∀h > 0. Systems (1) − (2) can be written in the following matrix form: dX dt = A(X)X, X(0) = X0, (5) Biomath 3 (2014), 1404212, http://dx.doi.org/10.11145/j.biomath.2014.04.212 Page 6 of 18 http://dx.doi.org/10.11145/j.biomath.2014.04.212 V Yatat et al., Mathematical Analysis of a Size Structured Tree-Grass... where X = (TNS , G, TS ) ∈ R3+ and A(X) = (Ai j)1≤i, j≤3 with A11 = −µNS , A12 = 0, A13 = ωS , A21 = 0, A22(X) = γG ( 1 − GKG ) − (σNS TNS + f ηG + µG), A23 = 0, A31(X) = γNS ( 1 − TS +TNSKT ) , A32 = 0, A33(X) = γS ( 1 − TS +TNSKT ) − (µS + ωS + σGG + f ηS w(G)). Using (5), the numerical scheme (4) can be rewritten as follows: B(Xk)Xk+1 = Xk, where B(Xk) = (Id3 −φ(h)A(X k)). (6) Thus B(Xk) = 1 + φ(h)µNS 0 −ωS φ(h) 0 1 −φ(h)Ak22 0 −φ(h)Ak31 0 1 −φ(h)A k 33  It suffices now to choose φ(h) such that the matrix B(Xk)) is an M-matrix for all h > 0, which implies that B−1(Xk) is a nonnegative matrix, for all h > 0. In particular, choosing φ such that 1 −φ(h)(γG − (µG + f ηG)) ≥ 0 1 −φ(h)(γS − (µS + ωS )) ≥ 0, (7) lead to positive diagonal terms and nonpositive off diagonal terms. We need to show that B(Xk) is invertible. Obviously 1 − φ(h)Ak22 is a positive eigenvalue. Let us define N, a submatrix of matrix B(Xk), as follows N = ( 1 + φ(h)µNS −ωS φ(h) −φ(h)Ak31 1 −φ(h)A k 33 ) . We already have trace(N) > 0. Then, a direct computation shows that det(N) > 0 if φ(h) is choosen such that 1 −φ(h) ( γS + γNS ωS µNS − (µS + ωS ) ) ≥ 0. Thus we have α(N) > 0, i.e. the eigenvalues have positive real parts, which implies that B(Xk) is invertible. Finally, choosing φ(h) = 1 − e−Qh Q , (8) with Q ≥ max ( γG − (µG + f ηG),γS − (µS + ωS ) + γNS ωS µNS ) , (9) matrix B(Xk) is an M-matrix. Furthermore, assum- ing Xk ≥ 0, we deduce Xk+1 = B−1(Xk)Xk ≥ 0. Lemma 1. Using the expression of φ defined in (8), the numerical scheme (4) is positively stable ( i.e for Xk ≥ 0, we obtain Xk+1 ≥ 0). An equilibrium Xe of the continuous model (1) verifies A(Xe)Xe = 0. Multiplying the above expression by φ(h) and summing with Xe yields (Id3 −φ(h)A(Xe))Xe = Xe, Thus, we deduce that the numerical scheme (4) and the continuous model (1) have the same equi- libria which are (assumed to be) hyperbolic. The dynamics of model (1) can be captured by any number Q satisfying Q ≥ max { |λ|2 2|Re(λ)| } , (10) where λ ∈ sp(J) with Ji j = ∂Ai ∂X j . We also have the following result: Lemma 2. If φ(h) is chosen as in Eqs. (8), (9) and (10), then the numerical scheme (4) is elementary stable ( i.e local stability properties of equilibria are preserved). The proof of Lemma 2 follows the proof of Theorem 2 in Dumont et al., 2010. B. NUMERICAL SIMULATIONS AND BIFUR- CATION PARAMETERS In literature we found the following parameters values Biomath 3 (2014), 1404212, http://dx.doi.org/10.11145/j.biomath.2014.04.212 Page 7 of 18 http://dx.doi.org/10.11145/j.biomath.2014.04.212 V Yatat et al., Mathematical Analysis of a Size Structured Tree-Grass... TABLE II Parameters values found in literature Parameters values References f 0-1 Langevelde et al. 2003 0-2 Accatino et al. 2010 γG 0.4(1) − 4.6(2) (1) Penning de Vries 1982 (2) Menaut et al. 1979 γS + γNS 0.456-7.2 Breman et al. 1995 µS + µNS 0.03-0.3 Accatino et al. 2010 0.4 Langevelde et al. 2003 µG 0.9 Langevelde et al. 2003 ηG 0.1(a)-1(b) (a) Van de Vijver 1999 (b) Accatino et al. 2010 ηS 0.02-0.6 Accatino et al. 2010 ωS 0.05-0.2 Walkeling et al. 2011 We now provide some numerical simulations to illustrate the theoretical results and for discussions. 1) Some monostability and bistability situa- tions: • Monostability. We choose γS γNS γG ηS ηG 1 2 3.1 0.5 0.5 µS µNS µG σNS σG 0.1 0.3 0.3 0.3 0.05 G0 ωS KT KG f 2 0.05 50 12 0.5 yr−1 R01 R 0 2 R G 1 R T NS 2 R 8.8889 5.6364 1.5006 1.2644 2.9424 With the chosen parameters, the savanna equi- librium is stable, i.e. sensitive trees, non- sensitive trees and grasses coexist. Figure 1 presents the 3D plot of the trajectories of system (1). It illustrates that the savanna equilibrium point is stable. Figure 1 also il- lustrates the monostability situation presented in Ligne 1 of Table I. 0 10 20 30 0.511.5 22.53 3.54 0 5 10 T S (t) T NS (t) G (t ) E TG • • • • • • • • Fig. 1. 3D plot of the trajectories of system (1) showing that the savanna equilibrium point ET G point is stable. The red bullets represent different initial conditions. • Bistability – Bistability involving forest and grassland equilibria. The state trajectories of the model will converge to a state depending of initial quantity. We choose γS γNS γG µS µNS µG f 0.4 2 2.1 0.1 0.3 0.3 0.5 yr−1 ηS =0.5, ηG = 0.5, KT =50, KG = 12 R01 R 0 2 R G 1 R T NS 2 R σG σNS 4.8889 3.8182 0.5731 0.9315 0.2958 0.1 0.3 The 3D plot of the trajectories of system (1) is depicted in Figure 2. It clearly ap- pears that the forest and grassland equi- libria are stables. Figure 2 illustrates the bistability situation presented in Ligne 7 of Table I. 0 5 10 15 20 25 30 35 5 10 15 20 0 5 10 T NS (t) T S (t) G (t ) Only grasses (grassland) • • • • • • • • • • • • Only trees (forest) • E T E G • • Fig. 2. 3D plot of the trajectories of system (1) showing that the forest (ET ) and grassland (EG ) equilibria are stable. The green bullets represent different initial conditions. – Bistability involving forest and savanna equilibria. The state trajectories of the model will converge to a state depending on initial quantity. We choose γS γNS γG µS µNS µG f 0.6 2 2.1 0.1 0.3 0.3 0.5 yr−1 ηS =0.5, ηG = 0.5, KT =50, KG = 12 R01 R 0 2 R G 1 R T NS 2 R σG σNS 6.2222 3.8182 1.1156 0.8942 1.2985 0.05 0.3 For these parameters, there exist two savanna equilibria but only one is stable as shown in Figure 3. Figure 3 also illus- trates the bistability situation presented in Ligne 2 of Table I. Biomath 3 (2014), 1404212, http://dx.doi.org/10.11145/j.biomath.2014.04.212 Page 8 of 18 http://dx.doi.org/10.11145/j.biomath.2014.04.212 V Yatat et al., Mathematical Analysis of a Size Structured Tree-Grass... 0 5 10 15 20 25 30 35 40 0 5 10 15 0 2 4 6 8 10 T NS (t) T S (t) G (t ) • • • • • • • • • • • • • Forest Savanna • Fig. 3. 3D plot of the trajectories of system (1) showing that the forest (ET ) and savanna (ET G ) equilibria are stable. The green and red bullets represent different initial conditions. Remark 5. Note also that we didn’t observed periodic behaviors in the previous simulations, considering the set of parameters presented in Table 2, while their existence cannot be completely ruled out by the analytical analysis. C. Some bifurcation parameters In this section we emphasize on some bifurca- tion parameters of system (1) which are such that the COFAC can converge to different steady state depending on the variation of these parameters. • The grass vs. sensitive-tree competition pa- rameter σG is a bifurcation parameter. Figure 4 presents how the system (1) changes from the savanna state to the grassland state as a function of the grass vs. sensitive-tree com- petition parameter σG. We choose γS γNS γG µS µNS µG ηS σNS f 0.4 1 4 0.1 0.3 0.1 0.5 0.3 0.2 yr−1 ηG = 0.5, KT =45, KG = 10 For these parameters values, system (1) un- dergoes a transcritical bifurcation. Indeed, we move from ligne 1 to ligne 5 of Table I. From left to right, (R01 = 3.7778, R 0 2 = 20, R G 1 = 2.2865, RT NS2 = 2.4721, R = 99.3058) → (R01 = 3.7778, R 0 2 = 20, R G 1 = 1.2943, R T NS 2 = 2.4721, R = 5.6803) → (R 0 1 = 3.7778, R02 = 20, R G 1 = 0.9026, R T NS 2 = 2.4721). For the last case, the savanna equilibrium ET G is undefined. 0 50 100 150 200 250 300 0 5 10 15 20 25 Time (year) (a) Sensitive Tree Non Sensitive Tree Grass 0 50 100 150 200 250 300 0 5 10 15 20 25 Time (year) (b) 0 50 100 150 200 250 300 0 5 10 15 20 25 Time (year) (c) Fig. 4. From savanna to grassland as a function of the grass vs. sensitive-tree competition parameter σG . From left to right, the fire period τ = 1f is fixed, while the grass vs. sensitive-tree competition parameter σG increases. In (a) (τ = 5, σG = 0), in (b) (τ = 5, σG = 0.02) and in (c) (τ = 5, σG = 0.04) • The fire period parameter τ = 1f is a bifur- cation parameter. Figure 5 presents a shift of the convergence of system (1) from the forest state to the grassland state as a function of the fire period τ. We choose TABLE III Parameters values for Figures 5 and 6 γS γNS γG µS µNS µG ηS σNS f 0.4 2 2.1 0.1 0.3 0.3 0.5 0.3 yr−1 ηG = 0.5, KT =50, KG = 12 For these parameters values, system (1) un- dergoes a forward bifurcation. Indeed, we move from ligne 3 to ligne 7 of Table I. From left to right, (R01 = 4.8889, R 0 2 = 1.0678, R G 1 = 1.3025, RT NS2 = 0.5720) → (R 0 1 = 4.8889, R02 = 1.3548, R G 1 = 0.5446, R T NS 2 = 0.6453, R = 0.0983) → (R01 = 4.8889, R 0 2 = 1.6154, RG1 = 0.5679, R T NS 2 = 0.6989, R = 0.1201). For the first case, the savanna equilibrium ET G is undefined. 0 50 100 150 200 0 5 10 15 20 25 30 35 Time (year) (a) Sensitive Tree Non Sensitive Tree Grass 0 50 100 150 200 0 0.5 1 1.5 2 2.5 3 3.5 4 Time (year) (b) 0 50 100 150 200 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Time (year) (c) Fig. 5. From forest to grassland as a function of the fire period τ. From left to right, the fire period τ increased, while the sensitive tree-grass competition parameter σG is fixed. In (a) (τ = 0.3, σG = 0.05), in (b) (τ = 0.4, σG = 0.05) and in (c) (τ = 0.5, σG = 0.05) Biomath 3 (2014), 1404212, http://dx.doi.org/10.11145/j.biomath.2014.04.212 Page 9 of 18 http://dx.doi.org/10.11145/j.biomath.2014.04.212 V Yatat et al., Mathematical Analysis of a Size Structured Tree-Grass... Suppose now that fire period is fixed and the grass vs. sensitive-tree competition parameter σG varies. Figure 6 illustrates a shift of the convergence of system (1) from the forest state to the grassland state through a savanna state as a function of the grass vs. sensitive- tree competition parameter σG. For the parameters values in Table III, sys- tem (1) exhibits to bifurcation phenomena: a pitchfork bifurcation and a transcritical bifur- cation. We move from ligne 3 (figure 6 (a), -(b)) to ligne 7 (figure 6 (e), -(f)) through ligne 2 (figure 6 (c), -(d)) of Table I. Indeed, in figure 6 (a), -(b) ET G is undefined, EG is unstable, ET is stable. In in figure 6 (c), -(d), EG remains unstable but we have bistability between ET G and ET : it is a case of pitchfork bifurcation. In figure 6 (d), -(f), ET G becomes unstable and we have a bistability between ET and EG: it is a case of transcritical bifurcation. Values of R01, R 0 2, R G 1 , R T NS 2 and R are given in Appendix F. 0 50 100 150 200 0 10 20 30 40 G (t ), T S (t ), T N S (t ) (a). σ G =0 Sensitive Tree Non Sensitive Tree Grass 0 50 100 150 200 0 10 20 30 40 (b). σG =0.02 0 50 100 150 200 0 5 10 15 (c). σ G =0.04 0 50 100 150 200 250 300 0 2 4 6 8 Time (year) G (t ), T S (t ), T N S (t ) (d). σ G =0.045 Sensitive Tree Non Sensitive Tree Grass 0 100 200 300 0 5 10 Time (year) (e). σ G =0.05 0 100 200 300 0 5 10 Time (year) (f). σ G =0.06 Fig. 6. From forest to grassland, with a transition through a savanna state, as a function of the sensitive tree-grass competition parameter. From left to right, the fire period τ is fixed at 4, while the grass vs. sensitive-tree competition parameter σG increased. • The grass vs. non sensitive-tree competition parameter σNS is a bifurcation parameter. A shift of the convergence of system (1) from the grassland state to the forest state as a function of the grass vs. non sensitive-tree competition parameter σNS is depicted in Fig. 7. We choose γS γNS γG µS µNS µG σG ηS f 0.4 1 4 0.1 0.3 0.1 0.05 0.5 yr−1 ηG = 0.5, KT =45, KG = 10 For these parameters values, system (1) ex- hibits a pitchfork bifurcation. We move from ligne 5 (figure 7 (a), -(b)) to ligne 7 (figure 7 (c), -(d)) of Table I. Indeed, in figure 7 (a), -(b) ET G is undefined, ET is unstable, EG is stable. In figure 7 (c), -(d), ET G exits but it is unstable and we have bistability between EG and ET : it is a case of pitchfork bifurcation. Values of R01, R 0 2, R G 1 , R T NS 2 and R are given in Appendix G. 0 100 200 300 0 5 10 15 20 25 Time (year) G (t ), T S (t ), T N S (t ) (a). τ =0.5, σ NS =0.5 0 100 200 300 0 5 10 15 20 25 30 Time (year) G (t ), T S (t ), T N S (t ) (b). τ =0.5, σ NS =0.6 0 50 100 150 0 5 10 15 20 25 30 Time (year) G (t ), T S (t ), T N S (t ) (c). τ =0.5, σ NS =0.65 0 50 100 150 0 5 10 15 20 25 30 Time (year) G (t ), T S (t ), T N S (t ) (d). τ =0.5, σ NS =0.7 Sensitive Tree Non Sensitive Tree Grass Sensitive Tree Non Sensitive Tree Grass Sensitive Tree Non Sensitive Tree Grass Sensitive Tree Non Sensitive Tree Grass Fig. 7. From grassland to forest as a function of the grass vs. non sensitive-tree competition parameter σNS . From left to right, the fire period τ is fixed, while the grass vs. non sensitive-tree competition parameter σNS increases. V. Conclusion and discussion In this work, we present and analyze a new mathematical model to study the interaction of tree and grass that explicitly makes fire intensity dependent on the grass biomass and distinguishes two levels of fire sensitivity within the woody biomass (implicitly relating to plant size and bark thickness). Fire was considered as a time- continuous forcing as in several existing models (Langevelde et al. 2003, Accatino et al. 2010, De Michele et al. 2011 and reference therein) with a constant frequency of fire return that can be interpreted as mainly expressing an external forcing to the tree-grass system from climate and human practices. What is novel in our model is that fire impact on tree biomass is modeled as a non-linear function w of the grass biomass. Using a non-linear function is to our knowledge only found in Staver et al. 2011. But this latter model made peculiar assumptions and does not predict grassland and forest as possible equilibria (only desert and savanna). The advantage of a non-linear function is that it can account for the absence of fire at low biomass. As a consequence and Biomath 3 (2014), 1404212, http://dx.doi.org/10.11145/j.biomath.2014.04.212 Page 10 of 18 http://dx.doi.org/10.11145/j.biomath.2014.04.212 V Yatat et al., Mathematical Analysis of a Size Structured Tree-Grass... although keeping the same modeling paradigm as in Langevelde et al. 2003, Accatino et al. 2010, De Michele et al. 2011, we reached different results and predictions. Distinguishing fire sensitive vs. fire insensi- tive woody biomass lead to three variables ex- pressing fractions of the above ground phy- tomass, namely grass and both fire-sensitive and - insensitive woody vegetation. It featured three cou- pled, non-linear ordinary differential equations.As several existing models (Baudena et al. 2010, Staver et al. 2011), our model acknowledges two major phenomena that regulate savanna dynamics, namely the fire-mediated negative feedback of grasses onto sensitive trees and the negative feed- back of grown-up, fire insensitive trees on grasses. We therefore explicitly model the asymmetric na- ture of tree-grass competitive interactions in fire- prone savannas. The analytical study of the model reveals three possible equilibria excluding tree-grass coexis- tence (desert, grassland, forest) along with equilib- ria for which woody and grassy components show durable coexistence (i.e. savanna vegetation). The number of such equilibrium points depends on the function used to model the increase of fire intensity with grass biomass(see Remark 2); for our model, we can have at most three savanna equilibria. We identified four ecologically meaningful thresholds that defined in parameter space regions of monos- tability, bistability as in Accatino et al. 2010, De Michele et al. 2011 and tristability with respect to the equilibria. Tristability of equilibria may mean that shifts from one stable state to another may often be less spectacular that hypothesized from previous models and that scenarios of vegetation changes may be more complex. The model features some parameters that have been analytically identified as liable to trigger bifurcations (i.e., the state variables of the model converges to different steady states), notably pa- rameters σNS and σG of asymmetric competi- tion that embody the depressing influence of fire insensitive trees on grasses and of grasses on sensitive woody biomass respectively. Since tree- grass asymmetric competition is largely mediated by fire, this finding of the role of those two parameters is not intuitive and is the result of the modeling effort and of the analytical anal- ysis. Since such parameters that quantify direct interactions between woody and grassy compo- nents appear crucial to understand the tree-grass dynamics in savanna ecosystems and for enhanced parameter assessment, they could be the focus of straightforward field experiments that would not request manipulating fire regime. Another bifur- cation parameter is the fire frequency, f , (or fire period parameter τ = 1f ) which has been assumed to be an external forcing parameter that integrates both climatic and human influences. Frequent fires preclude tree-grass coexistence and turn savannas into grasslands. In the wettest situations, or under subequatorial climates, very high fire frequencies (above one fire per year) seem to be needed to prevent the progression of forests over savannas (unpublished data of experiments carry out at La Lopé National Park in Gabon). However, it is questionable to model fire as a continuous forcing that regularly removes frac- tions of fire sensitive biomass. Indeed, several months can past between two successive fires, such that fire may be considered as an instantaneous perturbation of the savanna ecosystem. Several recent papers have proposed to model fires as stochastic events while keeping the continuous- time differential equation framework (Beckage et al. 2011) or using time discrete matrix models (Accatino & De Michele 2013). But in all those examples, fire characteristics remain mainly a lin- ear function of grass biomass. Another framework that we will explore in a forthcoming work in order to acknowledge the discrete nature of fire events is based on system of impulsive differential equations (Lakshmikantham et al. 1989, Bainov and Simeonov 1993). References [1] F. Accatino, C. De Michele, Humid savanna-forest dy- namics: a matrix model with vegetation-fire interactions and seasonality. Eco. Mod. 265, pp. 170-179, 2013. Biomath 3 (2014), 1404212, http://dx.doi.org/10.11145/j.biomath.2014.04.212 Page 11 of 18 http://dx.doi.org/10.11145/j.biomath.2014.04.212 V Yatat et al., Mathematical Analysis of a Size Structured Tree-Grass... [2] F. Accatino, C. De Michele, R. Vezzoli, D. Donzelli, R. Scholes, ”Tree-grass co-existence in savanna: interac- tions of rain and fire.” J. Theor. Biol. 267, pp. 235-242, 2010. http://dx.doi.org/10.1016/j.jtbi.2010.08.012 [3] R. Anguelov, Y. Dumont, J.M-S. Lubuma and M. Shillor, Comparison of some standard and nonstan- dard numerical methods for the MSEIR epidemiological model, In: T. Simos, G. Psihoyios, Ch Tsitouras (eds), Proceedings of the International Conference of Numer- ical Analysis and Applied Mathematics, Crete, Greece, 18-22 September 2009, American Institute of Physics Conference Proceedings-AIP 1168, Volume 2, 2009, pp. 1209-1212. http://dx.doi.org/10.1063/1.3241285 [4] R. Anguelov, Y. Dumont, and J.M.-S. Lubuma, On nonstandard finite difference schemes in biosciences. AIP Conf. Proc. 1487, pp. 212-223, 2012. http://dx.doi.org/10.1063/1.4758961 [5] R. Anguelov, Y. Dumont, J.M.-S. Lubuma, and E. Mure- ithi, ”Stability Analysis and Dynamics Preserving Non- Standar Finite Difference Schemes for Malaria Model”, Mathematical Population Studies, 20 (2), pp. 101-122, 2013. http://dx.doi.org/10.1080/08898480.2013.777240 [6] R. Anguelov, Y. Dumont, J.M.-S. Lubuma, and M. Shillor, Dynamically consistent nonstandard finite dif- ference schemes for epidemiological Models. Journal of Computational and Applied Mathematics, 255, pp. 161-182, 2014. http://dx.doi.org/10.1016/j.cam.2013.04.042 [7] D.D. Bainov, and P.S. Simeonov, Impulsive Differential Equations: Periodic Solutions and Applications. long- man, England, 1993. [8] M. Baudena, F. D’Andrea, A. Provenzale, An idealized model for tree-grass coexistence in savannas: the role of life stage structure and fire disturbances. J. Ecol. 98, pp. 74-80, 2010. http://dx.doi.org/10.1111/j.1365-2745.2009.01588.x [9] B. Beckage, L.J. Gross, W.J. Platt, Grass feedbacks on fire stabilize savannas. Eco. Mod. 222, pp. 2227-2233, 2011. http://dx.doi.org/10.1016/j.ecolmodel.2011.01.015 [10] W.J. Bond, G.F. Midgley, F.I. Woodward, What controls SouthAfrican vegetation-climate or fire? S. Afr. J. Bot. 69, pp. 79-91, 2003. [11] W.J. Bond, F.I. Woodward, G.F. Midgley, The global distribution of ecosystems in a world without fire. New Phytologist, 165, pp. 525-538, 2005. http://dx.doi.org/10.1111/j.1469-8137.2004.01252.x [12] H. Breman and J.J. Kessler. Woody plants in agroe- cosystems of semi-arid regions. With an emphasis on the Sahelian countries. Advanced series in Agricultural 23, Springer-Verlag, Berlin. 1995. [13] C. De Michele, F. Accatino, R. Vezzoli, R.J. Scholes, Savanna domain in the herbivores-fire parameter space exploiting a tree-grass-soil water dynamic model. J. Theor. Biol. 289, pp. 74-82, 2011. http://dx.doi.org/10.1016/j.jtbi.2011.08.014 [14] Y. Dumont, J.C. Russell, V. Lecomte and M. Le Corre, Conservation of endangered endemic seabirds within a multi-predatot context:The Barau’s petrel in Reunion island. Natural Ressource Modelling, 23, pp. 381-436, 2010. http://dx.doi.org/10.1111/j.1939-7445.2010.00068.x [15] Y. Dumont, J.M. Tchuenche, Mathematical Studies on the Sterile Insect Technique for the Chikungunya Dis- ease and Aedes albopictus. Journal of mathematical Biology, 65 (5), pp. 809-854, 2012. http://dx.doi.org/10.1007/s00285-011-0477-6 [16] P. D’Odorico, F. Laio, and L.A. Ridolfi, probabilistic analysis of fire-induced tree-grass coexistence in savan- nas. The American Naturalist, 167, pp. E79-E87, 2006. [17] C. Favier, J. Aleman, L. Bremond, M.A. Dubois, V. Freycon, and J.M. Yangakola, Abrupt shifts in African savanna tree cover along a climatic gradient. Global Ecology and Biogeography, 21, pp. 787-797, 2012. http://dx.doi.org/10.1111/j.1466-8238.2011.00725.x [18] http://www.scilab.org [19] http://www.mathworks.com [20] V. Lakshmikantham, D.D. Bainov and P.S. Simeonov, Theory of Impulsive Differential Equations, World Sci- entific, Singapore, 1989. [21] V.F. Langevelde, C. van de Vijver, L. Kumar, J. van de Koppel, N. de Ridder, J. van Andel, et al., Effects of fire and herbivory on the stability of savanna ecosystems. Ecology, 84 (2), pp. 337-350, 2003. [22] M.Y. Li and L. Wang. A criterion for stability of Matrices. J. Math. Ana. App. 225, pp. 249-264, 1998. [23] J.C. Menaut and J. Csar. Structure and primary produc- tivity of Lamto savannas, Ivory Coast. Ecology. 60, pp. 1197-1210, 1979. [24] F.W.T. Penning de Vries and M.A. Djiteye. La pro- ductivité des paturages sahéliens. Une étude des sols, des végétations et de l’exploitation de cette ressource naturelle. PUDOC, Wageningen, 1982. [25] M. Sankaran, J. Ratnam, and N. Hanan. Woody cover in African savannas: the role of resources, fire and herbivory. Global Ecology and Biogeography, 17, pp. 236-245, 2008. http://dx.doi.org/10.1111/j.1466-8238.2007.00360.x [26] M. Sankaran, N.P. Hanan, R.J. Scholes, J. Ratnam, D.J. Augustine, B.S. Cade, J. Gignoux, S.I. Higgins, X. LeRoux, F. Ludwig, J. Ardo, F. Banyikwa, A. Bronn, G. Bucini, K.K. Caylor, M.B. Coughenour, A. Diouf, W. Ekaya, C.J. Feral, E.C. February, P.G.H. Frost, P. Hiernaux, H. Hrabar, K.L. Metzger, H.H.T. Prins, S. Ringrose, W. Sea, J. Tews, J. Worden, N. Zambatis. De- terminants of woody coverin African savannas. Nature, 438, pp. 846-849, 2005. http://dx.doi.org/10.1038/nature04070 [27] M. Scheffer, and Stephen R. Carpenter. Catastrophic regime shifts in ecosystems: linking theory to observa- tion. TRENDS in Ecology and Evolution. Vol.18 No.12, pp. 648-656, December 2003. Biomath 3 (2014), 1404212, http://dx.doi.org/10.11145/j.biomath.2014.04.212 Page 12 of 18 http://dx.doi.org/10.1016/j.jtbi.2010.08.012 http://dx.doi.org/10.1063/1.3241285 http://dx.doi.org/10.1063/1.4758961 http://dx.doi.org/10.1080/08898480.2013.777240 http://dx.doi.org/10.1016/j.cam.2013.04.042 http://dx.doi.org/10.1111/j.1365-2745.2009.01588.x http://dx.doi.org/10.1111/j.1469-8137.2004.01252.x http://dx.doi.org/10.1016/j.jtbi.2011.08.014 http://dx.doi.org/10.1111/j.1939-7445.2010.00068.x http://dx.doi.org/10.1007/s00285-011-0477-6 http://dx.doi.org/10.1111/j.1466-8238.2011.00725.x http://www.scilab.org http://www.mathworks.com http://dx.doi.org/10.1111/j.1466-8238.2007.00360.x http://dx.doi.org/10.1038/nature04070 http://dx.doi.org/10.11145/j.biomath.2014.04.212 V Yatat et al., Mathematical Analysis of a Size Structured Tree-Grass... http://dx.doi.org/10.1016/j.tree.2003.09.002 [28] R.J. Scholes, Convex relationships in ecosystems con- taining mixtures of trees and grass. Environnemental and Ressource Economics, 26, pp. 559-574, 2003. [29] R.J. Scholes, S.R. Archer, Tree-Grass interactions in savannas. Annu. Rev. Ecol. Syst. 28, pp. 517-544, 1997. [30] I.P.J. Smit, G. Asner, N. Govender, T. Kennedy- Bowdoin, D. KNAPP, and J. Jacobson, Effects of fire on woody vegetation structure in African savanna. Eco- logical Applications, 20 (7), pp. 1865-1875, 2010. [31] A.C. Staver, S. Archibald, S. Levin, Tree cover in sub- Saharan Africa: Rainfall and fire constrain forest and savanna as alternative stable states. Ecology. 92(5), pp. 1063-1072, 2011. [32] D. Tilman, Competition and biodiversity in spatially structured habitats. Ecology, 75, pp. 2-16, 1994. [33] C.A. Van de Vijver, Foley and H. Olff, Changes in the woody component of an East African savanna during 25 years. Journal of Tropical Ecology 15, pp. 545-564, 1999. [34] J.L. Walkeling, A.C. Staver and W.J. Bond, Simply the best: the transition of savanna saplings to trees. Oikos, 120, pp. 1448-1451, 2011. [35] B. Walker, D. Ludwig, C.S. Holling, and R.M. Peter- man, Stability of semi-arid savanna grazing systems. Journal of Ecology, 69, pp. 473-498, 1981. [36] H. Walter, Ecology of tropical and subtropical vegeta- tion. Oliver and Boyd, Edinburgh, UK., 1971. Appendix A: Expressions of T∗S , T ∗ NS and G ∗ After straightforward but long computation, we show that T∗NS = γG σNS 1 − 1 R02 − G∗ KG  , T∗S = µNS ωS T∗NS , 1 KT ( 1 + ωS µNS ) T∗S = 1 − 1 R01 − µNS (σGG∗ + f ηS w(G∗)) γS µNS + γNS ωS , where G∗ is solution of w(G) = AG + B = F(G), (11) with A = 1 f ηS ( γG(ωS + µNS )(γS µNS + γNS ωS ) KG KT σNS µNS ωS −σG ) , B = (γS µNS + γNS ωS ) [ 1 − 1 R01 − γG (ωS +µNS ) KT σNS ωS ( 1 − 1 R02 )] f ηS µNS . We summarize the problem of existence of solutions of equation (11) in the following Table TABLE IV Existence of solutions of equation (11) A B Number of solutions > 0 > 0 0 or 2 solutions < 0 1 or 3 solutions < 0 > 0 1 solution < 0 0 solution Note that solutions G∗ of (11) that give rise to savanna equilibria must satisfy 0 < G∗ < KG 1 − 1 R02  . Appendix B: Proof of Theorem 1 let R00 = γS µS + ωS + µNS . In a matrical writing, System (1) reads as dX dt = A(X)X < Amax(X)X, (12) with X = (TS , TNS , G) ∈ R3+, A(X) = (Ai j)1≤i, j≤3 with A11 = γS ( 1 − TS +TNSKT ) − (µS + ωS + σGG + f ηS w(G)), A12 = γNS ( 1 − TS +TNSKT ) , A13 = 0, A21 = ωS , A22 = −µNS , A23 = 0, A31 = 0, A32 = 0, A33 = γG ( 1 − GKG ) −(σNS TNS + f ηG +µG). and Amax(X) =  γS −µS −ωS γNS 0 ωS −µNS 0 0 0 γG 1 − 1 R02   =( A B C D ) , with A = ( γS −µS −ωS γNS ωS −µNS ) , B =( 0 0 ) , C = ( 0 0 ) , and D = γG 1 − 1 R02  . Matrix Amax(X) is a Metzler matrix ( i.e all its off- diagonal terms are nonnegative) and α(Amax(X)) ≤ 0 if α(A) ≤ 0 and α(D) ≤ 0 where α denotes the stability modulus. Moreover, for matrix D, α(D) ≤ 0 if R 0 2 < 1. (13) Biomath 3 (2014), 1404212, http://dx.doi.org/10.11145/j.biomath.2014.04.212 Page 13 of 18 http://dx.doi.org/10.1016/j.tree.2003.09.002 http://dx.doi.org/10.11145/j.biomath.2014.04.212 V Yatat et al., Mathematical Analysis of a Size Structured Tree-Grass... For matrix A, α(A) ≤ 0 if trace(A) < 0 and det(A) > 0. trace(A) = γS −µS −ωS −µNS = γS 1 − 1 R00  . (14) det(A) = µNS (µS + ωS ) − (µNS γS + ωS γNS ) = µNS (µS + ωS )(1 −R01). (15) Furthermore, R00 < R 0 1 (16) Thus, from relations (13), (14), (15) and (16) we deduce that the desert equilibrium (0; 0; 0) is globally asymptotically stable whenever R01 < 1 and R02 < 1. Appendix C: Proof of Theorem 2 If R01 > 1, then the forest equilibrium ET exists. • Using the Jacobian matrix of system (1) at ET , one can prove that ET is locally asymp- totically stable if RT NS2 < 1. • The solution G of system (1) verify dG dt ≤ (γG − ( f ηG + µG))G, ≤ γG 1 − 1 R02  . (17) So, if R02 < 1, then lim t→+∞ G(t) = 0. (18) Moreover, the solutions TS and TNS of system (1) admit as a limit system, the system: dTS dt = (γS TS + γNS TNS ) ( 1 − TS + TNS KT ) −TS (µS + ωS ) = F1(TS , TNS ), dTNS dt = ωS TS −µNS TNS = F2(TS ; TNS ). (19) Now, let h(TS , TNS ) = T−1S . Then, one has ∂F1 h ∂TS + ∂F2 h ∂TNS = −γNS TNS T 2S ( 1 − TS +TNSKT ) − 1 KT TS (γS TS + γNS TNS ) −µNS T−1S . Furthermore, we have ∂F1h ∂TS + ∂F2h ∂TNS < 0 in Ω◦2 where Ω2 = {(TS , TNS ) ∈ R2+ | 0 ≤ TS + TNS ≤ KT}, and by the Bendixson-Dulac theorem, we deduce that system (19) don’t admits a periodic solution in Ω2. Moreover, the equilibrium (T̄S , T̄NS ) exists if R01 > 1 and using the Jacobian matrix of system (19), we deduce that (T̄S , T̄NS ) is locally asymptotically stable and then, glob- ally asymptotically stable since there is no periodic solution. Thus, if R02 < 1, then one has lim t→+∞ (TS , TNS , G)(t) = ET . • Suppose that R01 = γS µNS + γNS ωS µNS (µS + ωS ) > 1 and R02 = γG f ηG + µG > 1, then equilibria (T̄, T̄NS , 0), (0, 0, Ḡ) and (T ?S , T ? NS , G ?) are defined. The Jacobian matrix of system (1) at an arbitrarily equilibrium point is J =  J11 J12 J13 J21 J22 0 0 J32 J33  , where J11 = γS ( 1 − X+YKT ) − 1 KT (γS X + γNS Y ) −µS −ωS −σGZ − f ηS w(Z), J12 = γNS ( 1 − X+YKT ) − 1 KT (γS X + γNS Y ), J13 = −σG X − X f ηS w′(Z), J21 = ωS , J22 = −µNS , J32 = −σNS Z, J33 = γG − 2 γG KG Z −σNS Y − f ηG −µG. The second additive compound matrix of J is J[2] =  J11 + J22 0 −J13 J32 J11 + J33 J12 0 J21 J22 + J33  . (20) From the Jacobian matrix of system (1), the equilibria (0, 0, Ḡ) and (T ?S , T ? NS , G ?) are unstable if RḠ2 > 1 and R < 1. In the sequel, we suppose that RT̄NS2 < 1 to process with the discussion. Biomath 3 (2014), 1404212, http://dx.doi.org/10.11145/j.biomath.2014.04.212 Page 14 of 18 http://dx.doi.org/10.11145/j.biomath.2014.04.212 V Yatat et al., Mathematical Analysis of a Size Structured Tree-Grass... The second additive compound matrix (20) at the equilibrium (T̄S , T̄NS , 0) is J[2](T̄S , T̄NS , 0) = J11 + J22 0 σT̄S 0 J11 + J33 J12 0 J21 J22 + J33  (T̄S ,T̄NS ,0) . Let B = ( J11 + J33 J12 J21 J22 + J33 ) . Then, a simple calculation gives (J11 + J22)(T̄S ,T̄NS ,0) = γS ( 1 − T̄S + T̄NS KT ) − 1 KT (γS T̄S + γNS T̄NS ) −µS −ωS −µNS . Using the relations −µS −ωS = − ( γS + γNS ωS µNS ) ( 1 − T̄S + T̄NS KT ) and ( 1 − T̄S + T̄NS KT ) > 0, we have (J11 + J22)(T̄S ,T̄NS ,0) = − 1 KT T̄S ( γS + γNS ωS µNS ) −γNS ωS µNS ( 1 − T̄S + T̄NS KT ) −µNS < 0. Since J11 + J22 < 0 and R T̄NS 2 < 1, one has tr(B) = (J11 + J22 + 2J33)(T̄S ;T̄NS ;0), = J11 + J22 + 2γG 1 − 1 R T̄NS 2  < 0. Also, if RT̄NS2 < 1, one has J11 J33 = ( − γG KT T̄S ( γS + γNS ωS µNS ) −γGγNS ωS µNS ( 1 − T̄S +T̄NSKT )) × ( 1 − 1 R T̄NS 2 ) > 0, and J33(J22 + J33) = γG 1 − 1 R T̄NS 2  −µNS + γG 1 − 1 R T̄NS 2   > 0. With this in mind, we have det(B) = (J11 + J33)(J22 + J33) − J12 J21, = J11 J22 − J21 J12 + J11 J33 +J33(J22 + J33), = µNS KT T̄S ( γS + γNS ωS µNS ) + ωS KT (γS T̄S + γNS T̄NS ) +J11 J33 + J33(J22 + J33) > 0. Thus, if RT̄NS2 < 1, one has (J11 + J22)(T̄S ,T̄NS ,0) < 0, tr(B) < 0 and det(B) > 0. This implies that s(J[2](T̄S , T̄NS , 0)) < 0 where s denotes the stability modulus. Following Theorem 3.3 in Li and Wang 1998, we can deduce that there is no hopf bifurcation points for J(T̄S , T̄NS , 0). Since R T̄NS 2 < 1, the equilibrium point (T̄S , T̄NS , 0) is locally asymptotically stable and one can conclude that this equilibrium point is globally asymptotically stable if RT̄NS2 < 1,R Ḡ 2 > 1 and R < 1. This completes the proof. Appendix D: Proof of Theorem 4 Suppose that the savanna equilibrium ET G ex- ists. The Jacobian matrix of system (1) at ET G is J =  J11 J12 J13 J21 J22 0 0 J32 J33  , where J11 = γS ( 1 − TS +TNSKT ) − 1 KT (γS TS + γNS TNS ) −µS −ωS −σGG − f ηS w(G), J12 = γNS ( 1 − TS +TNSKT ) − 1 KT (γS TS + γNS TNS ), J13 = −σGTS − TS f ηS w′(G), J21 = ωS , J22 = −µNS , J32 = −σNS G, J33 = − γG KG G. Let A1 = −J11 J22 J33, A2 = J21 J12 J33, A3 = −J21 J32 J13, C1 = −J11 − J22 − J33, C2 = A1 + A2 + A3, C3 = J11 J33 + J11 J22 − J21 J12 + J22 J33 − C2 C1 . Biomath 3 (2014), 1404212, http://dx.doi.org/10.11145/j.biomath.2014.04.212 Page 15 of 18 http://dx.doi.org/10.11145/j.biomath.2014.04.212 V Yatat et al., Mathematical Analysis of a Size Structured Tree-Grass... Note that, by the Routh-Hurwitz theorem, the savanna equilibrium ET G is locally asymptotically stable if C1 > 0, C2 > 0 and C3 > 0. Moreover, components of the savanna equilibrium ET G satisfy TNS = ωS µNS TS , −µS −ωS −σGG − f ηS w(G) = − ( γS + γNS ωS µNS ) ×( 1 − TS +TNSKT ) , thus, C1 = −J11 − J22 − J33, = 1KT (γS TS + γNS TNS ) + γNS ωS µNS ( 1 − TS +TNSKT ) +µNS + γG KG G, > 0. C3 = J11 J33 + J11 J22 − J21 J12 + J22 J33 − C2 C1 , = 1KT (µNS + ωS )(γS TS + γNS TNS ) + γGG KG ( 1 KT (γS TS + γNS TNS ) +γNS ωS µNS ( 1 − TS +TNSKT ) + µNS ) − C2 C1 , = 1KT (µNS + ωS )(γS TS + γNS TNS ) − γG G KG KT (µNS +ωS )(γS TS +γNS TNS ) C1 + ωS σNS G(σG TS f ηS w′(G)TS ) C1 + γGG KG ( 1 KT (γS TS + γNS TNS ) +γNS ωS µNS ( 1 − TS +TNSKT ) + µNS ) , = 1KT (µNS + ωS )(γS TS + γNS TNS ) ( 1 − γG G KG C1 ) + ωS σNS G(σG TS f ηS w′(G)TS ) C1 + γGG KG ( 1 KT (γS TS + γNS TNS ) +γNS ωS µNS ( 1 − TS +TNSKT ) + µNS ) , = 1KT C1 (µNS + ωS )(γS TS + γNS TNS )×( 1 KT (γS TS + γNS TNS ) +γNS ωS µNS ( 1 − TS +TNSKT ) + µNS ) + ωS σNS G(σG TS f ηS w′(G)TS ) C1 + γGG KG ( 1 KT (γS TS +γNS TNS ) + γNS ωS µNS ( 1 − TS +TNSKT ) + µNS ) , C3 > 0. C2 = A1 + A2 + A3, = γGG KG KT (µNS + ωS )(γS TS + γNS TNS ) −ωS σNS G(σGTS + f ηS w′(G)), = ωS σNS GTS ( γG KG KT ωS σNS (µNS + ωS )×( γS + γNS ωS µNS ) −σG − f ηS w′(G) ) , = ωS σNS GTS (σG + f ηS w′(G))(R− 1). Thus, C2 > 0 if and only if R > 1. Finally, we deduce that the savanna equilibrium ET G, when it is unique, is locally asymptotically stable if R = R(G) > 1. The first part of Theorem 4 holds. One should note that C2 > 0 means that the slope of w (the sigmoidal function) is less than the slope of F where F is given by (11). Furthermore, by using relation (11) we deduce part 2 and part 3 of Theorem 4 graphically as follow Fig. 8. There exist two savanna equilibria but one is stable and the other is unstable. Fig. 9. There exist three savanna equilibria two are stable and one is unstable. Thus system (1) will converge to one of the two stable equilibria depending on initial conditions. Biomath 3 (2014), 1404212, http://dx.doi.org/10.11145/j.biomath.2014.04.212 Page 16 of 18 http://dx.doi.org/10.11145/j.biomath.2014.04.212 V Yatat et al., Mathematical Analysis of a Size Structured Tree-Grass... 0 5 10 15 20 25 30 35 40 45 50 −5 0 5 10 15 20 Time (year) G (t ), T S (t ), T N S (t ) Runge Kutta order (2,3): ode23 Sensitive Tree Non Sensitive Tree Grass 0 5 10 15 20 25 30 35 40 45 50 0 2 4 6 8 10 12 14 16 18 20 Time (year) G (t ), T S (t ), T N S (t ) Nonstandard Numerical Scheme Sensitive Tree Non Sensitive Tree Grass Fig. 10. For these figure, µG = 0.4. The ODE’s routine which is a standard numerical algorithms shows a spurious negative solutions. Appendix E: Spurious behaviors of Runge Kutta methods to approximate solutions of system (1) For the following figures we choose γS γNS γG ηS ηG 0.1 2 0.6 0.5 0.5 µS µNS µG σNS σG 0.3 0.3 0.02 0.05 G0 ωS KT KG f 2 0.05 50 12 0.5 Other examples of spurious solutions given by standard methods are also given in (Anguelov et al. 2009). Appendix F: Values ofR01, R 0 2, R G 1 , R T NS 2 andR in figure 6 • Figure 6 (a): R01 R 0 2 R G 1 R T NS 2 4.8889 4.9412 2.6928 0.9863 • Figure 6 (b): 0 50 100 150 −5 0 5 10 15 20 Time (year) G (t ), T S (t ), T N S (t ) Runge Kutta order (4,5): ode45 Sensitive Tree Non Sensitive Tree Grass 0 50 100 150 0 2 4 6 8 10 12 14 16 18 20 Time (year) G (t ), T S (t ), T N S (t ) Nonstandard Numerical Scheme Sensitive Tree Non Sensitive Tree Grass Fig. 11. For these figure, µG = 0.5. The ODE’s routine which is a standard numerical algorithms shows again a spurious negative solutions. R01 R 0 2 R G 1 R T NS 2 4.8889 4.9412 1.5813 0.9861 • Figure 6 (c): R01 R 0 2 R G 1 R T NS 2 R 4.8889 4.9412 1.1193 0.9861 1.3984 • Figure 6 (d): R01 R 0 2 R G 1 R T NS 2 R 4.8889 4.9412 1.0431 0.9861 1.2980 • Figure 6 (e): R01 R 0 2 R G 1 R T NS 2 R 4.8889 4.9412 0.9766 0.9861 0.5936 • Figure 6 (f): R01 R 0 2 R G 1 R T NS 2 R 4.8889 4.9412 0.8662 0.9861 0.5746 Biomath 3 (2014), 1404212, http://dx.doi.org/10.11145/j.biomath.2014.04.212 Page 17 of 18 http://dx.doi.org/10.11145/j.biomath.2014.04.212 V Yatat et al., Mathematical Analysis of a Size Structured Tree-Grass... Appendix G: Values ofR01, R 0 2, R G 1 , R T NS 2 andR in figure 7 • Figure 7 (a): R01 R 0 2 R G 1 R T NS 2 3.7778 3.6364 0.3840 1.1549 • Figure 7 (b): R01 R 0 2 R G 1 R T NS 2 3.7778 3.6364 0.3840 1.0162 • Figure 7 (c): R01 R 0 2 R G 1 R T NS 2 R 3.7778 3.6364 0.3840 0.9587 0.2066 • Figure 7 (d): R01 R 0 2 R G 1 R T NS 2 R 3.7778 3.6364 0.3840 0.9073 0.1453 Biomath 3 (2014), 1404212, http://dx.doi.org/10.11145/j.biomath.2014.04.212 Page 18 of 18 http://dx.doi.org/10.11145/j.biomath.2014.04.212 Introduction THE CONTINUOUS FIRE MODEL OF ASYMMETRIC TREE-GRASS COMPETITION (COFAC) MATHEMATICAL ANALYSIS Existence of equilibria, ecological thresholds and stability analysis Existence of equilibria Ecological thresholds interpretation Stability analysis Summary table of the qualitative analysis NUMERICAL SIMULATIONS A nonstandard scheme for the COFAC NUMERICAL SIMULATIONS AND BIFURCATION PARAMETERS Some monostability and bistability situations Some bifurcation parameters Conclusion and discussion References Bibliographie