A tribute to the use of minimalistic spatially-implicit models of savanna vegetation dynamics to address broad spatial scales in spite of scarce data

—The savanna biome encompasses a variety of vegetation physiognomies that traduce complex dynamical responses of plants to the rainfall gradients leading from tropical forests to hot deserts. Such responses are shaped by interactions between woody and grassy plants that can be either direct, disturbance-mediated or both. There has been increasing evidence that several vegetation physiognomies, sometimes highly contrasted, may durably coexist under similar rainfall conditions suggesting multi-stability or at least not abrupt transitions. These fascinating questions have triggered burgeoning modelling efforts which have, however, not yet delivered an integrated picture liable to furnish sensible predictions of potential vegetation at broad scales. In this paper, we will recall the key ecological processes and resulting vegetation dynamics that models should take into account. We will also present the main modelling options present in the literature and advocate the use of minimalistic models, capturing only the essential processes while retaining sufﬁcient mathematical tractability and restricting themselves to a minimal set of parameters assessable from the overall literature.

Dynamics of vegetation within the savanna biome has long interested ecologists as it clearly departs from the classical post-disturbance succession pathways that are expected to rapidly bring back closed canopy forest, as observed in most of temperate and wet tropical climatic zones   [22]). The last decades have witnessed burgeoning efforts of modelling as to account for the possibly long-lasting coexistence of grassy and woody components and try to predict potential shifts from two-phased vegetation physiognomies. These efforts have, however, not yet delivered an integrated picture liable to furnish at broad scales (i.e., for fractions of continents) sensible predic-tions of possible vegetation dynamics. Such a big picture is nevertheless desirable for figuring out the future of vegetation in the face of climate and anthropic change scenarios (Mayaux et al. (2004) [60],   [22], Archibald et al. (2009) [7], Accatino et al. (2010) [4], Favier et al. (2012) [36]). It is also necessary for applications to territories devoid of reference data and longterm observation sites, as it is the case for most of tropical Africa.
The objectives of the present contribution are fourfold. It first aims at recalling and synthetizing the main array of facts about ecological processes and resulting vegetation dynamics that models should aim to capture and render (see section II). Second, in order to claim genericity, we synthetize the main modelling options present in the literature, and put emphasis on minimalistic models, capturing only essential processes while retaining sufficient mathematical tractability and restricting themselves to a minimal set of assessable parameters (see section III). Thirdly, on this basis, we argue that such models have now become more comprehensive, and useful for meaningful predictions (see section IV). Finally, we discuss how those models may now help guiding data collection for improved calibration and testing of dynamical hypotheses (see section V).

TREE-GRASS INTERACTIONS MODELLING
A. Tree-grass coexistence and possible alternative stable states Over very large tropical territories, field observers have documented long lasting coexistence of notable levels of grass and woody biomass (Backéus (1992) [9]). The most frequently reported form of coexistence is observed locally through vegetation physiognomies that associate fairly continuous grassy cover and more or less scattered populations of trees and shrubs of varying clumping levels. This is referred to as savanna physiognomy (see Figure 1). Such vegetation types mixing both lifeforms are manifold and progressively merge in space or through time without clear-cut boundaries (Torello- Raventos et al. (2013) [95]). Another modality of long lasting association between herbaceous and woody lifeforms occur at landscape scale under the form of mosaics featuring forests (usually closed canopy ones) and open savannas or grasslands (e.g. Figure  1; Bond and Parr (2010) [21]). In those landscapes that pertain to moist-wet climates, normally seen as favourable to forests, the mosaics appear highly contrasted and among the most "emblematic vegetation transitions" in the world (Oliveras and Malhi (2016) [72]): outside the closed forest, woody vegetation is of low biomass and the dominant physiognomies relate mainly to grassland. Moreover, boundaries between forest and grassland are generally sharp (Hoffmann et al. 2012 [48], Cuni-Sanchez et al. (2016) [27]).
Our interpretations of those various physiognomies are limited by the length of the observation windows we can rely on for distinguishing trends against fluctuations. For field observations, this window length barely extend over some decades and this only for a very small number of sites where invaluable data have been gathered. At the scale of extensive territories, representativeness of those sites remains yet an open question. Remote sensing is progressively broadening our observational means. But the best nowadays space-borne sensors for estimating woody cover (Buccini and Hanan (2007) [26]) or biomass (Mermoz et al. (2014) [64], Bouvet et al. (2018) [24]) are recent and do not allow tracking changes far back. Moreover, the accuracy of those estimations, notably for woody cover is limited, due to the difficulty to separate grass vs. tree in signal responses in mixed stands. This is particularly true regarding long diachronic series that mainly feature optical images of insufficient spatial resolution. Apart from blatant changes, e.g., forest encroachment or recession, Mitchard and Flintrop (2013) [68], subtle evolution of the grass-tree balance in mixed physiognomies are still beyond reach. Remote sensing, however, recently brought two interesting contributions to the savanna debate. First, broad scale assessment of woody cover at regional (Central Africa, Favier et al. (2012) [36]) to continental/global scales (Hirota et al. (2011) [47]) clearly showed that contrasted levels of cover can coexist under the same ranges of climatic conditions, making the existence of multi-stable states at least plausible. Second, in both Central and West Africa, comparison between ancient air photographs from the 50s and satellite images from the 80-90s frequently evidenced a progress of forest over savannas/grasslands in landscape featuring contrasted mosaics of the type exemplified in Figure 1 (Youta Happi (1998) [114], Mitchard et al. (2011) [69]).
Even though there is still no conclusive evidence that alternative stable states may exist within the savanna biome, models should be able to account for them as plausible outcomes of tree-grass interactions. The same applies to savanna physiognomies locally associating trees and grasses that may be seen as either stable or transient twophase states. Since those mixed physiognomies are observable at broad scale, there is no reason to a priori rule out that some observed mixtures may be stable under their local environmental context. Indeed, hypothesis testing is a fundamental role of models, though this use is not so widespread in ecology. And to this aim, the wider the array of reasonable predictions the more relevant is the model.

B. Lines of thoughts
Most authors agree on the fact that soil water budget, herbivory (i.e. grazing and/or browsing) and fires are the principal factors influencing growth of woody and herbaceous plants and their dynamical interactions (Scholes and Archer (1997) [51]). Authors however diverge on the relative importance of those factors in shaping dynamical outcomes of tree-grass interactions. This is not surprising considering the broad extent of the savanna biome Airborne photo from N. Barbier, June 2017. and the variety of both environmental conditions and anthropogenic pressures that apply therein. A factor appearing pervasive in a given context is not systematically due to prevail elsewhere. One group of authors has been insisting on direct interactions among or between plant-types (i.e. tree-tree or tree-grass) such as competition for light or for soil limiting resources (often moisture via root systems) (e.g. Scholes and Archer (1997) [79], Scholes (2003) [78]). It is obvious that the treegrass interaction is highly asymmetric: trees have a strong competitive effect on grasses, but grasses have a weak competitive effect on mature trees, although they may have a strong effect on saplings that have not grown above the grass layer (Scholes (2003) [78], Figure 2-a).
Another group of authors has been emphasizing that woody vegetation would be likely to reach a closed canopy situation and suppress grasses in the absence of recurrent disturbances induced by fires or browsers (or both sources) that delay or block the build-up of woody biomass by destroying the aerial part of seedlings and saplings (e.g.   [22], Bond (2008) [17] Literature may sometimes overemphasize the distinction between 'interaction' (between plant types for limited resource) and 'disturbance' hypotheses (see Scholes and Archer (1997) [79]) as to make them appear as alternatives, though they are by no means mutually exclusive. It is widely acknowledged that to have notable impact on vegetation, fire disturbance requests sufficient intensity through enough dry grass biomass as main source of fuel. Under a certain level of grass biomass, owing to insufficient rainfall or intense grazing, fires tend to spread difficultly and, where occurring, have modest impacts on woody plants. Logically, most authors tend now to distinguish disturbance-limited (i.e., under moist climate) vs. water limited (i.e., arid) savannas (e.g. Bond et al. (2003) [20]). Inter-tree competition shapes the second type (Sankaran et al. (2005) [75]), while asymmetric and fire-meditated tree-   [31]) or because of modulation by edaphic conditions, grazing and anthropogenic pressures. Grazing may lead savannas toward physiognomies and functioning looking less fire-prone, i.e. more "arid-like", than expected from the only climate features as an emergent consequence of dynamical amplification of external forcing.

III. MAIN PUBLISHED MODELLING OPTIONS
The questions raised by observed or putative dynamics within the savanna biome have triggered an increasing interest in terms of modelling. Pioneering works (Walker et al. (1981) [103], Walker and Noy-Meir (1982) [104]) first used systems of ordinary differential equations (ODE) to address the particular case of arid, fire-immune savannas in which excessive grazing fosters bush encroachment (Skarpe (1990) [82]). This line of modelling featured grass and woody biomasses as state variables and aimed at explicitly depicting their interactions in relation to soil moisture dynamics. As such, it became a paradigm for 'interaction models' involving a limited resource, but the central assumption of soil niche partitioning between the two plant forms called Walter's (1971) hypothesis [105] has been ever since hotly debated and is obviously not verified in all ecological contexts where savannas, dry thickets or grasslands are observable.
Another line of ODE-based modelling built on the application to savannas of the initial concept of asymmetric competition of (Tilman (1994) [94]) through a simple framework that allows considering both direct and disturbance-mediated plant interactions. Tilman's framework reinterpretation (see Accatino et al. (2010) [4]), De Michele et al. (2011) [28] used two states variables, namely cover-fractions of grass (G) and tree (T ) assumed exclusive and summing between zero and one. It modelled their interacting dynamics in a system of two ODE.
where, T and G are dimensionless and denote the fractions of sites occupied by tree and grass respectively. c T and c G are the colonisation rates of tree and grass respectively. δ T and δ G represent the mortality rates of tree and grass respectively. In the sequel, we refer to system (1) as Tilman's model. Logistic growth of the inferior competitor (grasses plus herbs) is bounded and depressed by the cover of the superior competitor (woody plants) which logistic growth is not directly affected by grasses (asymmetric competition). In system (1) where, δ F represents the trees vulnerability to fire, f is the fire frequency (inversely proportional to fire return time period) and ω(G) is a function of grass biomass that represents the fire impact.
Through ω(G), there is thus indirect, fire-mediated negative feed-back of grass cover onto tree cover that counterbalance direct, tree-grass asymmetric interactions. A larger array of models (see Tables I and II) took a leaf from the previous modelling framework (system 1 and equation 2). Main sources of variations between models were: (1) nature of the equations and temporal treatment of fire disturbance (time-continuous forcing, i.e. ODE vs. time-discrete or impulsive occurrences); (2) nature of the function expressing grass-fire feedback on trees (linear vs. nonlinear); (3) integration of herbivory in addition to fire; (4) facultative explicit treatment of water availability through models with one (and sometimes more) additional state variables expressing water resource in interaction with vegetation variables. We will refer to such models as 'ecohydrological' (see Table I), among which is system (3) proposed by Accatino et al. (2010) [4] that features a first equation devoted to the dynamics of a soil moisture variable (S) where p w 1 (per year) represents the rainfall rate normalized with respect to root zone capacity, ε (per year) is the evaporation, τ T and τ G (per year), are water uptake parameters for tree and grass respectively. c T , c G , δ T , δ G , δ F and f are defined as in system (1) and equation (2). Note that setting ω(G) = 0 in the second equation of (3) makes Tilman's model (1) analogous to the system coupling the second and the third equations of (3). Moreover, if the S variable is held constant, the main difference between systems (1) and (3) is that Accatino et al. (2010) [4] considered ω(G) = G (i.e., impact of fire on trees as a linear function of grass biomass) while in Tilman's model, this function is equal to zero (no impact of fire).
Taking ω(G) as any increasing function of the grass cover, provides a more general expression of the fire impact on trees. Without loss of generality, we referred to Holling type functions (Holling (1959) where, G in tons per hectare (t.ha −1 ) is grass biomass, α is the value takes by G when fire intensity is half its maximum, and the integer θ determines the steepness of the sigmoid. Nonlinear response was retained by some other authors Considering a nonlinear (sigmoidal) shape for ω(G) allows for the existence of up to three treegrass coexistence (i.e. savanna) equilibria, while two of them may be simultaneously stable (i.e. bistability) and forest-savanna-grassland tristability is reachable (Yatat et al. (2014) [112], (2018) [109], Tchuinté Tamen et al. (2014) [91], (2017b) [88]). Conversely, we proved that for linear ω(G) functions, tristability is unreachable (Yatat et al.  [36] obtained along a general climatic transect over central African (latitude in the range of 3 -4 • north). These results concerned a very large range of woody cover (wc) variations (from very low values approaching grassland physiognomies to nearly 80% cover (i.e. forest) through wc of 40%, i.e. savanna) which suggests grassland/savanna/forest tristability as at least plausible. ODE models have been criticized on the basis that they only predict abrupt transitions between two alternative stable states (Accatino and De Michele (2016) [3]) that are deemed unrealistic. However, tristability of equilibria as well as bistability of two savanna equilibria suggests that shifts from one stable state to another may be less spectacular than hypothesized from previous models and that models may render more complex pathways of vegetation changes (see   [109] for further discussion).  [85]) that considered cover fractions. Modelling biomasses help accounting from the fact that plant types are not mutually exclusive at a given point in space since field studies suggested that grass often develop under tree crowns (Belsky et al. (1989) [16], Belsky (1994) [15], Weltzin and Coughenour (1990) [96]). Some models used more than two size classes through matrix population models (e.g. Accatino et al. (2013) [2], (2016) [5] and references therein). Simpler models used only two sizerelated variables in addition to grass and simply account for the asymmetric nature of tree-grass interactions as discussed previously (see also Scholes (2003) [112], (2017) [110]). This distinction stems from field observations (Trollope (1984) [97], Trollope and Trollope (1996) [99]) that evidence rapid decline of percent topkill with tree height (see Table I: Comparison of several models of tree-grass dynamics with respect to some modelling options. Walter's hypothesis refers to the differences in root depth of herbaceous and woody vegetation in water seeking while ecohydrological frameworks stand for models that consider additional state variables expressing water resource in interaction with vegetation variables. From Tchuinté Tamen (2017) [87] and Yatat (2018) [113]. The symbol * means that we refer to system (1 [112], (2017) [110]) as long as other complexities were not introduced.
A strong objection against ODE models is that fire is not a forcing that continuously removes a small fraction of biomass through time as per the previous ODE equation systems. Instead, fire actually suppresses a substantial fraction of biomass at once through punctual outbreaks that shape ecosystem aspect and immediate post-fire functioning ( Figure 4). This principle was implemented Table II: Summary of the characteristics of tree-grass interactions models with respect to fire modelling options (continued Table I). From Tchuinté Tamen (2017) [87] and Yatat (2018) [113]. Some references (unticked) do not model fire. The symbol * means that we refer to system (1).     [88], see also Figure 7) based on two important parameters of strong intuitive meaning, namely mean total annual rainfall and fire frequency. They thereby achieved delineation of domains in which main physiognomies (i.e., grassland, savanna, forest) can be expected as stable. Bistability situations (forest-grassland and forest-savanna) were also highlighted for sufficiently high fire frequencies. For low fire frequencies, a sensible gradient of increasing woody cover with increasing annual precipitation was found. But one may note that only dense woodlands (i.e. still two-phase vegetation) were obtained even for the highest rainfall range, while forest stricto -sensu (mono-phase, with no tall, light-demanding savanna grasses in the understory) is widely observed for the corresponding ranges of precipitations. Moreover, transitions between vegetation types in relation to fire frequencies proved tricky. Indeed, in the high rainfall range, increasing fire frequency leads from the aforementioned woodlands to forest. In the intermediate rainfall range, increasing fire frequency leads from savanna mono-stability to forest-savanna bistability (see Figure 7). Analogously, for fairly low rainfall, grassland stability shifts to grassland-forest bistability. Hence, all over the rainfall gradient it looks as if sufficient frequency of fire were a necessary condition to reach forest (bi)-stability. This is contradicted by empirical knowledge according to which frequent fires are known to jeopardize or at least delay woody biomass build-up, but never favour it   Where did this critical problem come from? Most of subsequent papers barely evoked the question. A large share of them investigated different modelling options, often more complex and/or less tractable; or they assumed particular biogeographic conditions. In a further contribution, Accatino and De Michele (2016) [3] argued about intrinsic limitations of ODE-based modelling. They also put forward that there is no evidence according to which observed vegetation physiognomies may be close to a stable equilibrium point. It is in fact undisputable that climate is likely to vary through time, and there is no guaranty that woody vegetation can track such variation with enough celerity. They also underline as questionable the assumption according to which the parameter f of fire frequency should be treated as a constant forcing, independently of vegetation characteristics. All these arguments brought them to propose a 'non-equilibrium model', based on stochastic difference equations, as alternative to the timecontinuous model of Accatino et al. (2010) [4] referred to as 'equilibrium-model' (EM). In their non-equilibrium (NEM) model, fire occurrence is a stochastic event all the more likely to occur in a given dry season that ignitable dry grass biomass abundantly built-up in the foregone rainy seasons.
Accatino and De Michele (2016) [3] compared the predictions of their (EM) vs. (NEM) models. They argued that separation of fire-immune vs. fire-prone savannas as an indirect consequence of rainfall is an emergent property with their NEM while it is artificially induced by the choice of the f parameter with EM. They also pointed out that when considering high rainfall, the NEM is able to reproduce the "bimodality" of woody cover extent observed in remote sensing studies. But in fact, their NEM was not a straightforward time-discrete analogue of the ODE based EM of Accatino et al. (2010) [4], since it features several novelties. Therefore, the differences they reported between EM and NEM are not a simple consequence of time-continuous fire forcing vs. time-discrete fire occurrences. In fact, several aspects altogether contribute to the more satisfactory results obtained with the NEM by Accatino and De Michele (2016) [3]. We will subsequently illustrate the fact that predictions that are qualitatively satisfactory can be obtained by directly improving the ODE based 'EM' framework, notably regarding the firemediated feedback of grass onto tree dynamics.

B. Fire frequency, grass biomass and fire impact
The way in which fire impact is modelled in the previous equation systems (1 and 3) is obviously a crucial question. As stated by Scholes (2003) [88] [78], modelling savanna dynamics in fire-prone contexts actually requires introducing an "equation predicting the effect of grass biomass, via fire intensity, on tree biomass".
In fact, non-linearity in ω(G) may be justified on various, non-exclusive grounds, since what is important is to properly model, as a whole, the causal chain that leads from grass abundance and ignition regime to woody biomass suppression. As steps in this chain we may identify: (i) fire frequency to be seen as an external forcing upon the tree-grass system (think about a targeted fire regime in a managed area such as a ranch or a protected area); (ii) actual yearly fire probability (or frequency) of occurrence in any arbitrary small piece of land once (i) has been set; (iii) fire potential impact (intensity and flame height) on woody vegetation; (iv) fire actual impact that also depends on features intrinsic to woody vegetation (see below section IV-C).
Fire intensity which is strictly speaking a quantity of energy released (Bond and Keeley (2005)  In our modelling, the f parameter was kept as constant multiplier of w(G), but we interpret it as a man-induced "targeted" fire frequency (as for instance in a fire management plan), which will not translate into actual frequency of fires of notable impact as long as grass biomass is not of sufficient quantity. With this interpretation, the actual fire regime may substantially differ from the targeted one, as frequently observed in the field. And, whatever f values, any hypothetical piece of land will actually be fire-prone only if other forcing factors (climate, herbivory, etc.) allow for sufficient grass biomass. Most previous modelling papers including those from our group did not elaborate much regarding the successive steps involved in the grass-fire feedback. Distinction between (i) and (ii) may appear subtle and to our knowledge was never emphasized before. It directly results from space-implicit savanna modelling. Fire regime, which is nowadays overwhelmingly maninduced (Govender et al (2006) [44], Archibald (2009) [7]) is a forcing at landscape scale since people do not go and set fire in every piece of land. They instead count on fire propagation that depends on abundance and spatial evenness of dry grass. In presence of low and unevenly distributed grassy fuel, fire will barely propagate leaving a large share of the area unburnt. This makes the difference between steps (i) and (ii), as frequently observed in the field (see Diouf et al. (2012) [31]; Figure 6). The response of percent area burnt to grass abundance is likely to be sharply nonlinear, as suggested by the impressive results reported by McNaughton (1992) [61] at the scale of the entire Serengeti complex in Tanzania. In this remarkable study, local fire frequency dwindled over a decade following grass biomass suppression by soaring herbivore populations, while the ignition regime by people dwelling around the park likely remained more or less the same. In fact, since we are here dealing with mean-field models the ω(G) function is also due to embody the difficult spreading of fire in presence of fuel of overall low quantity keeping in mind natural spatial variability of grass biomass ( Figure 6). Non-linearity of ω(G) seems therefore a necessary feature for adequately capturing the fire-mediated grass-tree feedback.

C. Tree survival
To be relevant, the most parsimonious models featuring just grass and tree state variables must overcome the limitation pointed out in sub-section IV-A for the precursor model of Accatino et al. (2010) [4]. All things being equal, any increase in fire frequency should never increase the woody component of vegetation. Fire, if any, is expected to be of no substantial consequence over the driest stretch of the rainfall gradient while for the moister part, it is widely observed that extending the average time between successive fires (decreasing frequency) favours the building up of woody vegetation. Accounting for that proved to be a challenge for minimal two-variable models that do not distinguish between fire sensitive and fire insensitive woody fractions. Non-linearity of the ω(G) function, though important proved not sufficient to overcome this problem. Tchuinté Tamen et al. (2017) [90] further introduced a second non-linear decreasing function, which directly expresses that high woody biomasses, corresponding to tall trees proportionately experience far less firerelated losses than low woody biomasses relating to seedlings, saplings and shrubs (see Figure 2). We hence proposed the following function to denote the fire-induced tree/shrub mortality:  [90] designed a model that also implemented punctual fire events through impulsive differential equations. But predictions of the ODE model itself were already satisfactory.

D. Relating to water resource
Several savanna models have explicitly modelled soil water resources, via a dedicated equation as in system (3). But the soil moisture dynamics is very rapid compared to change in vegetation. Soil moisture variations linked to a given rainfall event are damped within a few days (Barbier et al. (2008) [10]), while vegetation growth proceeds over months for grasses and even years for woody plants. It therefore makes sense to consider vegetation dynamics in relation to a level of soil water resource that is approximately constant for a given level of total annual rainfall or, ideally, water deficit (rainfall minus evapotranspiration). Martinez-Garcia et al. (2014) [59] reached similar conclusions for a partial differential equation model of non-local plant-plant interactions for water. This justifies letting parameters in the vegetation dynamics equations directly depend on climatic parameters. This principle sustains the model expressed through the following system: where, • G and T are grass and tree biomasses respectively; • W is the mean annual precipitation (in millimeters per year, mm.yr −1 ); are annual productions of grass and tree biomasses respectively, where γ G and γ T (in yr −1 ) express maximal growths of grass and tree biomasses respectively, half saturations b G and b T (in mm.yr −1 ) determine how quickly growth increase with water availability; • K G (W ) = c G 1 + d G e −aGW , and K T (W ) = c T 1 + d T e −aT W are carrying capacities of grass and tree respectively, where c G and c T (in t.ha −1 ) denote the maximum values of the grass and tree biomasses, a G and a T (mm −1 yr) control the steepness of the curves of K G and K T respectively, and d G and d T control the location of their inflection points; • δ G and δ T respectively express the rates of grass and tree biomasses loss by herbivores (grazing and/or browsing) or by human action; • η T G denotes the asymmetric influence of trees on grass for light (shading) and resources (water, nutrients) which relate to competitive or facilitative influences; • λ f G is the specific loss of grass biomass due to fire; • f = 1 τ is the fire frequency, where τ is the fire return period.
Submitting model (6) to bifurcation analysis provides Figure 7-(b) that is to be compared to Figure 7-(a) from Accatino et al. (2010) [4], which has been reobtained using Matcont (see Govaerts (2000) [42], Dhooge et al. (2003) [29], Govaerts et al. (2007) [43] and references therein). Both Figures 7-(a) and 7-(b) are sensible regarding low fire frequencies for which increasing MAR leads to a sequence of physiognomies of increasing woody biomass (i.e., grassland, savanna, forest). But in Figure 7-(b), the improvement resulting from introducing ω(G) and ϑ(T ) is apparent when increasing fire frequency (f ) at different levels of the rainfall gradient. For high MAR values, the expected physiognomy shifts from monostable forest to forest-grassland bistability. Indeed, in presence of high MAR, it is known that grasslands are due to be encroached by forest under fire prevention or even just because of decreasing fire frequencies (Jeffery et al. (2014) [51]). Our model accords with field observations in that a high fire frequency is indeed a necessary condition to perpetuate the grassland (or savannas of low woody biomass) physiognomies. Moreover, largescale observations of bimodality between high and very low woody cover situations (Hirota et al. (2011) [47], Favier et al. (2012) [36]) can be accounted for by the forest-grassland bistability, though the converse is not necessarily true. In fact bimodality may stem from either transient situations or topographical heterogeneity and does not automatically implies bistability. For low to intermediate MAR values, say 600 − 1000 mm, fire is known to be less pervasive, though field observations or experiment results depict woody vegetation thickening in case of fire frequency decrease (Brookman-Amissah (1980) [25]). The model is able to render such thickening as a shift from the grassland to the savanna stability domain. The model also predicts forest-savanna or savanna-grassland bistability, and even tristability thereof for restricted domains in the MAR-fire frequency plane that were situated around 1000 mm of MAR. This value is of course dependent on parameter values used for computations underlying figure 7-(b). Refined calibrations relating to a specific regional context may displace the thresholds. Notably, the parameter expressing the influence of woody biomass on grasses proved to be influential on the thresholds between vegetation states (Tchuinté Tamen et al. (2018) [88]).

E. Impulsive time-periodic occurrences of fire events
In previous models, the traditional timecontinuous fire forcing formalism is often used. However, it is questionable to model fire as a permanent forcing that continuously removes fractions of fire sensitive biomass all over the year. Indeed, several months and even years can pass between two successive fires, such that fire may be considered as an instantaneous perturbation of the savanna ecosystem (Yatat et al. (2017) [110], Tchuinté Tamen (2016) [89], (2017) [90]; see also Table IV, [13]) is that they barely lend themselves to analytical approaches.
Based on Table IV, page 18, we further consider in our group (Yatat et al. (2017) [110], (2018) [109],   [111], Tchuinté Tamen (2016) [89], (2017) [90], (2018) [88]) impulsive time-periodic fire events which is an approximation that keeps the potential of analytical investigation as large as possible while modelling discrete fires. An impulsive differential equations system can be used to express fire through impul-sive periodic occurrences (e.g. system (7), below): ∆T (nτ ) = −ϑ(T (nτ ))ω(G(nτ ))T (nτ ), where, • for π ∈ {G, T }, ∆π(nτ ) = π(nτ + ) − π(nτ ) and π(nτ + ) = lim θ→0 + π(nτ + θ); • τ = 1 f is the period between two consecutive fires; • N f is a countable number of fire occurrences; • nτ , n = 1, 2, ..., N f , are called moments of impulsive effects of fire, and satisfy 0 ≤ τ < 2τ < ... < N f τ . Properties of models (6) and (7) have been analysed in Tchuinté Tamen et al. (2018) [88]. Below we provide some numerical simulations as to illustrate the bifurcations between the stable domains delineated in Figure 7-b as a consequence of increasing the fire frequency for two particular values of MAR, i.e. W = 946 mm.y −1 and W = 1003 mm.y −1 . For each of the two MAR values, we compare the consequences of increasing f with the ODE (system (6)) and the IDE (system (7)) frameworks, by comparing Figure 8 against Figure 9 and Figure 10 against Figure 11. For the lower MAR situation, ODE and IDE frameworks qualitatively agree and show the bifurcation from monostable savanna to monostable grassland through an intermediate bistable situation (Figure  8-b; Figure 9-b). For the higher MAR case, both frameworks also show the transition from monostable forest to forest-grassland bistability through tristability involving savanna. Qualitatively, the   [41] predictions of the two frameworks agree about the predicted sequence of vegetation physiognomies when increasing the fire frequency. However, the IDE model systematically predicted bifurcations for lower values of f than for the ODE. This indicates that shifting to the conceptually more satisfactory IDE framework will introduce specificities in forthcoming stages concerning refined calibration and comparison with real-world observations.

V. DISCUSSION AND PROSPECTS
In the present paper we emphasize that ecologists did probably not yet exploit all the potential of simple ODE systems for modelling vegetation dynamics in the savanna biome to which most seasonal tropical ecosystems pertain. We showed that reasonable, non-trivial predictions can be obtained in reference to hypothetical situations directly relating to rainfall and fire frequency gradients. Application to specific contexts and locations would request refined calibration for the parameters of the generic minimalistic model. But it may also invite to better address specific processes deemed influential in a particular situation under study. This would mean complexifying the model to match a specific piece of reality. Though this is actually a natural and sensible trend in science, parsimony is an opposing principle that tells us to keep complexification under control. A meaningful, balanced modelling approach should restrict to what we strictly need to account for a well-defined array of empirical facts in a particular situation.
On the empirical side, the ongoing development of remote-sensing techniques and derived products is providing avenues to better depict the spatiotemporal variation of environmental factors, such as rainfall (e.g. via CHELSA, Karger et al. (2017) [52]), topography (via the SRTM, Farr et al. (2007) [35], published at increased spatial resolution by NASA in 2013) or fires (http://modis-fire.umd.edu; e.g. Archibald et al. (2009) [7], Diouf et al. (2012) [31]). This also applies to the monitoring of some vegetation variables, though disentangling effects on most remotely-sensed signals from grasses vs. trees in mixed savanna physiognomies is still challenging. Improving and diversifying sources of remote sensing information will obviously help sorting out relevant predictions from unrealistic ones and refine the benchmarking of models. But most of the parameters expressing vegetation dynamics will remain out of reach of remotelysensed assessment and will remain dependent on field information. An increased effort of field data and f = 0.7, respectively when using the ODE model (system (6)). collection is obviously desirable but insufficient means for research in most tropical countries is enduring reality. That strong data limitation is alas probably here to stay finally pleads for parsimony in modelling. It also underlines the importance for modelling to be sufficiently convincing and accessible to ecologists as to guide data acquisition and orient scarce resources towards assessing parameters proven as the most influential by sensibility analyses.
On the modelling side, within the class of models that distinguished size classes for the woody component of vegetation, some used more than two size classes through matrix population models (e.g. Accatino et al. (2013) [2], (2016) [5]). However, such models remain generally simulationbased and usually involve a fairly large number of parameters. Thus, it is not easy/possible to assess how model parameter variations may influence the model outcomes   [109]). ODE and IDE models featuring two woody variables, in addition to grass, proved analytically tractable (Beckage et al. (2009) [14], Staver et al. (2011) [83], Yatat et al. (2014) [112], (2017) [110]) as long as other complexities were not introduced. For example, Yatat et al. (2014) [112] (resp. Yatat et al. (2017) [110]) studied a ODE-like (resp. IDE-like) tree-grass interactions model where in addition to grass they considered two classes of woody plants: fire-sensitive like seedlings and fire insensitive. But, based on recent publications of our group, we found that even with models that feature only one state variable for woody component, meaningful results are obtained and, to some extent, are qualitatively similar to those obtained with models that used two size-related variables for woody component (Tchuinté Tamen et al. (2014) [91], (2016) [89], (2017) [90], (2018) [88],   [109]).
The IDE framework is an obvious improvement that expresses a reasonable trade-off between increased realism and decreased analytical tractability. Within the framework of IDE, a further step in that direction could be considering patterns of fire occurrences featuring stochastic components instead of the deterministic periodic regime we used (Yatat et al. (2017) [110], (2018) [109], Tchuinté Tamen et al. (2016) [89], (2017) [90], (2018) [88]). But one may note here that periodicity of fire outbreaks is not that unrealistic in the context of subequatorial humid savannas for which fires can only occur at the end of dry seasons which are of short duration. Here the annual climatic cycle strongly defines the temporal window for fires while in less humid savannas fire is simultaneously less frequent and liable to occur all over extensive dry seasons (Diouf et al. (2012) [31]). We believe that while mathematical tractability or theoretical study of a model is not an absolute requirement or is not always possible, it remains nevertheless desirable at least for two reasons. First, it can appear as a kind of guarantee that numerical simulations displayed by the model are not the result of some numerical artifacts. In other words, the choice or the construction of a suitable algorithm to solve a given (complex) model strongly relies on its qualitative study or, when this study is not possible, on the analysis of some sub-models, that can be mathematically tractable. Nowadays, there are more and more works that point out some spurious behaviors that may appear when using some 'classical' schemes for model simulations (see for example Anguelov et al. (2012) [6]). Second, any theoretical analysis may provide useful informations about the role of some particular parameters in the dynamics of the system.
We have here focused on spatially-implicit models because we believe that such models have still important insights to provide and also because spatially-explicit models are far more demanding in terms of parametrization and more difficult to study theoretically. Substantial efforts to design and run spatial models of savannas have however been made during the last decade (Borgogno et al. (2009) [23]). Most of them relied on individual-based models such as cellular automata. At this step of the discussion, it seems meaningful to point out that there are some authors who pleaded for a mutualistic or complementary relationship between mathematical tractabilitybased and simulation-based formalisms (Omohundro (1984) [73], Wolfram (1985) [108], Weimar (1997) [106], Narbel (2006) [71], Dietrich et al. (2014) [30],   [34]) and we also agree with that. As an illustration, it may be more difficult to achieve a very deep theoretical analysis of a partial differential equations model when taking into account spatial heterogeneity while it seems more easy to handle it when using for example an individual-based model or cellular automaton formalism. Independently, partial differential equations (PDE) have also widely been used to account for the genesis of vegetation regular spatial patterns (bare soil vs. dense shrubby cover) in the particular context of arid savannas (Lefever and Lejeune (1997) [56], Klausmeier (1999) [53], Sherratt (2005) [ [90].
The analysis of this model shows that there exists monostable or bistable travelling solutions, related to the boundary movements in the forestgrassland mosaic. And the authors showed that depending on fire-return time as well as difference in diffusion coefficients of woody and herbaceous vegetation, fire events are able to greatly slow down or even stop the progression of forest in the humid part of the savanna biome. This kind of results, obtained from theoretical analysis, are of great interest for practical needs or management policies. However, as an improvement of   [109], there are also some ongoing works that aim to deal with existence of travelling waves for system of Impulsive Partial Differential Equations (IPDE) that model savanna dynamics (e.g.   [111]). This type of modelling is of great interest considering that the forest-savanna ecotone is the most widespread in the tropics and that forests have been encroaching during the last decades in many humid savannas of West and Central Africa, and to a lesser extent in other regions of the world (Oliveras and Malhi (2016) [72]). On the other hand, forest encroachment, which is of great consequence for the global carbon balance of terrestrial ecosystems proved heterogeneous in space and time (Oliveras and Malhi (2016) [72]). Hence, modelling should help better understand the hierarchy of processes and forces accounting for such heterogeneity. This makes spatially-explicit modelling desirable. But, interpretation and calibration of local, diffusion operators as used in   [109] cannot rely on much empirical knowledge in plant ecology. And one may note that this also applies to colonization rates between adjacent cells that are central to cellular automata models. In the case of PDEs, modelling non-local plant-plant interaction processes (e.g. through interaction kernels) is mathematically more challenging, though attempts in Martinez-Garcia et al. (2014) [59] and in Lefever et al. (2009) [55] or Lefever and Turner (2012) [57] nevertheless provide sources of inspiration.
A recent line of criticism questions the relevance of reasoning in reference to equilibrium states. There is indeed no compelling evidence that physiognomies presently observable correspond or even are close to predictable equilibria determined by current environmental conditions. In fact, parameters reflecting environmental variables, notably climate are due to fluctuate or even change through time. And vegetation, especially its woody component may be unable to keep pace with such variations and rather track them with delay thereby remaining distant from any equilibrium state. Long-lasting consequences of past climate periods probably still mark present vegetation. For instance, in Central Africa, concomitant to, drier period occurred some centuries ago (Europe's "Little Ice Age", 500-200 years BP), which probably provoked forest cover recession and fragmentation (Oslisly et al. (2013) [74]). The trend of widespread forest boundaries displacement within savannas, as observed during the last decades may be a delayed recovery after this past drier period that is progressing at slow and unequal pace owing to the counteracting influence of fires. Increased CO 2 availability that favours C3 woody plants against C4 savanna grasses (Bond and Midgley (2000) [19]) may also reinforce this trend.

VI. CONCLUSION
Minimal savanna models using ODE systems have been criticized from different standpoints. The first one was the poor realism of the overall picture made by the predicted stable equilibria. The present paper shows that some unrealistic predictions are not a direct drawback of the ODE framework, but rather derive from inadequate modelling of the crucial fire-mediated negative feedback of grassy biomass onto woody vegetation. Using nonlinear functions such as ω(G) and ϑ(T ) (as in equations 5) is not only justified by the nature of the mechanisms at play, but also proved sufficient to get a meaningful "big picture" of vegetation physiognomies predicted as stable for varying mean annual rainfall and fire frequency (as in Figure 7-(b)). Another argument against ODE models is that they predict too contrasted stable states meaning that shifts between them would appear as more catastrophic than actually observed. But some strong contrasts such as landscape mosaics of forest and grassland are indeed observable in the field (see Figure 1). Moreover equilibria that have attraction domains "adjacent" in Figure 7 do not systematically show contrasted biomass values. Transitions may actually be progressive in terms of state variables. Indeed, at low fire frequency, the transition from savanna to forest along the MAR gradient corresponds to a continuous increase of tree biomass with concomitant decrease of grass biomass. The same applies to the transition from savanna to grassland via increased fire frequency (see Figures 8 and 9). Here, the shape of the nonlinear functions embodying the grass-fire feedback matters as shown in previous works (Yatat et al. (2014) [112], Tchuinté Tamen et al. (2014) [91]). Strongly nonlinear shapes (e.g. Holling functions of higher order see Table III) allow for multiple coexistence equilibria (i.e. multiple savanna physiognomies) of different woody biomass values, which may be seen as "stepping stones" between grassland and forest. Thus, the ODE framework does not automatically imply abrupt changes between equilibria of very contrasted biomass values.
On the other hand, it is undisputable that modelling fire as an external forcing continuously suppressing small amount of biomass through time is not satisfactory. Models based on punctual fires impacting large shares of biomasses are more relevant. This is implemented in time-discrete stochastic models which are however of limited analytical tractability. Impulsive differential equation systems are a good compromise since they permit to model time-discrete fire impact while keeping ODE for modelling vegetation growth. As such, they remain analytically tractable to a large extent while being more realistic.
In this paper, we show that minimal savanna models are able to provide a wide array of meaningful and relevant predictions of savanna dynamics while retaining sufficient mathematical tractability and restricting themselves to a minimal set of parameters assessable from the overall literature. Moreover, simplicity is overarching whatever the level of tractability. With a simple model, simulations can claim a thorough exploring of all parameters space. Conversely, it is difficult to be sure that sufficient exploration of model behaviors has been carried out for overcomplicated models which tend to flourish in ecology. Moreover, because there is naturally substantial uncertainty for many parameter values, it is difficult to conclude whether results are due to the ranges taken for parameters or to the structure of the model itself. Therefore, using complex models, it becomes even more illusory to test hypotheses, while this is one of the fundamental roles of modelling.