Pattern Formation in Mesic Savannas

We analyze a spatially extended version of a well-known model of forest-savanna dynamics, which presents as a system of nonlinear partial integro-differential equations, and study necessary conditions for pattern-forming bifurcations. Homogeneous solutions dominate the dynamics of the standard forest-savanna model, regardless of the length scales of the various spatial processes considered. However, several different pattern-forming scenarios are possible upon including spatial resource limitation, such as competition for water, soil nutrients, or herbivory effects. Using numerical simulations and continuation, we study the nature of the resulting patterns as a function of system parameters and length scales, uncovering subcritical pattern-forming bifurcations and observing significant regions of multistability for realistic parameter regimes. Finally, we discuss our results in the context of extant savanna-forest modeling efforts and highlight ongoing challenges in building a unifying mathematical model for savannas across different rainfall levels.


Introduction
Savanna and tropical forest are two of the most crucial biomes with regard to current conservation efforts, with savannas covering approximately one eight of the Earths total landmass and tropical forests a key carbon sink.Savannas are typically defined by coexistence of grass and trees with an open canopy.Their emergence has proven notoriously difficult to model and predict and relies on multiple complex processes, including fire disturbance and herbivory [49].Savannas are also diverse in nature; they exist over a wide range of mean annual rainfall (MAR), considered to be a strong determinant of vegetative cover, with the driest savanna found below 1000mm MAR and the wettest occurring above 2000mm MAR.Recent empirical findings suggest that savanna and tropical forest constitute stable alternatives at intermediate rainfall levels [37,39,40].This bistability suggests the possibility of switching or even irreversibly tipping between these states.Hence there is intense interest in developing appropriate models to capture this bistability, and thereby properly evaluate the stability and resilience to perturbations of tropical forest-savanna ecosystems.
There a burgeoning literature regarding the mathematical modeling of savannas and adjacent biomes such as woodland and tropical forest.Typically, these models involve interactions between one or more species of trees, shrubs and grasses, with each functional type of vegetation having different tolerances or reactions to fire and herbivory [41].Some of these models of savannas are spatially implicit and describe only proportions of a landscape occupied by each cover type [19,42].These models are accurate for describing well-mixed regions of vegetation, which is typically true at small length scales.However, when considering larger-scale systems that may include multiple types of vegetation, frameworks with explicit spatial dimensions are required.In recent years, such models were proposed either relying on partial differential equations (PDEs) [10,52], cellular automata [18,53] or spatial stochastic models [12,13].Much of this work has focused on the ranges and mechanisms for coexistence or multistability between biomes such as grasslands, savanna, woodland or forest in the tropics.
In spite of the significant modeling effort outlined above, there has been somewhat less attention paid to the spatial structure within the ecosystems predicted by these mathematical models, at least for mesic savannas, i.e. those at intermediate to high rainfall.In contrast, there is an extensive literature regarding pattern formation in vegetation models of dryland and semi-arid ecosystems, with the tiger bush model being among the best known examples [22].This work principally considers low rainfall ecosystems (less than 500mm MAR) in which water limitations highly constrain vegetative growth, and hence the models necessarily involve explicitly hydrological modeling of ground and soil water content, in addition to the vegetation state variables [20,28,50,36].Despite their more arid setting, the majority of these models share the qualitative feature of bistability with forest-savanna models, although in this case the bistability is between a state with non-zero vegetative biomass and a bare ground state.One of the key lessons from this body of work is that spatial patterning can have significant implications for the stability of ecosystems and may even qualitatively change the types of transitions that occur between alternative stable states.In particular, pattern-forming instabilities in a homogeneous stable state may not be a precursor to ecosystem collapse, but may instead allow the system to maintain a state of higher vegetative biomass for a much wider region of parameter space [32].It is evidently important to understand the extent to which these conclusions from pattern formation in arid ecosystems carry over to forest-savanna systems at intermediate rainfall, which have very different dominant interactions with different (relative) spatial scales.
Various empirical studies have reported evidence for spatial patterning in savanna ecosystems, particularly in nutrient deficient savannas [23] and at the boundaries between savanna and forest systems [33].This work inspired several branches of mathematical research seeking to explain these observations with different patterning mechanisms proposed for dry versus wetter savannas.Many models of spatial patterning at lower MAR are adapted from the arid ecosystem modeling paradigms and focus on patterns emerging from water resource competition and scarcity [4,17,48].These models typically involve systems of PDEs (sometimes with nonlocal operators for seed dispersal) with state variables for woody and herbaceous cover, in addition to soil and groundwater quantities.They produce a wide array of patterned steady states via Turing-type bifurcations.In addition, numerical continuation has proven an invaluable tool in this work, allowing detailed descriptions of the dynamics away from local bifurcations and thus highlighting possible paths for regime shifts.However, these works typically do not emphasize fire or herbivory effects, which are likely to be more significant at intermediate to higher MAR, and appear more suited to study ecosystem transitions and stability at the boundary between savanna and semi arid states with low vegetative cover.Other researchers have proposed spatially explicit mathematical models of mesic savannas focused on disturbance and facilitative interactions appropriate for higher MAR settings that can similarly produce patterned steady states [23,25,43].These works serve as an important proof of principle in terms of the mechanisms which can lead to stable heterogeneous savanna ecosystems, but have certain drawbacks that limit their applicability as general models for studying transitions between savanna and forest states in the bistable MAR range.For instance, some of these models are posed with a single state variable for biomass [23,25], affording mathematical tractability, but limiting their ability to reflect the functional properties of ecosystems by discriminating between different cover types.Other models require quite specific spatial kernels to establish stable patterns, raising questions about generality and robustness [43].
In this work, we aim at an intermediate complexity description of mesic savannas using a spatially explicit nonlocal mathematical model based on a well-studied forest-savanna framework often referred to as the Staver-Levin model [13,42].We demonstrate that stable spatially patterned ecosystems are possible in this framework under realistic assumptions on the lengthscales of the dominant spatial interactions at play.In order to retain a relatively tractable mathematical model, we don't explicit add hydrological constraints to our model, but we do include a nonlocal resource limitation on forest tree growth, in addition to the spatial fire spread and seed dispersal processes, since this turns out to be crucial for the formation of stable heterogeneous solutions.A particular advantage of our framework is that we can capture transitions between different biomes in much greater detail than previous modeling efforts in this domain.For example, our model draws clear distinctions between grassland, woodland, savanna and closed canopy tropical forest ecosystems, making it appropriate for more detailed future studies of ecosystem resilience and transitions between these, potentially multistable, alternative states.Moreover, we model spatial interactions via nonlocal operators, which is considered particularly appropriate for potentially long-range processes such as seed dispersal [29,45], affording us considerable flexibility to adapt our model to different scenarios and empirical data, as well as giving us robust and general qualitative conclusions.
The paper is organized as follows: In Section 2, we introduce a spatially extended version of the Staver-Levin model with nonlocal spatial interactions between the various vegetation types.We then perform a linear stability analysis to show that the steady-states of the original nonspatial model remain stable in the spatially extended model under very mild assumptions.In Section 3, we introduce a new version of the Staver-Levin model with resource limitation on the growth of forest trees and show that in this more realistic model homogeneous steady-state solutions may lose stability via Turing bifurcations, leading to the emergence of heterogeneous steady states.We investigate the structure of the bifurcations and the emerging heterogeneous solutions in detail for the grass-forest model as a function of system parameters and the different length scales in the model.We further show that in the full resource limited spatial Staver-Levin model, linear stability analysis once more confirms the presence of Turing bifurcations and the emergence of heterogeneous steady states is confirmed by numerical simulations of the system.Section 4 provides a discussion of the main results and their implications for tropical forest-savanna ecosystems and their management.

A spatially extended forest-savanna model 2.1 Model overview
The Staver-Levin model (SL model) describes the interactions between savanna trees (T), savanna tree saplings (S), forest trees (F), and grass (G), and was originally motivated by empirical evidence supporting forest-savanna bistability at intermediate rainfall in the tropics [40].In particular, there is evidence that frequent fires can prevent the formation of closed canopy forests and help to maintain savannas in regions that would otherwise be suitable for forest establishment based on their climate [7].
The original SL model [39,42] incorporates spatial extent implicitly by assuming that all vegetation types are spatially well-mixed and hence it tracks the proportion of space occupied by each vegetation type.In this framework, grass represents an "open" patch on which new trees can grow, but grass patches also carry fires that limit the expansion of both savanna and forest trees, albeit via different mechanisms.Forest trees can grow on grass, sapling or savanna occupied patches, at a rate associated with the amount of seeds available on that patch coming from adjacent forest trees.Similarly, savanna tree saplings grow on grass patches at a rate depending on the prevalence of adult savanna trees in their vicinity.Grass represents the primary flammable cover that carries fires which kill forest trees.Forest tree mortality to fire thus depends on the level of grassy cover and this mortality (burning) rate is denoted by ϕ(G).Empirically, fire is very frequently observed in savanna systems below approximately 40% tree cover, but is very rarely observed in savanna systems with more than 40% tree cover [39].Moreover, percolation models of fire spread dynamics support a sharp threshold for fire activation as a function of flammable cover [35].Hence the forest tree burning rate, ϕ(G), depends on grass levels and is typically chosen to be an increasing sigmoidal function with a sharp transition from low-fire to fire-prone regimes.
Fire carried by the grassy layer also impacts the maturation of savanna saplings to adult savanna trees.However, savanna saplings are "top killed", rather than totally destroyed, by fire and their recruitment to adult savanna trees is merely delayed as they are able to resprout later.Thus the sapling-to-savanna recruitment rate, ω(G), depends on grass levels, and is chosen to be a decreasing sigmoidal function with an abrupt transition from high recruitment in the absence of fire to low recruitment in the fire-prone regime.Adult savanna trees are adapted to fire and do not suffer significant excess mortality in a fire-prone environment.Finally, since grass grows much more quickly than the other vegetation types, we assume that whenever savanna trees, savanna saplings or forest trees die, either via natural mortality or fire, they immediately revert back to the grass state.
The interaction rules described above generate the following system of ODEs for the proportions of space covered by each of the four functional types of vegetation: where G denotes the grass cover proportion, S is the savanna saplings proportion, T is the savanna tree cover proportion and F is the forest tree cover proportion.The positive constants µ and ν are the mortality rates of saplings and savanna trees respectively, while α is the forest tree birth rate and β the savanna tree birth rate.The algebraic constraint given by the final equation in (1) ensures that all space is filled by one of the four types of vegetation.Table 1 summarizes the parameters of the SL model, along with their ecological interpretations and default numerical values.
In addition to the four-functional-type model described above, numerous other related models with the similar underlying interaction rules, some including explicit spatial extent, have been studied in the literature [13,12,19,24,35,46,52].
A key feature of the SL model is the threshold response to fire as a function of flammable cover.The nonlinear functions ϕ and ω represent how fire affects tree demography, via forest tree mortality and via the maturation of saplings into savanna trees, and are crucial to the emergence of savanna and forest as alternative stable states in the model.In our theoretical analysis, we assume that ϕ and ω are Heaviside step functions, retaining the key qualitative properties of the fire spread process and simplifying our stability calculations.However, ϕ and ω have smooth, yet sharp, sigmoidal profiles for all numerical investigations (see Table 1).
In recent work, the authors have extended the SL model to a spatially explicit setting by considering interacting particle systems based on the interaction rules outlined above and proving the convergence of these processes to more tractable mean-field limiting processes [31].The distribution of the mean-field limiting processes is governed by the following system of nonlinear integro-differential equations: for each (x, t) ∈ Ω × R + for some Ω ⊂ R 2 with α, β, µ, and ν positive constants, J F , J T , w ∈ L 1 (Ω; R + ) and ϕ, ω ∈ C(R + ; R + ).In our analysis, we will take Ω = R, but when the system is considered on a compact domain for numerical experiments, we employ periodic boundary conditions.Since the system of equations given by (2) describes the evolution of probability densities we retain a normalization condition similar to the space filling constraint for the nonspatial model (1).In particular, G(x, t) + S(x, t) + T (x, t) + F (x, t) = 1 for each (x, t) The kernel function w measures the ability of fire to spread spatially from a point that is already burning.The constants α and β account for the strength of forest-tree and savanna-tree invasion via seed dispersal, with the spatial distribution of these seeds captured by the kernels J F and J T .The inclusion of nonlocal or long-range interactions is considered most appropriate for spatial vegetation models as dispersal of seeds is often long range or even heavy-tailed [29,45].The spatial interaction (fire spread and seed dispersal) are assumed isotropic so all kernels are of convolution type and the model ( 2) is thus posed on a homogeneous spatial domain.In all numerical results, we use zero mean Gaussian kernels with different standard deviations (to reflect the relative length scales of the different spatial process).In particular, we have with w(x) = G(x, σ W ), J T (x) = G(x, σ T ), and J F (x) = G(x, σ F ).Our theoretical results hold for a broad class of kernels obeying some mild assumptions (see below).  Figure 1 shows a two-parameter bifurcation diagram for the nonspatial system (1) with the forest tree birth rate, α, and the savanna tree birth rate, β, chosen as bifurcation parameters.Most parameter regimes have one or more stable equilibrium solutions, ranging from all-grass states to savanna and forest states.Moreover, there is significant multistability for large regions of the parameter space and stable oscillations are also possible in the nonspatial model (yellow region in Figure 1).The interested reader may consult [46] for a more detailed bifurcation analysis of the nonspatial SL model as a function of all system parameters, but Figure 1 provides a representative summary of the possible dynamics.

Pattern formation in the spatially extended SL Model
A key goal of the present work is to explore the dynamics of the spatially extended SL model given by (2), and a natural first step is to understand the stability of the homogeneous equilibria shown in Figure 1 in the presence of spatial interactions.If Ω = R and the kernels J F , J T and w are appropriately normalized in the spatial model, a steady-state solution to (1) is also a solution to (2); it is then evidently of interest to determine if a stable steady-state of (1) remains stable in (2).We study the local stability of an arbitrary homogeneous steady-state by linearizing the system about this fixed point and determining whether or not a perturbation about this solution will be damped or amplified by the system dynamics.In particular, we decompose the perturbation over an appropriate Fourier basis; let ξ ∈ R denote the wave number (Fourier variable) and suppose λ i,ξ for i = 1, ..., 4 denotes the eigenvalues associated to the wave number ξ of the linearized system about some homogeneous steady-state.The steady-state solution will be linearly unstable if for some i = 1, ..., 4: there where ℜ(λ) denotes the real part of λ.Since the homogeneous steady-state was assumed stable in the nonspatial system (1), we have λ i,0 < 0 for i = 1, ..., 4 and to avoid degeneracy (blow-up of solutions) we must also have that lim sup |ξ|→∞ ℜ(λ i,ξ ) < 0. Under these conditions, we expect a branch of stable heterogeneous solutions to emerge as the homogeneous steady-state loses stability, typically referred to as a pattern-forming or symmetry-breaking instability.Since its introduction by Turing [47], this type of mechanism has found applications in a myriad of fields, including neuroscience, ecology, material science, and various branches of biology, chemistry and physics [11,14,26].In particular, there is an extensive literature in ecology on spatial patterning in arid ecosystems and associated modeling efforts to explain these patterns mechanistically [15,16,22,27,50].Here we investigate if the spatial SL model, which is intended as a model for ecosystems at intermediate Points labeled GH denote codimension 2 Bautin (or Generalized Hopf) bifurcation points at which the Hopf bifurcations change criticality.The all-grass state is the only stable state in the dark green shaded region; the all-grass state and a forest dominated state are bistable in the orange shaded region.
rainfall levels, also predicts patterned or heterogeneous vegetation distributions in the absence of any externally imposed spatially heterogeneous structure.
We choose Ω = R for ease of exposition, but the calculations and conclusions are analogous for other standard domains of interest, for example, Ω = R 2 and compact subsets of R 2 with periodic boundary conditions.Furthermore, we restrict attention to a class of kernels satisfying the following assumptions: .) J is an even function.
The normalization condition in (H1.) serves to separate the intensity of spatial processes (such as seed dispersal and fire spread) from their dispersion or variance, and also guarantees that equilibrium solutions of the nonspatial model ( 1) are spatially homogeneous solutions of (2).Condition (H2.) preserves the isotropy of the problem, i.e. we do not consider prevailing wind effects, sloped terrain or other environmental features which might bias the directionality of the fire or seed dispersal processes.For the fire threshold functions, ϕ and ω, we assume that each function is a scaled and shifted version of the Heaviside step function, i.e. (H3.)For some threshold parameter θ 1 ∈ (0, 1) and level parameters ω 0 , ω 1 such that 0 < ω 1 < ω 0 , Similarly, for some θ 2 ∈ (0, 1) and level parameters ϕ 0 , ϕ 1 such that 0 < ϕ 0 < ϕ 1 , In particular, (H3.) implies that ϕ ′ (x) = ω ′ (x) = 0 for almost every x ∈ (0, 1).Carrying out the requisite linear stability analysis under the assumptions above yields the following result; the supporting calculations are deferred to Appendix A.1.
Proposition 1 states that a stable steady-state solution to the nonspatial model will be locally asymptotically stable in the spatially extended Staver-Levin model given by ( 2), regardless of the particular distributions chosen for the seed dispersal kernels of forest and savanna trees, and irrespective of the fire spread kernel (as long as they obey (A1.-A2.)).This result encompasses all symmetric spatial kernels that are probability density functions, regardless of the distribution of their mass, or whether they are thin or heavy-tailed.Hence, on a homogeneous spatial domain, stable solutions of the ODE model ( 1) are likely to dominate the dynamics.Moreoever, extensive numerical simulations of the spatial model ( 2) on a homogeneous domain for a range of parameters and initial conditions did not reveal the emergence of any stable nonhomogeneous solutions.
Many standard pattern formation paradigms, going back even to the seminal work of Turing [47], have the feature of a long-range inhibitory process and a short-range excitation (activation) process (see also [1]).In the present context, seed dispersal, a long-range excitatory process, and fire spread, a relatively shorter-range inhibitory process, are the spatial processes considered and hence we may have intuited the inability of the system to form patterns, as elucidated in Proposition 1.However, Proposition 1 in fact gives the stronger conclusion that pattern formation can be ruled out even without the need to specify the structure of the kernels, and hence the relative spatial scales of the processes (since only (H1.-H2) are assumed and this does not specify which processes are short or long range).The reason that spatial instabilities are not present, even if one were to unrealistically assume that fire was longer range than seed dispersal, is that the spatial impact of fire is essentially a higher order effect due to the structure of our model.Although fire does provide an inhibitory interaction that limits tree growth, its (nonlocal) impact is neglected upon linearization and thus is not sufficiently strong to induce an instability.This is due in particular to the steep nature of the sigmoidal fire onset functions ϕ and ω, which both have the property that their derivatives are almost everywhere zero (see Appendix A.1 for further details).We can therefore conclude that our model must take into account additional spatial processes that limit (inhibit) tree growth to make pattern formation possible via a Turing-type mechanism.Remark 1.In numerical experiments with fire threshold functions ϕ and ω are smooth steep sigmoids, we did not observe any pattern forming regions apart from those predicted by the linear stability analysis with Heaviside fire threshold functions.

Resource limitation effects
In many spatially extended ecological models resource limitation is an essential ingredient for the existence of pattern-forming instabilities [15,16,22,27,50].Pattern formation in nutrient deficient savanna ecosystems has been studied from both empirical and theoretical perspectives [5,23] with heterogeneous patterns reported as being particularly common at the interface of savannas and tropical forests [33].Empirical studies support the hypothesis that forests and savannas are alternative stable states at intermediate mean annual rainfall (MAR) in the tropics [40], and this potential for bistability is well captured in the standard SL model interactions (see (1) and Figure 1).However, closed canopy forest ecosystems (i.e.steady states with high values of F ) should be more favored, even within the bistable range of MAR, as rainfall increases (or the terrain becomes more hospitable with regard to the availability of other resources).This effect is particularly important if we are interested in modeling the forest-savanna transition as it will have strong implications for estimates on the resilience of forest or savanna ecosystems, as well as their abilities to potentially invade one another (cf.[52]); introducing resource limitation to the model serves to account for this mediation of the competition between the vegetation types.With regard to the resources necessary to support savanna or forest trees, we have in mind water and soil nutrients, such as nitrogen, which forest trees typically require in more abundance than savanna trees, saplings or grass.
Consider the grass-forest SL model (ignoring saplings and savannas trees for now) without spatial extent and add a resource limitation effect on forest trees to obtain the following nonspatial model: where r ∈ [0, 1] denotes the resource constraint parameter and G + F = 1.In Figure 2, we have a numerical bifurcation analysis of system (5); the system exhibits either a stable all-grass solution (green region), bistability between a grassland and a forest dominated solution (blue region), or is monostable (white region) with a solution which contains both grass and forest in proportions that depend on both α and r.The regions in which each of these possibilities occurs can be viewed in α-r space in the two-parameter bifurcation diagram shown in Figure 2 C. We now evaluate the stability of each of these steady-state solutions in a spatially extended model.Analogous to the nonlocal spatial extensions shown before, the spatially extended version of ( 5) is given by the following system of integro-differential equations: where the kernel R is a probability density function summarizing the degree of competition that trees exert on one another at various distances.We assume that R obeys (A1.-A2.) and once more fix Ω = R for definiteness, but the following calculations and conclusions can of course be extended to more general domains.If R is a Dirac delta distribution, then the resource constraint is analogous to that of the classic Fisher-Kolmogorov equation.Since F (x, t) = 1 − G(x, t) for each (x, t) ∈ R × R + , the grass-forest system reduces to the single integro-differential equation: The grass-forest SL model with resource limitation on forest tree growth given by ( 8) can exhibit pattern-forming instabilities and, via numerical bifurcation analysis, we can further study the structure of these bifurcations and their implications for the ecosystems under consideration.When the resource limitation on forest trees is added to the 4-species model, the only possible spatially induced bifurcations of equilibria are those coming from the forest-grass subsystem (see Section 3.2 for details); hence the forthcoming detailed analysis of the forest-grass subsystem allows us to understand the source of instabilities in the simplest setting possible.Our analysis shows that threshold fire feedback or disturbance is compatible with pattern formation in resource limited biomes and that there are multiple different pattern formation paradigms in different parameter regimes of α-r-space and for different relative spatial scales of seed dispersal and resource competition.In particular, we find several different bistable regimes resulting from backward (subcritical) bifurcations of the homogeneous steady states.

Stability of steady-state solutions with resource limitation
We introduce one additional assumption in order to classify the parameter space into regions that can or cannot exhibit pattern formation in the system (8).The key condition for our calculations will be whether or not the Fourier transforms of the kernel functions are nonnegative or not.A simple sufficient condition for a dispersal kernel J to have a nonnegative Fourier transform is that: (H4.) J ∈ L 1 (R; R + ) attains its maximum value at zero.Hypothesis (H4.) captures virtually all of the kernels typically used in vegetation modeling, including the Gaussian, exponential, and Cauchy kernels.A necessary and sufficient condition for a function to have a nonnegative Fourier transform is that the kernel be positive definite [6]: Definition 1 (Positive definite function).A function f : R → C is positive definite if for any real numbers x 1 , . . ., x n the n × n matrix A whose entries are given by A i,j = f (x i − x j ) is positive semi-definite.
In most practical cases (H4.) is sufficient and much easier to verify but the forthcoming analysis is unchanged and the results remain valid for the more general class of positive definite kernels.The following result rules out spatially induced interactions destabilizing the all-grass state under the conditions above; the supporting calculations are contained in Appendix A.2.
Proposition 2. If the all-grass state is a (locally asymptotically) stable equilibrium of the nonspatial resource limited SL model given by (5), then it is also a stable spatially homogeneous equilibrium of the spatially extended resource limited SL model given by (8) with Ω = R with positive definite kernels obeying (H1.), (H2.), and ω and ϕ obeying (H3.).
Having ruled out bifurcations of the all-grass solution, we now evaluate the stability of mixed equilibria in which both the grass and forest types are present.Define the quantities The dispersion relation for any homogeneous equilibrium Ḡ ∈ (0, 1) is then given by where we recall that ξ denotes the wave number (see Appendix A.2 for the derivation of this expression).For Ḡ to lose stability and potentially generate new stable heterogeneous solutions the following conditions must hold: (i.) λ 0 < 0, i.e.Ḡ is stable in the ODE ( 5), (ii.) there exists ξ ̸ = 0 such that λ ξ > 0, (iii.)λ ξ < 0 for all |ξ| sufficiently large.
Since Ḡ ∈ (0, 1) obeys condition (10) is equivalent to the inequality which is true for every equilibrium Ḡ ∈ (0, 1) since the steady state equation ( 11) cannot be satisfied if 1 − Ḡ ≥ r due to the strict positivity of ϕ.Therefore the non-degeneracy condition (10) is satisfied for all equilibria Ḡ ∈ (0, 1) and we are guaranteed that condition (iii.)holds.
Due to (12), we have Hence (13) gives us a necessary condition for (ii.) to hold, namely that for a given homogeneous solution Ḡ.Note that since N (α, r) does not depend on the kernel functions, our necessary condition for instabilities does not depend on the nature of the spatial interactions, but only on their overall strengths through the parameters α and r. Figure 4 A shows the values of of α and r for which N (α, r) > 0 for at least one equilibrium that is stable in the ODE model.We only consider values of α > 0.9 because, for α < 0.9, the all-grass state is the only stable solution of the nonspatial model.Since N (α, r) is a multi-valued function in bistable regions, we take the maximum of N (α, r) across the two stable equilibria in the bistable regions.Instability of the homogeneous solutions is ruled out (regardless of the nature of the spatial interactions) in the white regions which span most of the parameter space.However, there are significant regions highlighted in red where pattern-forming instabilities are possible.
In Figure 3 B, we take Gaussian kernels and plot the dispersion relation for the equilibria which are stable for the ODE (5) for two sets of parameters, one in the monostable and one in the bistable regime.For the bistable parameter set (blue curves), the dispersion relation of the grass-dominated equilibrium with Ḡ ≈ 0.62 is lower blue curve and hence this homogeneous equilibrium does not lose stability in the presence of spatial interactions.However, the dispersion relation of the corresponding forest-dominated equilibrium (with Ḡ ≈ 0.24) becomes positive for a range of wave numbers and hence does lose stability.The dispersion relation in red corresponds to a mixed forest-grass equilibrium in a monostable regime which also loses stability upon the introduction of appropriate spatial interactions.
Figure 3 C shows the pattern forming instability region across a large region of α-r parameter space; this amounts to calculating the dispersion relations at each point and recording the maximum achieved across both equilibria and all wave numbers.The red shaded regions are the so-called Turing space where a stable homogeneous solution has lost stability in the presence of spatial interactions due to a real eigenvalue becoming positive.The plots shown here are with Gaussian kernels for seed dispersal, fire spread and resource limitation.The standard deviation of the Gaussian fire spread kernel (σ W ) is smaller than that of seed dispersal (σ F ), which is in turn smaller than that of the resource competition kernel (σ R ).As we show presently, pattern-forming instabilities persist as the relationships between these parameters are varied over relatively wide ranges and the instabilities do not necessarily depend on having σ R < σ F .
In Figure 4 we confirm the pattern formation predicted in Figure 3 by solving the spatial model with periodic boundary conditions.Figure 4 B and C show stable heterogeneous steady state solutions for the system (7) for parameter values predicted to exhibit pattern formation in panel B of Figure 3.In both cases, a homogeneous solution loses stability and gives birth to a new stable heterogeneous solution; the solution shown in panel B is in the bistable regime and the second homogeneous equilibrium retains stability, while the solution in panel C is not in the bistable regime and appears to be the only stable solution in this parameter regime (up to translation).Panel D shows some of the complex patterns that arise on a 2 dimensional spatial domain when the homogeneous solution Ḡ ≈ 0.24 undergoes a symmetry breaking bifurcation; the rich diversity of possible patterns is due in part to the periodic boundary conditions used in these simulations.Figure 5 shows the Turing regions in panels A and D for two different values of σ R , the standard deviation of the resource competition kernel, with the kernel much more localized in D than in A, leading to a much smaller Turing region.The Turing regions are considerably larger when σ R is larger (panel A vs panel D), but the situation in panel A is only realistic when we consider resource competition between lower lying vegetation, since this case requires resource competition to be more nonlocal than seed dispersal.The horizontal dashed magenta line in Figure 5 A indicates the one-parameter slice of the diagram corresponding to the bifurcation diagram drawn in panel B and the magenta star indicates the position in α-r space of the solution shown in panel C. Panel B of Figure 5 is a one-parameter bifurcation diagram of the spatially extended model 8 as the seed dispersal rate, α, is varied; the L 1 norm of the grass component of the solution shown on the yaxis, stable homogeneous steady states are in red, unstable homogeneous steady states in black, stable heterogeneous steady states and unstable heterogeneous steady states are in green and blue respectively.The homogeneous all-grass solution is not shown but is always unstable for the ranges of α depicted here.
In Figure 5 B, we observe changes in stability in the homogeneous steady states curves around α ≈ 4.35 and α ≈ 4.75 (at the intersections of the red and black curves).These two bifurcation points are distinct from the saddle node points of the nonspatial model, indicating that these are instabilities induced by the spatial interactions, matching exactly the predictions of the linear stability analysis in Figure 5 A. Numerical continuation reveals that the bifurcations are in fact subcritical and that there are multiple stable heterogeneous solutions.In the bifurcation analysis presented here we show the two stable heterogeneous solution curves which span the widest ranges of α while remaining stable and which also appear to have the largest basins of attraction (based on solving the initial value problem via timestepping).Due to the periodic boundary conditions, there appear to be multiple other unstable solution branches (not pictured), some of which become stable for relatively narrow ranges of the bifurcation parameter.Further analytical work and different methods would be required to understand this behavior systematically, but we emphasize instead the more ecologically relevant aspect of the analysis, i.e. the subcritical nature of the bifurcations and the width of the multistable region.There is a significant region of multistability between the stable heterogeneous solutions (green curves) shown and the stable homogeneous solutions (red curves).instability in the homogeneous solutions, and persists until α ≈ 8.44, when it disappears in another saddle-node bifurcation.The second heterogeneous solution curve is a three-bump solution which similarly appears and disappears in what appear to be saddle-node type bifurcations.Figure 5 D, E and F are analogous to the top row of the same figure in content but focus on a scenario likely to be more realistic for the types of vegetation present in forest-savanna ecosystems, i.e. short scale fire lengthscale (σ W = 0.025), intermediate scale resource competition (σ R = 0.05) and longer scale seed dispersal (σ F = 0.1).As Figure 5 D shows, the Turing regions are relatively small in this case but the subcritical nature of the bifurcations and the significant early onset of patterning before the bifurcations remains present, as we can see in Figure 5 E. Indeed, there is still a small region of bistability beyond the point at which the homogeneous branch regains stability around α ≈ 6.08 in panel E, although it is harder to observe here than in panel B. In this case, the heterogeneous branch appears to both gain and lose stability in saddle-node bifurcations, with the unstable branches connecting to the homogeneous branches at the Turing bifurcation points predicted by the linear stability analysis.The solution shown in Figure 5 F now has a much higher wavelength than before due to the more localized competition (smaller σ R ).
Our model has three spatial scales: the fire spread scale, the seed dispersal scale and the resource competition scale.We next briefly explore the relationship between the emergence of patterns and the relative values of these scales.Since the dispersion relation that characterizes the onset of pattern formation is not analytically tractable, we resort to numerical analysis to understand the persistence of patterns as we vary the relationships between the spatial scales.As we have throughout, we assume centered Gaussian kernels so that the length scales are summarized by the standard deviation parameters σ W , σ F and σ R .We take the fire dispersal parameter σ W = 0.01 and we fix the seed dispersal parameter σ F = 0.1, since forest tree seeds should typically disperse Figure 6 shows the maximal value of the dispersion relation across all homogeneous steady states of the nonlocal model ( 8) in r-σ R space for fixed values of σ F = 0.1 and σ W = 0.01 and various values of α.As usual, pattern-forming parameter regions are highlighted in red, with darker red indicating faster growing instabilities.These plots reveal that the pattern-forming instabilities persist for a much larger range of values of σ R than we may have anticipated and are still present even when the σ R is significantly less than σ F (σ F = 0.1 in all panels).In fact, pattern forming instabilities persist for σ R as low as ≈ 0.02 across a wide range of values of the resource constraint and forest birth rate parameters (r and α respectively).In panels B and C, the nonspatial system is in monostable parameter regimes, but the system is bistable for part of the range of r values in panel A (between the saddle nodes bifurcations indicated by the vertical dashed green lines).
From an ecological perspective, the large range of pattern forming values of σ R is crucial, because it means that patterns may be found even where the distance-scale of resource competition is smaller than that of dispersal; this may occur widely in nature, since resource competition between trees is often considered to be relatively localized compared to dispersal [30].Thus, although the Turing space grows with σ R , it is not crucial to our conclusions σ R be as large as σ F , or that resource limitation be of a similar scale to seed dispersal.Other classical models typically require resource competition to be more nonlocal than growth (or seed dispersal) processes have been to pattern in arid or ecosystems with lower lying plant species that disperse very locally [22].

The four-type SL model with resource limitation
The four-type nonspatial SL model with resource limitation is given by Figure 7 shows a two-parameter bifurcation diagram for the system (14) as a function of the forest tree birth rate (α) and the savanna tree birth rate (β) for a fixed value of the resource constraint (r = 0.84).The presence of the resource constraint pushes the stable forest dominated regions of parameter space to higher values of α when compared to Figure 1, but many of the qualitative features of the model remain similar to the non-resource limited case.We neglect higher values of α in Figure 7 since the dynamics become confined to the forest-grass subsystem which was studied extensively in Section 3.1.Stable oscillations are still present in the resource-limited model (yellow shaded region), but in the analysis that follows we continue to focus on homogeneous equilibrium solutions and their stability within the spatially extended version of the model.We did not observe a loss of stability of any homogeneous periodic solutions due to the presence of spatial interactions for the kernels and parameters studied here, but this would be an interesting phenomenon to study in more detail and there are techniques available suited to this problem [34].
Remark 2. Figure 7 shows a number of bifurcations around α ≈ 1.65 and β ≈ 1.65.While bifurcation curves and codimension two points may appear to overlap, no additional degeneracy was observed in the bifurcations.Stable oscillations disappear in this region through a complex sequence of homoclinic and fold of limit cycles bifurcations.While this phenomenon is not of practical interest for this paper, especially since it is confined in a small region of the parameter space, the geometry of the flow in this region is complex and potentially of independent interest.The dynamics of the spatial version of the resource limited model are governed by the following system of equations: for each (x, t) ∈ Ω × R + .
To investigate instabilities of homogeneous equilibria of ( 15), consider a stable equilibrium of the nonspatial model given by ( 14), Ḡ, S, T , F , which is also a homogenseous equilibrium of (15).Suppose that hypotheses (H1.-H3.)hold and assume that the spatial kernels are positive definite functions.Upon carrying out the standard linearization procedure (see Appendix A.1 for details), we are left to study the eigenvalue problem λ ξ ⃗ u = A ξ ⃗ u for ⃗ u = (ĝ, ŝ, τ , f ) where A ξ is given by and ξ denotes the wavenumber of the perturbation.Thus, by expanding the determinant of A ξ along the fourth row, we can read off that the first set of eigenvalues are given by The expression (17) is exactly what we previously obtained when studying the forest-grass subsystem under resource limitation (cf.equation ( 9)), and only relates to interactions and parameters involving forest and grass alone.Hence our previous analysis of ( 17) applies once more.In particular, the homogeneous all-grass solution remains stable when it is stable in the nonspatial model and stable patterns can be found in parameter regimes where the forest-grass subsystem is stable with respect to perturbations (S = T = 0 for all x ∈ Ω).
The other components of the linearized spectrum will be given by the eigenvalues of the reduced coefficient matrix: Case (i.): S = T = 0 If we consider stability of homogeneous equilibria within the forest-grass subsystem, then S = T = 0 and the reduced coefficient matrix above becomes Thus we can read off that the second set of eigenvalues are given by λ ξ,2 = −α F 1 − F /r < 0 for all ξ ∈ R but it is constant in ξ and thus not the source of a potentially instability in the homogeneous solution.In fact, λ ξ,2 is negative by hypothesis since we assume the equilibrium to be stable for the mean-field model.The final two sets of eigenvalues are those of the coefficient matrix If the trace of the matrix above is negative and the determinant positive for all ξ ∈ R, then there can be no bifurcation arising from these eigenvalues.The trace is negative by inspection (it does not depend on ξ) and the determinant is given by where ) .Therefore the forest-grass subsystem solutions can only lose stability if there exists an ξ > 0 such that λ ξ,F > 0 where λ ξ,F are the eigenvalues associated with the forest/grass types given by ( 17) and there will be no other unstable modes arising from the rest of the system.As we have seen before, this typically leads to stable heterogeneous steady states via (subcritical) Turing bifurcations.Case (ii.): S > 0 or T > 0 It remains to check whether there could be an instability of a homogeneous steady state in which savanna and saplings are present arising from the eigenvalues associated with the reduced coefficient matrix (18).By using the third order Routh-Hurwitz criteria on (18), it can be shown that no eigenvalues of this matrix can have positive real-part when F ≤ r, even when we are not necessarily in the forest-grass subsystem, i.e. S > 0 or T > 0. However, we cannot have a homogeneous equilibrium with F > r since it will not be able to satisfying the steady state equation for forest trees (see equation (14d)).Therefore pattern-forming bifurcations only arise in the four-type resource limited model given by (15) when an eigenvalue λ ξ,F from equation (17) gains positive real-part.Hence we have reduced the analysis of the highly complex dispersion relation associated with the matrix (16) to the relatively simpler dispersion relation given by (17), which we previously studied in detail in Subsection 3.1.

Discussion and conclusions
We have outlined a spatially explicit modeling framework for savanna-forest ecosystems at intermediate rainfall and shown that this system admits pattern-forming instabilities leading to stable heterogeneous landscape structures.Moreover, this is only possible when we augment the Staver-Levin model with an additional resource constraint on forest tree growth.The lack of patterns without this additional feature is not necessarily unexpected to those familiar with typical patternforming motifs in nonlinear PDEs, but does serve to emphasize the complexity inherent in building a realistic spatial mesic savanna model.As our linear stability analysis shows, the spatial inhibition of the fire process is effectively higher order in our model and thus not sufficiently strong to induce instabilities in the homogeneous solutions, regardless of the assumptions on the relative spatial scales of fire spread and seed dispersal.In effect, our addition of a resource limitation term that lowers the growth rate of forest trees adds a strong inhibitory interaction to the model and thus adds a key ingredient of many standard pattern formation paradigms.
From an ecological perspective, the resource constraint in our model is currently incorporated in a purely phenomenological way, but there are several candidate spatial processes that may limit forest tree growth in practice.Two possible distinct spatial processes generating this effect are nutrient and water competition, with nutrient competition referring to soil nutrients such as nitrogen.Trees in nutrient poor savannas have been empirically observed to exhibit regular spatial patterns [23] and patterns are a key feature of hydrologically coupled vegetation models for semi-arid ecosystems [28].
Another key spatial process in savanna ecosystems, which remains relatively less explored from a modeling perspective, is herbivory by large mammals [38,41].In the case of savanna-forest mosaics, browsing by herbivores like elephants [9] or nyala [21] along the forest edge could reduce tree growth, thereby preventing forest patches from expanding into savannas.For herbivory to have the type of effect that would stabilize a savanna-forest mosaic, browsing would have to intensify with increasing forest patch size; this is counterintuitive since, given the same number of herbivores concentrated on a smaller patch would tend to intensify their effect [2,49].Consider, however, a savanna-forest mosaic in a system like Lopé National Park in Gabon.Hypothetically, the existence of the savanna-forest mosaic with an extensive network of edges likely allows the persistence of a larger elephant population than would occur otherwise [9], which in turn could potentially slow the recruitment of forest trees [44] and thus help prevent forest from taking over savanna patches.In effect, coupling herbivore population dynamics to landscape vegetation composition could produce emergent dynamics that are currently neglected in most theoretical models for herbivore impacts in savannas.This potential novel mechanism for stabilizing vegetation patterns would benefit from further theoretical (as well as empirical) investigation.
Numerical simulations and continuation analysis revealed the subcritical nature of the patternforming bifurcations in our model, giving rise to the onset of stable heterogeneous solutions well before the spatially homogeneous steady states lose stability.These numerical results emphasize the necessity of going beyond near equilibrium analysis to understand the possible transitions that might occur as systems are subjected to exogenous perturbations (see Figure 5).We showed two heterogeneous solution curves in our continuation results, but there are many more unstable solutions curves not shown, some of which appear to snake back and forth as the bifurcation parameter is varied and may gain stability over relatively narrow parameter ranges.We did not show these curves for visual clarity, but it would be interesting to investigate the nature of this behaviour and to understand to what extent it persists when the symmetry induced by the periodic boundary conditions is broken.
The main conclusions of the bifurcation analysis are robust to variation in the relative spatial scales in the problem, particularly the length scales of the seed dispersal and resource limitation processes.Although pattern formation occurs over wider ranges of parameter space when resource competition is more nonlocal than seed dispersal, there are still significant regions of parameter space that omit patterns when seed dispersal is considered to be more long range than competition for resources.This is an important consideration for real-world applications as, in contrast with low lying vegetation patterns such as the famous tiger bush [22], we expect wind dispersal of forest tree seeds to be the longest range spatial process in our system.From an ecological standpoint, the robust persistence of patterned vegetation in mesic savanna/forest systems (referred to in the empirical literature as a "savanna-forest mosaic" [30]) is also significant, and suggests that the combination of fire, dispersal limitation, and resource constraints may contribute to stabilizing mosaics in intermediate rainfall landscapes.
A central goal for the ecological and mathematical modeling community interested in savannaforest ecosystems is a unified mathematical model able to explain observed savanna-forest configurations over a broad range of mean annual rainfall (and subject to local topography).Although such a unified forest-savanna model is crucial for prediction and conservation efforts, no such model currently exists in the literature [32].Such a unified model would be required to reproduce observed spatial patterns, as shown here and in other papers focusing on drier savannas [4,17,48], and provide realistic ecosystems descriptions by reflecting the diversity of functional types within different tropical ecosystems.We believe the model proposed here is somewhat unique in satisfying both of these criteria while retaining sufficient tractability to allow detailed analysis of solutions.
Finally, while this study is limited to a spatially homogeneous domain, and hence most appropriate for modeling at relatively small spatial scales, we note that there is growing acknowledgment in the modeling community of the key role that heterogeneity (or underling spatial pattern of the domain) plays in determining the final spatial patterns that we observe empirically [3,32].This consideration is particularly relevant in the present context as savannas are subject to significant continental scale rainfall gradients in Africa and Amazonia [8,51], and thus we view the present work as a prerequisite step for a more realistic mesic savanna model that would produce patterns via an interplay between emergent effects and externally imposed spatial pattern owing to domain heterogeneity.Significant challenges remain as we progress towards a unified savanna-forest modeling framework that could be empirically validated, and ultimately this is likely to require integrating hydrology and heterogeneity, along with the functional vegetation types and spatial interactions presented here.equilibrium solution ( Ḡ, S, T , F ) and truncate at first order in the perturbations to obtain the linearisation of the system about the chosen spatially homogeneous steady-state.After eliminating the exponential prefactors, we see that the eigenvalues λ ξ associated to the wave number ξ obey: Thus for each ξ ∈ R, we have the eigenvalue problem where the coefficient matrix A is given by At this point, no nonlocal effects of fire remain in the model.Expanding along the fourth row of the matrix above, we immediately identify that the component of the spectrum relating to forest trees can be considered in isolation and is given by This component of the spectrum is stable in the mean-field model if and only if Condition (20) corresponds to asking that this part of the spectrum is stable for the zero Fourier mode, i.e. λ F (0) < 0. In order for a Turing bifurcation to arise here we assume that (20)  Hence this component of the spectrum is stable in the spatially extended model whenever it is stable for the mean-field model and therefore cannot destabilize the spatially homogeneous steady-state.
Next consider the eigenvalues relating to the reduced coefficient matrix, i.e.

A.2 Stability of homogeneous equilibria with resource limitation
Once more we linearize about a fixed point of the mean-field version of the model, i.e. we consider a stable equilibrium solution Ḡ of ( 5) and then take a perturbation of this equilibrium solution of the form G(x, t) = Ḡ + ϵ g(x, t) = Ḡ + ϵ R ĝ(ξ, t)e iξx dξ, where ĝ denotes the Fourier transform of the perturbation term g and ĝ(ξ, t) = e λ ξ t ĝ(ξ, 0) ≡ e λ ξ t ĝ(ξ) denotes the amplitude associated with the wave number ξ ∈ R. Inserting the perturbed solution from ( 23) into ( 8) and discarding second order (and higher) terms in ϵ shows that to leading order λ ξ obeys: for each ξ ∈ R. Our kernels are all even functions by assumption, so we obtain the simplified dispersion relation for each ξ ∈ R. If we consider the all grass state Ḡ = 1 (which is stable for all values of α before the transcritical bifurcation denoted by the vertical blue line in Figure 2 C), then equation ( 24) reduces to λ ξ = −ϕ (1) + α ĴF (|ξ|) for each ξ ∈ R. Since the Fourier transform of a nonnegative function achieves its maximum at zero, sup ξ λ ξ = −ϕ (1) + α ĴF (0) = −ϕ(1) + α and hence the stability criterion for the all grass state is the same as it was for the mean-field model (1).Thus if Ḡ = 1 was stable for the nonspatial model ( 5), it will remain stable for the spatially extended model given by (7).This means that in r-α space, the all grass steady-state is stable in the spatial model in the green shaded region of parameter space from Figure 2 D if we start sufficiently close to Ḡ = 1, regardless of the nature of the spatial interactions (i.e. for any kernels obeying assumptions (H1.) and (H2.)).

Figure 1 :
Figure 1: Two-parameter bifurcation diagram in α (forest tree birth rate) and β (savanna tree birth rate) for the nonspatial Staver-Levin model (1).Transcritical bifurcation curves in blue, saddle node curves in magenta, supercritical Hopf curves in purple and subcritical Hopf curves in dark green.Points labeled GH denote codimension 2 Bautin (or Generalized Hopf) bifurcation points at which the Hopf bifurcations change criticality.The all-grass state is the only stable state in the dark green shaded region; the all-grass state and a forest dominated state are bistable in the orange shaded region.

Figure 2 :
Figure 2: A/B: One-parameter bifurcation diagrams for the resource limited model for various values of the resource constraint r; stable fixed points in red and unstable fixed points in black.C: Two-parameter bifurcation diagram in r and α (forest tree birth rate); curve of transcritical bifurcations in blue, curves of saddle node bifurcations in magenta.

Figure 3 :
Figure3: A: Red shaded regions satisfy the necessary condition for an instability of a spatially homogeneous solution, i.e. max (N (α, r)) > 0, where the max is taken over all stable equilibria of the nonspatial system.B: Plots of the dispersion relations for the spatial model for two parameter sets varying α and r (the system is bistable in one case and monostable in the other).C: Maximum of the principal eigenvalue of the linearized system with red regions denoting pattern-forming parameter regimes.Parameters: Gaussian kernels with σ W = 0.05, σ F = 0.1 and σ R = 0.4.

Figure 5 C-Figure 4 :
Figure5shows the Turing regions in panels A and D for two different values of σ R , the standard deviation of the resource competition kernel, with the kernel much more localized in D than in A, leading to a much smaller Turing region.The Turing regions are considerably larger when σ R is larger (panel A vs panel D), but the situation in panel A is only realistic when we consider resource competition between lower lying vegetation, since this case requires resource competition to be more nonlocal than seed dispersal.The horizontal dashed magenta line in Figure5A indicates the one-parameter slice of the diagram corresponding to the bifurcation diagram drawn in panel B and the magenta star indicates the position in α-r space of the solution shown in panel C. Panel B of Figure5is a one-parameter bifurcation diagram of the spatially extended model 8 as the seed dispersal rate, α, is varied; the L 1 norm of the grass component of the solution shown on the yaxis, stable homogeneous steady states are in red, unstable homogeneous steady states in black, stable heterogeneous steady states and unstable heterogeneous steady states are in green and blue respectively.The homogeneous all-grass solution is not shown but is always unstable for the ranges of α depicted here.In Figure5B, we observe changes in stability in the homogeneous steady states curves around α ≈ 4.35 and α ≈ 4.75 (at the intersections of the red and black curves).These two bifurcation points are distinct from the saddle node points of the nonspatial model, indicating that these are instabilities induced by the spatial interactions, matching exactly the predictions of the linear stability analysis in Figure5A.Numerical continuation reveals that the bifurcations are in fact subcritical and that there are multiple stable heterogeneous solutions.In the bifurcation analysis presented here we show the two stable heterogeneous solution curves which span the widest ranges of α while remaining stable and which also appear to have the largest basins of attraction (based on solving the initial value problem via timestepping).Due to the periodic boundary conditions, there appear to be multiple other unstable solution branches (not pictured), some of which become stable for relatively narrow ranges of the bifurcation parameter.Further analytical work and different methods would be required to understand this behavior systematically, but we emphasize instead the more ecologically relevant aspect of the analysis, i.e. the subcritical nature of the bifurcations and the width of the multistable region.There is a significant region of multistability between the stable heterogeneous solutions (green curves) shown and the stable homogeneous solutions (red curves).Figure5C shows two bistable heterogeneous solution for α = 4.5 and r = 0.84 with the color of the star in panel B indicating which solution belongs to each stable branch.A two-bump solution curve appears in a saddle-node bifurcation around α ≈ 2.92, significantly before the linear stability analysis predicts an

Figure 5 :
Figure 5: Top row parameter regime: σ F = 0.1, σ W = 0.025, σ R = 0.15.Bottom row parameter regime: σ F = 0.1, σ W = 0.025, σ R = 0.05.A/D: Heatmap of the maximum dispersion relation obtained by linearizing about each (nonspatially) stable homogeneous equilibrium.The dashed horizontal magenta line corresponds to the one-parameter bifurcation diagrams in B/E and the magenta star indicates the position of the solutions shown in C/F respectively.B/E: Oneparameter bifurcation diagrams for the spatial model varying the forest tree birth rate (α); the positions of the solutions shown in panel C are indicated with a star (in panel B the color of the star matches of the color of the solution curve in panel C).C/F: Steady-state profiles of stable heterogeneous solutions with two bistable solutions shown in panel C.

Figure 6 :
Figure 6: A/B/C: Heatmaps showing the regions where the dispersion relation is positive for at least one homogeneous equilibrium in the absence of spatial interactions, i.e. the pattern forming region, as a function of σ R (the standard deviation of the resource competition kernel) for different values of α (the forest tree birth rate).The green dashed lines in panel A indicate the position of saddle node bifurcation points (and hence the bistable region) in the non-spatial model.Other parameter values: σ F = 0.1 and σ W = 0.01.

Figure 7 :
Figure 7: Two-parameter bifurcation diagram in α (forest tree birth rate) and β (savanna tree birth rate) for the nonspatial resource-limited Staver-Levin model (14) with r = 0.84.Transcritical bifurcation curves in blue, saddle node curves in magenta, supercritical Hopf curves in purple and subcritical Hopf curves in dark green.Points labeled CP denote codimension 2 Cusp bifurcations and points labeled BT denote codimension 2 Bogdanov-Takens bifurcations.Points labeled GH denote codimension 2 Bautin (or Generalized Hopf) bifurcation points which are points at which the Hopf bifurcations change criticality.The all-grass state is the only stable state in the dark green shaded region (bottom left corner).

Table 1 :
Summary of parameters and their interpretations for the SL model