Emergent structure and dynamics of tropical forest-grassland landscapes

Significance Tropical forests may face the risk of abrupt dieback due to amplifying feedback between forest loss and fire spread. Considering the patch-scale rules of forest and fire spread, we find that forest expands into grassland at a rate proportional to its perimeter, while it recedes at a rate proportional to its perimeter and the area of adjacent grassland. Looking at the landscape-scale balance of changes in forest area, we find that these two quantities respectively appear in the gain and loss parts of the equation. Such a relation between spatial structure and expected change of forest area can help identify which parts of the landscape are best targeted for conservation or restoration to avert forest dieback.


INTRODUCTION
Satellite [1,2] and ground observations [3,4] show that tropical forest (high tree cover) and tropical savanna (low tree cover) can exist under the same environmental conditions, making the distribution of tree cover bimodal.On the one hand, fire exclusion experiments have shown that fire can maintain low tree cover [5].On the other hand, fire occurs almost exclusively below a tree cover threshold of about 40% [1,[6][7][8][9], which is consistent with fire being a contagion process on grass patches [10,11], while tree patches block fire.Such a highly nonlinear response of fire to grass together with an empirically consistent response of vegetation to fire was shown to be sufficient for inducing bistability in simple models [12].Taken together, the bimodality, the mutual interaction between fire and vegetation, and the availability of a plausible underlying mechanism suggest that tropical forest and savanna are alternative stable states, maintained by a feedback between vegetation and fire [1], and between which transitions would neither be gradual nor easily reversed [13,14].
Bistability of forest and savanna has been studied with a variety of modelling approaches, which can be classified as microscopic versus mean-field models.The underlying processes concern the spatiotemporal population dynamics of discrete vegetation patches, which can spread or block fire.These can be most realistically modelled by microscopic models, such as interacting particle systems [15] or cellular automata [16], which consider the stochastic dynamics of such discrete constituents interacting in a spatial domain according to simple rules.However, as microscopic models are hard to analyse, one usually looks for a coarse-grained approximation that permits analysis.Mean-field models provide such an approximation, typically in the form of a small number of differential equations that describe the average properties of the considered populations through time, such as cover fractions of each species.If the averages are taken over the whole landscape, the resulting mean-field model is non-spatial and describes macroscopic dynamics via ordinary differential equations (ODEs) [9,12,17].If averages are taken over a neighbourhood, the mean-field model is spatial and describes the dynamics on a mesoscopic scale, via partial differential [18,19], spatial difference (20; spatial mean field in 21), or partial integro-differential equations [mean field in 15].Mean-field models owe their simple closed form to an assumption of statistical independence between species' occurrences in space [e.g.22,23], which permits writing the interaction between any two species as the product of their occurrences.However, assuming statistical independence in space implies neglect of spatial structure.
Despite their disregard for spatial structure and resulting biases [e.g.22,24], mean-field models have been indispensable tools for gaining theoretical insight in alternative stable tree cover states in the tropics.The Staver-Levin model of tropical tree cover bistability [12] is a non-spatial mean-field model in which the variables represent grass and tree cover fractions in the landscape, with interaction between species captured as the product of their cover fractions.Fire spread is not included explicitly.Instead, the effects of fire on vegetation are implicitly accounted for by making the relevant conversion rates a threshold function of grass cover, where the threshold corresponds to the point where large contiguous grass patches emerge, also known as the percolation threshold [25,26].The Staver-Levin model has provided a first proof of principle for alternative stable tree cover states in the tropics, and showed additional complex behaviours, such as cycles and stochastic resonance [12,17].Spatial mean-field models of the Staver-Levin model fur-ther showed emergent phenomena due to spatial interaction on mesoscopic scales, such as travelling and pinning fronts between states (8, 18, 19, 21; spatial mean field in 15), front curvature effects [19,21], pattern formation [27], and coexistence states [28].Even though they are spatial, they are still mean-field models, as they do not consider the fundamental spreading processes of forest and fire on patches, but use equations, with implicit assumptions on the spatial structure of the patches at finer scales than those modelled.
The effect of this fine-grained spatial structure can only be studied via microscopic models.Schertzer et al. (2014) [10] proposed a cellular automaton implementation of the Staver-Levin model in which the effect of fire is still captured implicitly, as a threshold function of flammable vegetation.The form of this vegetationfire relation was obtained from separate simulations of fire spread as a standard percolation process.The cellular automaton and its mean-field approximation were shown to exhibit bistability.Thereby, Schertzer et al. [10] provided the first mechanistic explanation of the role of fire as a percolation process in bistability of tropical tree cover.It also justifies the qualitative form of the fire-vegetation dependence assumed in mean-field models.The more recent interacting particle system by Patterson et al. (2021) [15] follows a similar approach, by implementing fire as a threshold function of neighbourhood grass cover, where the threshold is assumed to match with that of site percolation [25,26].However, standard percolation theory assumes that the occurrences of spreading cells at different points in space are statistically independent (Section 1.1 in 25).Thus, if fire spread is approximated as a standard percolation process, one disregards the spatial structure of flammable vegetation.Hence, although the microscopic Staver-Levin models [10,15] consider the fine-grained patch structure, they still rely on a mean-field assumption in their implicit treatment of fire, making them prone to biases in regimes with spatial structure.To avoid these biases, microscopic models require explicit consideration of fire spread in interaction with the vegetation landscape, such as in the cellular automaton by Hébert-Dufresne et al. [16] [see also 29].In this cellular automaton, forest bistability emerges only from simple microscopic rules of vegetation and fire spread, i.e. without assuming equations or thresholds for the effects of fire.Note that larger-scale forest transitions have also been modelled with a cellular automaton, with the effects of climate and fire as spatially heterogeneous parameters [30].
In this work, we examine the spontaneous emergence of nonlinear dynamics and bistability of tropical forest from the patch-scale rules of forest and fire spread.We first use the cellular automaton of Hébert-Dufresne et al. [16] to observe the emergent structure and bistability in simulations.Next, based on the observations that forest and fire spread occur near the forest perimeter and on separated timescales, we set up a macroscopic balance equation of forest area change (Eq.9).This enables us to analyse the emergent dynamics as a function of the relevant structure, and will show that the nonlinearity is caused by the forest geometry.Then, we derive a forest resilience indicator based on our balance equation, providing a proposed link between the geometry and resilience of tropical forest.Finally, we compare our results against mean-field approximations.This will show that the assumption of absence of spatial correlations is strongly violated, particularly near the tipping threshold of forest dieback, while mean-field models still permit accurate expressions for the spatially uncorrelated regime.

RESULTS
The FGBA probabilistic cellular automaton The FGBA probabilistic cellular automaton (adapted from [16] -see Fig. 1 and Methods) models the stochastic dynamics of tropical vegetation and fire on a square lattice and in continuous time.The key empirical facts of tropical forest and fire dynamics captured by the FGBA automaton are: (i ) fires only naturally ignite in grasslands but they can spread into forest, (ii ) fires spread more easily in grassland than in forest, such that forests suppress fires, albeit imperfectly, (iii ) forest dynamics occur on a strongly separated timescale from fire spread and grass regrowth.
This results in the following reaction rules in the cellular automaton.At any time, each lattice cell can be in one of four states: F -forest, G -grass, B -burning, A -ash.Conversions between these states can occur spontaneously or due to spread to neighbouring cells (Fig. 1a and Table M1).The spontaneous conversions are: forest recruitment on grass or ash cells due to long-distance seed dispersal or from a homogeneous seed bank (G → F or A → F at rate β), forest mortality (F → G at rate γ), fire ignition on grass cells (G → B at rate ϕ), and grass regrowth on ash cells (A → G at rate λ).The conversions due to spread to neighbours are: forest recruitment due to short-distance seed dispersal on grass or ash cells (GF → FF, AF → FF at rate α), fire spread on grass (GB → BB at rate ρ g ) or on tree cells (FB → BB at rate ρ f ).Chosen parameters are in the ranges empirically justified by [16] for a square domain of 100x100 cells, with cell size ∆x=∆y=30m.The timescale separation between fire and forest dynamics implies that the rates satisfy ρ g , ρ f , µ, λ ≫ α, β, γ.In particular, we choose So, fire spreading and extinction ρ g , ρ f , µ occur on the scale of hours, while grass regrows on ash over months (λ) and forest spread, growth and mortality α, β, γ occur over decades.We take fire ignition rate ϕ ∼ 1/N such that fires spontaneously occur about once per year in the modelled area.fractions of cells in each state during a simulation with low fire ignition rate ϕ, starting from an all-grass state.Due to the low fire ignition rate, the only stable steady state is a nearly closed canopy (reached after 300 years, Fig. 1h).Before canopy closure, brief events of rapidly spreading fire counteract a gradual spread of forest.After canopy closure, fires are unable to spread.Timescale separation of forest dynamics (Fig. 1b), grass regrowth (Fig. 1c), and fire spread (Fig. 1d) shows clearly.
Figure 2a shows a bifurcation diagram of steady state forest area in the FGBA cellular automaton, denoted by [ F ] (see Eq. M1), versus fire ignition rate ϕ.Unstable steady states (saddles) were obtained by applying feedback control to the simulations (see Methods).Bistability occurs above a critical ignition rate ϕ, with alternative stable states grassland ( [ F ] ≈0) and forest ( [ F ] ≈0.83).Simulations initiated at the saddle will tip randomly up or down (Fig. 2b).Near the lower end of the bistability range, the saddle solution is fairly homogeneous, but for higher ϕ values, a single hole of grass in forest arises (insets in Fig. 2a).

Fast and slow subprocesses
The timescale separation (Eqs. 1 and 2) permits treatment of the joint vegetation and fire dynamics as a fastslow system.Fire spread occurs on the fast timescale, where the vegetation landscape is treated as constant.Forest dynamics occur on the slow timescale, where the effects of fire are a steady state function of the vegetation landscape.
a. Fast process: fires spreading in a given landscape On the timescale of a single fire event, forest dynamics are negligible (α, β, γ ≪ 1/day) such that we can consider the total landscape of forest patches as fixed.For each ignition event, this results in the following dynamics.A fire ignites on a grass cell, then spreads across its grassland cluster at a rate ρ g per BG pair, after which it reaches the interface with adjacent forest, where it starts intruding the forest at a rate ρ f per BF pair.At any time, a burning cell can stop burning spontaneously, converting to ash at a rate µ.The probabilities of fire spreading into a neighbouring grass or forest cell before the originating cell stops burning are given by where we have shown the chosen values in our simulations [adopted from 16].Since regrowth of grass and ignition of new fires occur at a much slower rate than fire spread (ϕN ≲ λ ≪ ρ f , µ, ρ g ) and our domain is relatively small (see Section S2B), we observe repeated fire spreading events well separated in time (Fig. 1b), each ending with spontaneous extinction.When a fire in grassland cluster with index j reaches its interface with adjacent forest, the resulting forest loss due to this single fire event can be approximated as (see Methods): where [ FG ] j counts the number of forest cells adjacent to grassland cluster j (with both sides of the equation optionally normalised by N ).This approximation relies on the assumptions that the fire reaches the whole interface with forest (i.e.p g → 1) and only once per fire (i.e.ρ g ≫ λ ≫ ϕN ), and that p f is small.b.Slow processes: forest demography and fire damage Forest demography and loss due to repeated fires occur on the slow timescale.Writing the number of forest-grass neighbour pairs as [ FG ] (divided by N , equivalently the total perimeter of forest or grass patches, see Eq. M1), the dynamics for tree recruitment and mortality result in an expected rate of change for [ F ] : In Eq. 5, the rates of change are β [ G ] for spontaneous forest growth on grass, γ [ F ] for spontaneous forest mortality, and α [ FG ] for spread of forest into grass at its perimeter.The rate of forest erosion at its perimeter due to fire damage over many fire events is the weighted sum over all grass clusters j=1, ..., n c , i.e.
where [ G ] j is the fraction of G cells in grass cluster j (so, [ G ] j ), ϕN [ G ] j is the rate at which fires spontaneously ignite in grass cluster j (ϕ is the rate per cell and FIG. 3. Rate of change according to Eq. 9. (a,b) Spatial snapshots at indicated times, (c,d) time series of [F], (e,f) time series of ∆ gain F , ∆ loss F , (g,h) time series of the right-hand side of Eq. 9 (gain minus loss).At t=0, the simulation is started on the saddle (on the left, [F](0) ≈ 0.2 and on the right, [F](0) ≈ 0.7, see blue and red circles in Fig. 2).Towards the left (along t + ), a simulation that tips up and towards the right (along t − ) a simulation that tips down is shown.Parameters are shown in Table M1.Columns correspond to leftmost and rightmost vertical dashed lines in Fig. 2 (ϕN =0.257 and ϕN =1.32).Domain size: 100x100 cells.
N [ G ] j is the area of the cluster), and ∆ loss F,j =p f [ FG ] j is the conversion of forest to ash caused by each fire event (see Eq. 4) (note also that [ FG ] = nc j=1 [ FG ] j ).By defining the grassland-weighted forest perimeter as the expression for forest loss becomes The grassland-weighted forest perimeter ⟨ [ FG ] ⟩ cg is the average perimeter of forest clusters weighted by the relative size of their adjacent grass cluster.

Emergent slow dynamics
We now form the balance between the slow processes discussed above, assuming fire converts trees immediately to grass (i.e.λ ≫ ϕN ).The resulting expected rate of forest area change during a short time interval is where we used Eqs. 5 and 8, and assumed on the lefthand side that N is sufficiently large, such that, via the law of large numbers, ⟨d [ F ] /dt⟩ ≈ d [ F ] /dt.Equation 9 can be understood intuitively as forest and grass competing for space within clusters (spontaneous terms) and at their interface (interaction terms).A larger interface [ FG ] leads simultaneously to faster forest spread (proportional to its perimeter [ FG ] ), and to increased exposure to fires (proportional to its grassland-weighted perimeter ⟨ [ FG ] ⟩ cg ).Fires are most damaging to forest when [ G ] forms a single cluster, i.e.
, such that each fire reaches the whole interface.Conversely, when forest patches break [ G ] into several clusters ⟨ [ FG ] ⟩ cg is smaller than [ FG ] , such that several ignitions are required to have the same effect, slowing forest erosion down.Additionally, the total amount of grass N [ G ] determines the number of ignitions and hence the rate at which grass spreads into forest.The parameters determine the relative weight of each of the discussed effects.Figure 3 shows example simulations along trajectories starting from the saddle equilibria of Fig. 2, showing forest area [ F ] in space and time (a-d), the gain/loss terms ∆ gain F and ∆ loss F defined in Eqs. 5 and 8 (e-f), and the right-hand side of Eq. 9 (gain minus loss, g-h).The left column of Fig. 3 shows simulations for low fire ignition rate ϕ and low [ F ] (0), and the right column for high ϕ and high [ F ] (0).Each column shows two realisations, both starting from the same saddle steady state.One realisation evolves toward high forest cover, shown on axis t + (increasing to the left from t=0), the other realisation evolves toward low forest cover, shown on axis t − (increasing to the right from t=0).In the stable steady states, gain (green) and loss (red) terms vary around the same mean.On the saddle (at t=0), gain and loss functions cross, indicating that the steady states and changes are accurately captured by Eq. 9.The largest changes in forest cover [ F ] occur when there are large changes in forest loss due to fire.The snapshots in Fig. 3a,b show that the high-cover state changes as an expanding/contracting hole in the forest, whereas the low-cover state is more homogeneous.

Emergent nonlinear relations
Equation 9 explains how the rate of change of [ F ] depends on the perimeter quantities [ FG ] and for three different values of ϕ, and for an ensemble of simulations starting from the saddle in Fig. 2a, with each point being a value observed at a discrete observation time.Remarkably, we observe that [ FG ] and FIG. 4. Emergent relations between perimeter quantities and forest area [F]: (a-c) forest perimeter [FG]   ϕ), which implies that [ FG ] , ⟨ [ FG ] ⟩ cg are changing on a much faster timescale, making them slaved to [ F ] .Figure 4d-f shows the terms on the right-hand side of Eq. 9 depending on [ F ] , splitting between gain and loss terms ∆ gain F , ∆ loss F , as defined in Eqs. 5 and 8. Steady states occur when gain equals loss (∆ gain F =∆ loss F ).The resulting plot of d [ F ] /dt versus [ F ] in Fig. 4g-i clearly shows the bistability of [ F ] .
Replacing the quantities [ FG ] and ⟨ [ FG ] ⟩ cg by their steady-state functions [ where [ FG ] * , ⟨ [ FG ] ⟩ * cg are functions of [ F ] and ϕ (as shown in Fig. 4a-c), and [ G ] =1 − [ F ] .With these functions [ FG ] * and ⟨ [ FG ] ⟩ * cg , the observed bistability is caused by a classic double-well potential of the gradient system Eq.10.In this ODE, nonlinearities emerge due to the equilibrium dependence of the interface on forest area (affecting [ FG ] * and ⟨ [ FG ] ⟩ * cg ), due to the segmentation of grass cells near and below the percolation threshold (affecting ⟨ [ FG ] ⟩ * cg ) and due to dependence of the ignition rate on grass patch size (multiplying ⟨ [ FG ] ⟩ * cg with [ G ] ).In Fig. S1, we show the roots of Eq. 10 using a nonparame- . These match well with the steady states obtained via control (dotdashed red).
If there is only one connected component of grass cells, we have ⟨ [ FG ] ⟩ cg = [ FG ] , such that Eq. 10 simplifies to For homogeneous initial conditions (i.e. [F ] is about the same in different large subsections of the domain), this approximation is expected to be valid for small [ F ] , where most grass cells belong to the giant connected component.Figure S1 shows the resulting steady states of Eq. 11 as a function of ϕ and [ F ] when only using the fit of [ FG ] * ( [ F ] ; ϕ) (dashed blue).The approximation is good for landscapes with low forest cover ( [ F ] ≲ 0.2).
Above [ F ] ≈ 0.2, it fails because the grassland breaks up into multiple clusters and fires are smaller than in case of a single cluster.Figure 4a-c already indicated that the single-cluster approximation is accurate for low forest cover since for low [ F ] in the scatterplots.

Resilience to perturbations
One can evaluate the right-hand side of Eq. 9 for a landscape before and after application of a small perturbation, to determine if the perturbation will be dampened or amplified under the dynamics.More precisely, we may define the sensitivity as for a landscape X and a perturbation δX, where [F ] is the right-hand side of Eq. 9. Negative values of λ F correspond to dampening (negative feedback) and positive values to amplification (positive feedback).Given that the dynamics of Eq. 9 are an approximate function of [ F ] only, the average of λ F (X+δX) over naturally expected perturbations δX (call this λF (X); see Methods, Eq.M9) can be interpreted as the approximate derivative which fully characterises the local stability of the landscape.Therefore, the sign and magnitude of λF (X) are indicators for the stability or criticality of a landscape.The dependence on only [ F ] also implies that Eq. 9 is a gradient system, such that λF (X) is the concavity of the potential energy function at X, corresponding to the classic potential-well metaphor of local resilience [14,31].Figure 5d shows λF (X) for the traversed landscapes when tipping up to the forest or down to the grassland state (same landscapes as in Fig. 3b,d,f,g).Comparison of the magnitude of λF (X) in the alternative stable states reveals that (for parameters of Fig. 3d,f,g) the grassland state is more resilient than the forest state.Figure 5a-c shows the positive feedback for a forest landscape with a ).(d) Sensitivity to perturbations (Eq.12) for all the traversed landscapes when tipping from the saddle: down to the grassland state (t − ), or, up to the forest state (t + ).The used landscapes correspond to the red vertical dashed line in Fig. 2 and the rightmost columns in Figs. 3 and 4.
hole of critical size (same landscape as shown at t=0 in Fig. 3b).Panel b shows which cells along the perimeter of the largest grass cluster contribute to the loss term ∆ loss F (red) and the gain term ∆ gain F (green).Panel a shows the effect of the perturbation obtained by converting the green cells to forest, which causes an increase in [F ] (more green in panel a).Panel c shows the effect of the perturbation converting the red cells to grass, which causes a decrease in [F ] (more red in panel c), illustrating the spatial distribution of the positive feedback.One could, in principle, also test the effect of large perturbations, and whether they will induce a transition to an alternative stable state, but it is not in general clear which perturbations are to be expected.However, in the special case of β=γ=0, when the large perturbation is a single hole in contiguous forest, only the size of the hole matters, and a simple expression for the critical size required to tip abruptly to a non-forest state can be derived (see Methods, Eq.M8).

Comparison to mean-field approximations
Our analysis of Eq. 9 enabled us to obtain macroscopic steady states and dynamics without making mean-field assumptions.Sections S3 and S3 derive a hierarchy of mean-field models for which we compare their predictions to our results to examine their validity.The simple mean field (Section S3) is unable to capture repeated fire extinction on a fast timescale and nearest-neighbour spreading of forest and fire, leading to severe bias.When instead assuming timescale separation between forest and fire dynamics and treating fire as a site percolation process in landscapes with uniform random (i.e., spatially uncorrelated) placement of forest, we obtain the meanfield approximation of Eq. 9: In Eq. 13 we substituted the a-priori unknown forest perimeter [ FG ] and the grassland-weighted perimeter ⟨ [ FG ] ⟩ cg by their expressions assuming absence of correlations: ) is the quantity given in Eq. 7 for uniformly randomly placed forest (see orange curve in Fig. S3a-c).Eq. 13 is bistable with stable high and low forest cover steady states over a wide range of parameters.It accurately predicts the location of both stable steady states (blue line in Fig. S2).Yet, its prediction of the threshold (unstable) steady state and the dynamics remains strongly biased.Indeed, away from the stable steady states, the interplay of forest demography and fire-induced forest erosion creates landscapes in which forest is spatially aggregated, strongly violating the assumption of absence of correlations.This can be seen in Fig. S3a-c.which shows that for given forest area [ F ] , the forest perimeter of simulations, [ FG ] , lies below that predicted by the mean-field, [ FG ] implying that forest is more spatially aggregated than assumed in the mean field.Aggregation results from forest spreading close to existing forest and from lower survival of forest cells that are more exposed to fire (i.e., less aggregated).Forest gain is smaller when aggregated (Fig. S3d-f) due to the smaller perimeter.Aggregation reduces forest loss at low cover while it increases forest loss at high cover (Fig. S3d-f).This is so because aggregation makes forest cells individually less exposed to fire, but collectively less effective at blocking fires, where the individual effect is dominant at low cover and the collective effect is dominant at cover values near and above the fire percolation threshold.
For our choice of parameters, the stable steady states and the lower saddle-node bifurcation contain negligible spatial structure, such that mean-field predictions are accurate.We derive their expressions below for β≈0 .
a. Stable steady states By Eq. 13, the low-cover steady state [ F ] * − has to be approximately zero for β≈0.At high forest cover, loss due to fire is negligible, such that the high-cover steady state [ F ] * + can also be obtained from Eq. 13 (using β≈0): For the chosen parameters, we have [ F ] * + =0.83, which is in excellent agreement with simulations (Fig. 2).
b. Onset of bistability At low forest cover, grass consists of a single cluster, such that we can write in Eq. 13 An expression for the lower bifurcation point can be obtained by finding the root of the derivative of the righthand side of Eq. 16 with respect to [ Rearranging this relation for ϕN , we can obtain the fire ignition rate above which tropical forests are bistable with grasslands: For the chosen parameters, (ϕN ) min =0.25, which agrees well with simulations (Fig. 2) and which corresponds to a maximum fire return interval of (ϕN ) −1 min =4y.

DISCUSSION
In this paper, we showed how nonlinear dynamics and bistability of tropical forest emerge spontaneously from the patch-scale rules of forest and fire spread, without assuming equations or thresholds for the effects of fire as in previous work.Below, we first summarise our main results on structure and dynamics.Then we discuss the importance of the emergent structure as indicated by comparison with mean-field approximations.Finally, we highlight the potential practical implications of our results for resilience assessment and conservation.
a. Emergent structure and dynamics Our simulations showed that spatial structure emerges due to forest expansion and fire-induced damage at the forest perimeter.As a consequence, the forest perimeter appeared in both the gain and loss side of our landscape-scale balance equation of forest area, Eq. 9, where losses require weighting by adjacent grassland area.Remarkably, when plotting the changes predicted by our balance equation versus forest area, using landscapes from the simulations, we found that they lie on an approximate curve (Fig. 4gh).As this curve shows the change of forest area as a function of forest area, this means that the emergent macroscopic dynamics can be described by a simple ordinary differential equation, Eq. 10.In this emergent closed form of our balance equation, the perimeter quantities determine the nonlinearities.Therefore, Eq. 10 elucidates how forest dynamics and bistability are linked to the forest geometry that emerges from the patch-scale spreading rules.Note that, as in previous work, Eq. 10 does not include fire explicitly, because it does not contain equations for fire.This follows from timescale separation between fire and vegetation dynamics, an assumption that was already implicit in mean-field models (12,(17)(18)(19)21; mean field in 10, 15) and microscopic models (microscopic models in 10, 15) focusing on alternative stable states.However, previous work derived the implicit effect of fire in closed form by relying on standard percolation theory, which assumes that occurrence of flammable patches is spatially uncorrelated [25].As we did not rely on percolation theory, but observed the closed form emerging in simulations (Fig. 4), we could avoid the biases that affect previous work.
b. Evaluation of mean-field models We compared mean-field models against the emergent closed form of our balance equation to assess their validity (Fig. S3) and to show where spatial structure is important.This showed that mean-field models are in qualitative but not quantitative agreement with simulations: existence of bistability, but not its parameter range is robust to mean-field assumptions.In particular, the simple mean field (Eq.S7) is strongly biased due to its failure to account for two phenomena that are present in the microscale model: (i ) spontaneous fire extinction on the fast timescale, leading to separated rapid fire spreading events, (ii ) nearest-neighbour spreading of fire and forest, leading to emergent aggregation of forest patches away from the steady states.The former violates the meanfield assumption of large system size (N →∞) and the latter that of absence of correlations.That spatial structure affects steady states and dynamics is well known [e.g.22,24].Even when addressing timescale separation and using results from percolation theory for the effect of fire (Eq.S18 or Eq.S22), a large bias remains except near the alternative stable states Fig. S3.This is because standard percolation theory only considers lattice configurations with uniform random (i.e., spatially uncorrelated) placement of flammable sites while our forest-grass landscapes are shaped by past fires and vegetation dynamics.As Fig. 3a,b (at t=0) shows, forest aggregation is particularly strong at the tipping threshold for forest collapse, implying that mean-field models cannot be used to study abrupt forest dieback.Despite their severe bias concerning forest dynamics and tipping, mean-field models are still useful for studying regimes with little structure, such as near the stable equilibria or for dynamics with uniform seed dispersal.This enabled us to derive expressions for these equilibria (Eq.15) and the point of onset of bistability (Eq.17).The latter result was not obtained by previous mean-field models because they did not include parameters that relate directly to fire [10, 12, see Section S7 for a suggested modification] or they did not account for timescale separation in finite domains [16].

c. Implications for resilience assessment and conservation
The link between geometry and dynamics implies that tropical forest resilience can be empirically estimated from its spatial structure.The spatial structure, as captured by the perimeter quantities [ FG ] and ⟨ [ FG ] ⟩ cg , can hence be treated as a measurable parameter additional to the microscopic parameters.Microscopic parameters (given in Table M1) can be inferred from remote-sensed data (as in 16 or 30) or from experiments [as for fire spread in 11], while the perimeter quantities can be calculated for any observed landscape.In regimes with negligible spatial structure, one can assess stability or resilience from the microscopic parameters alone, based on mean-field results.E.g., if the onset point of bistability at low tree cover lies in the regime without spatial structure, as in our simulations, one can directly estimate the minimum fire ignition rate for onset of bistability from the microscopic parameters (Eq.17).This expression then shows which natural or abandoned degraded areas of low tree cover with fire ignition rate beyond this point will not spontaneously recover to closed tropical forest.In regimes with spatial structure, the mean field is highly inaccurate (Fig. S3), such that spatial structure needs to be considered in addition to the parameters.In particular, in our simulations, the tipping threshold obtains spatial structure at higher fire ignition rates (Fig. 3a,b at t=0) and approaches the stable forest equilibrium much more closely than in the mean field (Fig. S2).While this makes the mean field unsuitable for studying forest resilience and dieback, our balance equation (Eq.9) does not have this limitation because it makes no assumption on spatial structure.We demonstrated how Eq. 9 permits estimation of the resilience of a landscape to perturbations, via λ F (Eq. 12).In contrast to generic indicators of resilience [31,32], λ F is an indicator that can be obtained from a single landscape, and for which the contribution of each relevant spatial process can be examined.Furthermore, landscape perturbations by human intervention can be evaluated, through sensitivity λ F , for how they will amplify or mitigate fire-vegetation feedback.Forest conservation/restoration may introduce targeted perturbations that most efficiently prevent resilience loss of high-cover states or induce resilience loss of low-cover states.For instance, in Fig. 5a, forest dieback is averted by a perturbation that divides the largest grass cluster into smaller ones.It may thus be anticipated that maintenance or creation of barriers to fire spread will be essential here.

METHODS
a. Details of the FGBA probabilistic cellular automaton The FGBA probabilistic cellular automaton is a minimal spatial stochastic process that models the joint dynamics of tropical vegetation and fire.It is an adapted version of the BGT(A) model of ref. [16].The modifications compared to ref. [16] are: (i ) it runs in continuous time, (ii ) it includes a spontaneous forest mortality rate γ, (iii ) species T is labelled as F, consistent with other models of tropical vegetation dynamics [12,15,18].Note that according to some definitions, probabilistic continuous-time cellular automata are considered as interacting particle systems.In general, when studying the stochastic dynamics of a number n of interacting species on a square lattice with N cells, the state of the system can be represented as where X i is the label of the species that occupies cell i.Each cell is occupied by exactly one of four possible species: grass, forest, burning and ash, with labels G, F, B and A. Transitions between states (species) occur in continuous time, resulting in a continuous-time Markov chain with a state space of size n N .The reaction rules for transitions between states are shown in Table M1, where spontaneous conversions are shown on the left and conversions due to nearest-neighbour interactions on the right (see also Fig. 1a).
The latter type of interaction occurs over each four nearest neighbour connections of the indicated type.E.g., fire will spread into a given grass cell with a rate ρ g for each burning neighbour.For realistic timescales, our parameters satisfy Eqs. 1 and 2, which were empirically justified in [16].We borrow our notation from the moment closure literature [e.g.[22][23][24]41], writing the global fraction of species with label x and the interface between species with label x with label y respectively as where both are normalised by N , δ is the Kronecker delta function (δ x (y)=1 if y=x and 0 otherwise) and A ∈ {0, 1} N ×N the adjacency matrix.We simulated the cellular automaton via a Gillespie algorithm [42] and used a domain of N =100×100 (N =200×200 in Fig. 1) cells with periodic boundary conditions.
b. Noninvasive feedback control To study steady states regardless of their stability in a simulation, we apply noninvasive feedback control [43][44][45][46].To obtain the dependence of equilibria of [ F ] on fire ignition rate ϕ, we introduced an artificial stabilizing feedback loop of the form The factor g is called the feedback control gain and is problem specific.The property of noninvasiveness means FIG.M1.Feedback control applied to the CA without spontaneous mortality (γ=0): (a) the unstable steady state of the bifurcation diagram (dashed) was derived via feedback control by letting ϕ be a function of [F] (blue line) such that it is stabilised, then obtaining (ϕ, [F]) by averaging and repeating this for many [F] ref (with appropriate g), (b) a regular simulation with the same ϕ value (solid grey in (a)) and starting from the final state of the controlled simulation tips up or down depending on initial perturbations, (c) snapshot of the domain at the saddle for the control indicated in (a) (black: forest, white: grass).For other parameters, see Table M1.
that the controlled simulations have the same equilibria as regular simulations [47][48][49].This implies that if one extracts the equilibrium values of the controlled simulation (ϕ * , [ F ] * ), one can use them to plot a 1-parameter bifurcation diagram of the simulation without control.Fig. M1 shows the control graphically.The feedback Eq.M2, indicated in blue in Fig. M1a, stabilises a steady state that is unstable in a regular simulation.This can be seen in Fig. M1b, where the unstable steady state is first stabilised via control, after which the control is removed and a regular simulation is started with the effective rate and the landscape (inset) obtained from the controlled simulation.Depending on initial perturbations, the regular simulation gets either attracted to the 100% forest state or to the low tree cover state.When the controlled simulation is in equilibrium (steady part of the blue curve in Fig. M1b), the steady state values of [ F ] are obtained via taking the time average, i.e.
[ F ] = 1 where t 0 is the time after which the dynamics have settled to a steady state and T the averaging time.If n G→B is the number of ignition events between t = t 0 and t = t 0 + T , the steady state of ϕ is obtained by calculating the mean ignition rate as n G→B /T and dividing this by the mean number of grass cells, such that where [ G ] is obtained as in Eq.M3.When repeating this exercise for many [ F ] ref values, one can get multiple points on the unstable branch.Points on the stable branches can be obtained with regular simulations.On the final selection of points, we applied Gaussian process regression to obtain smooth curves and used moving block bootstrapping [50] to obtain confidence intervals.One of the advantages of applying control is that one can obtain states for which one would have to wait prohibitively long in a regular simulation due to their instability.c.Forest loss due to a single fire A fire in grassland cluster with index j that reaches its interface with adjacent forest induces a forest loss that can be approximated as follows.Consider a single forest cell i located at the interface with grassland cluster j with [ FG ] i,j number of neighbouring grass cells.When assuming that spreading events are independent, the probability that the forest cell gets burnt is the complement of the probability that none of its neighbouring grass cells in cluster j spread the fire to the forest cell: where the approximation on the right is valid for small p f .Summing over all forest cells at the interface of grassland cluster j, we obtain the expected loss of forest per fire event as shown in Eq. 4: This approximation also assumes that burning forest cells at the interface do not spread the fire further, which also relies on p f being small.For an evaluation of the validity of this approximation in case of landscapes without spatial structure, see Fig. S7.d.Critical hole size for an abrupt shift when γ=β=0 When there are no spontaneous transitions and we perturb a fully closed forest of 100% cover by creating a hole with grassland, an expression can be obtained for the critical hole size beyond which fire causes an abrupt shift to grassland.Using that grassland is a single cluster, Eq. 9 becomes which has two absorbing steady states [ F ] * l =0 and [ F ] * h =1, and an unstable steady state at [ F ] * c =1 − α ϕN p f .The critical hole size is then the complement of the unstable steady state: which can also be written as [ G ] c = ϕ 1 /ϕ, where ϕ 1 is the value of ϕ for which [ G ] c = 1, or also the lower limit of the bistablity range.e. Sensitivity to perturbations We estimated λF in Fig. 5d for a given landscape by averaging Eq. 12 over realisations of different types of perturbations: (Eq.13) (Eq. 10) ,mf ,mf ,mf ,mf FIG.S3.Emergent relations between key quantities and forest area [F] compared to the mean field with site percolation (dots: simulations, lines: corresponding mean field quantities from Section S4B): (a-c) forest perimeter [FG]

S2. RELEVANT CHARACTERISTICS OF THE FIRE SPREADING PROCESS
Before obtaining the mean-field equations for coupled vegetation and fire dynamics, we analyse the fire spreading process in isolation.The insights from this section will enable us to set up a mean-field model that constitutes the fairest comparison against the analysis in the main text.

A. Definition and mean field
When we remove state F and its conversion rates to/from other types (α, β, γ, ρ f ) from the FGBA process, the dynamics show fire spread alone.We call this the GBA process.Writing x i as shorthand for δ x (X i ) (equalling 1 if X i = x and 0 otherwise) and taking expectations in each cell i, we obtain equations for the rate of change of the expectation that cell i is occupied by species x ∈ {G, B, A}, where ⟨•⟩ are ensemble averages, [ x ] the domain fraction of species x, and [ xy ] the total number of neighbouring xy pairs divided by N , later referred to as the xy interface or xy perimeter.This set of equations can be derived rigorously from the master equation [e.g.51].To go from individual (left) to population level (right), we summed over i and divided by N , using Eq.M1.Equation S1 is not a closed system.To close the system, we need to determine all undetermined terms [ xy ] on the right-hand side without creating new unknowns.The simplest way to do this is to assume absence of pairwise correlations, i.e.
We take the additional assumption of N → ∞, such that the law of large numbers applies and [ x ] → ⟨ [ x ] ⟩.These assumptions are valid when all cells in an large domain interact with each other at uniform contact rates of order 1/N .This results in the simple mean-field approximation of the GBA process: where we also used the dot notation for time derivatives.Substituting [ and taking only the independent equations, we finally obtain We further focus on the case ϕ = 0, the reason for which will become clear below.When ϕ = 0, Eq.S3 has two steady states, a trivial one at 1+µ/λ ).The eigenvalues of the Jacobian of Eq.S3 show that for 4ρ g /µ > 1, the trivial steady states is a saddle and the non-trivial a spiral sink.For 4ρ g /µ < 1, the trivial state state is a stable node and the only physical solution.Hence, the steady states exchange stability at the transcritical bifurcation at 4ρ g /µ = 1.The GBA process with ϕ = 0 is equivalent to the SIRS spreading process in epidemiology [41], which represents spread of a disease in a population with waning immunity.Fire B plays the role of infected individuals.Infections spread through a population of susceptibles G at rate ρ g per GB link.They subsequently acquire a state of immunity A at rate µ, which can be lost at rate λ.The non-trivial steady state corresponds to the endemic equilibrium and the transcritical bifurcation to the epidemic threshold R 0 .A well-known characteristic of this spreading process with ϕ = 0 and [ B ] 0 > 0 is that in finite systems, it goes extinct in finite time, even when R 0 > 1 [52].This is so because stochastic excursions away from the non-trivial equilibrium will eventually reach the absorbing trivial state.When the spontaneous ignition rate ϕ > 0 and the time to extinction is much smaller than the typical waiting time between ignition events, there are repeated fire events separated by extinction events.The dynamics then effectively behave as a series of GBA processes with ϕ = 0 and [ B ] 0 > 0. This is what we observe in cellular automaton simulations on a square lattice of 100×100 cells for realistic parameters (Fig. S4b, right panel).The mean-field approximation (Eq.S3), on the other hand, does not show extinction due to its assumption of N → ∞.Instead, it shows a single pulse (Fig. S4a, left panel) after which a high-ash low-grass and non-zero fire steady state (the endemic equilibrium) is reached (Fig. S4a, right panel).In the case ϕ = 0, the required lattice size to avoid extinction with high probability depends on the initial conditions [52], but for realistic parameter ranges, it is unrealistically large.This can be understood as follows.

B. Extinction in finite systems
• When the initial condition is a single fire, at least one cell has to keep on burning until the density of grass has regrown to a level sufficient for a new wave to propagate.This translates into the condition (L/∆x) 2 exp(−µ/λ) ≥ O(1), such that L ≥ O(10 4•10 4 ) for our parameters (taking a grid size of ∆x = 0.03km as in [16]).
• When initial conditions are such that a band of the domain is immune at the start, a single fire can keep on burning by crossing the domain repeatedly [52].When assuming ρ g ≫ µ and using that waiting times between spreading events are exponentially distributed with mean 1/ρ g , a fire will spread throughout the domain in a time of the order τ ≈ L/(ρ g ∆x).For there to be sufficient regrowth of grass on this time scale, we need L/(ρ g ∆x) ≈ 1/λ, or L ≈ ρ g ∆x/λ.For the parameters we have chosen, this means L = O(10 4 )km, i.e. the order of magnitude of the earth's circumference, which is drastically smaller than the above estimate yet still impractically large.
Taking more conservative estimates for fire spreading rates or taking account of a small positive fire ignition rate ϕ = O(λ/N ) for the initial condition with a single burning cell, this may be decreased by an order of magnitude, i.e. the size of a continent or country.Still, in reality, extinction will occur on smaller scales due to spatiotemporal heterogeneity of forcing parameters as a consequence of climatic seasonality or existence of natural or artificial boundaries (such as forests), leading to a lower effective system size.Hence, in any real system, repeated extinction and system-scale oscillations are to be expected.

C. Percolation analysis of a single fire event
Therefore, a single fire in realistically sized systems corresponds to the case ϕ = 0, starting with a single burning cell.Using that the regrowth of grass occurs on a much slower time scale, we can further also set λ = 0 in our following analysis.The GBA process with ϕ, λ = 0 is equivalent to susceptible-infected-recovered (SIR) epidemic spreading [41].The final size of the epidemic in SIR epidemic spreading on a lattice shows a continuous phase transition (CPT) at a critical spreading rate ρ g and scaling laws near the critical point obey those of the ordinary percolation universality class [53].Figure S5 shows mean quantities for SIR epidemic spreading on a square lattice in a range of infection probabilities and initial number of immune individuals, which are spatially uniformly distributed.In particular, we show that SIR epidemic spreading is a type of mixed site-bond percolation, with bond occupation probability given by p b := p g = ρ g /(ρ g + µ) (which is fixed at 0.9 in the main text, as in ref. [16]) and site occupation probability given by p s := [ G ] 0 , i.e. the initial fraction of cells that are grass in fire spreading, or the complement of the initial fraction of immune individuals in epidemic spreading (with the rest being susceptible).In Fig. S5, we record the mean cumulative probability of being burnt ⟨Q⟩ and the susceptibility χ of ⟨Q⟩, defined as where [ • ] 0 denotes initial value and [ • ] * final value.For N →∞, ⟨Q⟩ converges to the percolation probability P ∞ , which is the probability that a grass cell belongs to the giant connected component.We use the location where χ peaks as an estimate of the percolation threshold [25,56].When p s = [ G ] 0 = 1, we have pure bond percolation and when ρ g /µ → ∞, we have pure site percolation.The percolation threshold for standard mixed site-bond percolation (from [54]) is shown in Fig. S5 with a dot-dashed blue curve.The pure bond percolation threshold (when [ G ] 0 = 1) of SIR epidemic spreading occurs at higher p b than in standard bond percolation (p b ≈ 0.538 > 0.5) because the possibility of spreading to multiple neighbours makes the bond occupation probability spatially autocorrelated, as shown by [53].The pure site percolation limit shows the classical value (for the square lattice) of [ G ] 0 ≈ 0.593.From Eq. S3 (with ϕ = λ = 0), we can obtain the mean-field percolation threshold p s,mf by finding where the trivial state becomes unstable in Eq.S3, which is given by As shown in Fig. S5c, the mean-field approximation shows a large bias towards lower values.Eq.S6 (dashed black).The dash-dotted blue line indicates the location of the infinite-size percolation threshold for uncorrelated mixed site-bond percolation (taken from [54]).The GBA model's percolation threshold lies at higher values (b) due to spatial correlation of pg as explained in the text.The colour scale was taken from [55].See Fig. S11 for more detail.
a. Implications for the FGBA process On landscapes with forest, fires can be blocked (albeit imperfectly) by forest cells.These landscapes obtain a steady state shape due to the shaping processes of forest demography and fire.Hence, results from the spatially uniform [ G ] 0 above do not apply to percolation effects in the full FGBA process.That is, when fire spreads on landscapes with forest, the critical point for pure site percolation (ρ g /µ → ∞, or p g → 1) will in general depend not only on the site occupation probability [ G ] 0 but also on the spatial correlation function of site occupation.For the idealised case of fire-proof forest (p f = 0), fire percolation on real landscapes is then equivalent to correlated mixed percolation, where correlations in bond occupation probability occur due to the spreading process, and correlations in site occupation probability occur due to the nonrandom spatial structure of the landscape.When p f > 0, the spreading process becomes a heterogeneous (correlated) bond percolation processes, i.e. a percolation process in which fire spread on grass occurs with bond occupation probability p g and on forest with bond occupation probability p f .The possibility of spreading on forest decreases the percolation thresholds compared to the correlated mixed percolation limit of p f → 0. This decrease is expected to be small because forests do not spread fires well (p f ≈ 0).

S3. SIMPLE MEAN FIELD OF JOINT FOREST AND FIRE SPREAD
When we follow the same steps as in S2A, we obtain the simple mean-field approximation of the FGBA process: The simple mean field shows bistability of tree cover (first shown by ref. [16] for γ = 0) in ranges of all parameters that are expected to show considerable spatial heterogeneity in a given ecosystem: α, β, γ, ϕ (Fig. S6).However, despite being qualitatively correct, it shows a large bias compared to simulations.For the parameter ranges of our simulations, it has no non-trivial solution for positive fire ignition rate ϕ, so we had to choose different parameter values to find its bistability range.This bias is due to the inability of the simple mean to capture two effects: repeated fire extinction on a fast timescale and the spatial nature of the two spreading processes.In the following section, we derive alternative mean-field models that partially correct for these biases.

S4. TWO-TIMESCALE MEAN FIELD OF JOINT FOREST AND FIRE SPREAD
Here, we derive an alternative mean-field model that takes account of separation of fire events in systems with realistic sizes, assuming that fire spread occurs on a much faster time scale than forest spread.This means that we FIG.S6.Steady states and bifurcations of forest cover [F] in the simple mean field of the FGBA process (blue: steady state manifold, green/orange contours at fixed axis values, red: saddle-node bifurcations): (a) versus fire ignition rate ϕ and forest spreading rate α, (b) versus fire ignition rate ϕ and spontaneous forest growth rate β, (c) versus fire ignition rate ϕ and spontaneous forest mortality rate γ.Due to its large bias, the simple mean field shows different bistability ranges than the simulations.We set pg=0.25 (requiring a ρg that is 27x smaller than in simulations) such that bistability ranges are visible (remaining parameters are as in Table M1).
can consider the fire spreading process in isolation with ϕ = λ = 0 (as argued in Section S2), and take the asymptotic amount of forest burnt by a single fire before extinction on the fast time scale as forest mortality per fire event on the slow time scale.

A. Well-mixed fire and forest
We start with the simplest case, where both vegetation and fire mix uniformly, which is one way to conform with the mean-field assumption of absence of correlations.
1. Fast process: forest loss due to a single fire On the fast time scale, we can set all small parameters related to forest demography, grass regrowth and fire ignition to zero (α=β=γ=λ=ϕ=0), such that we obtain where the products arise from the well-mixedness assumption as before.By rewriting the equations for we can obtain [ G ] and [ F ] as a function of [ A ] via separation of variables and integration: Substituting Eq.S10 into the equation for d dt [ A ] in Eq.S8 and setting the time derivative to zero, we obtain an implicit relation of the asymptotic amount of vegetation burnt: When taking an initial state consisting of only grass and forest, that is * can be found numerically as a function of ρ g /µ, ρ f /µ and [ F ] 0 .This can in turn be used to obtain the total amount of forest lost due to a single fire, from Eq. S10, where subscript mf denotes mean field.A plot of Equation S12 versus [ F ] is given in Fig. S7 (solid purple curve), which shows that the percolation threshold -where forest starts blocking fire -lies at unrealistically low grass cover, as was also the case for perfectly blocking forest (Fig. S5c, dashed line).

Slow processes: forest demography and fire damage
Now, we can define the mean-field forest gain and loss terms on the slow time scale as such that the final mean-field model becomes The steady states are shown in Fig. S9 in solid purple for the same parameters as those used in the simulations of the main text.For low and high tree cover, it reproduces the steady states fairly accurately, but unlike the simulations, it shows no wide saddle in between.Hence, while this is an improvement compared to the simple mean field, there is still a large bias at intermediate tree cover.To address this bias, we need drop the assumption of uniform mixing for fire spread. ( FIG. S7.Grass burning probabilities and forest loss per fire in landscapes without spatial structure: (a) probability that a grass cell burns ⟨Qp g ⟩, for pg = 1 ('+') and for pg = 0.9 ('•'), together with fits to logistic functions; (b) loss per fire estimated from grassland-weighted forest ('×', Eq.S16), from ⟨Q0.9⟩ ('×', Eq.S20), from ⟨Q1⟩ ('+', Eq.S20), by assuming uniform mixing (purple curve, Eq.S12), and measured in fire simulations with uniform random placement of forest ('•').
B. Spatial fire percolation and uniformly randomly placed forest While the uniform mixing assumption may be ecologically justified for forest spread in case of species with longrange seed dispersal, it is much harder to justify for fire spread, which is fundamentally a local contagion process.Therefore, we aim to take into account the effects of fire as a percolation process while still assuming absence of spatial correlations between forest cells.Because the percolation process affects forest loss, this only affects the loss function.Earlier mean-field models [10,12,15] accounted for the effects of fire percolation by making fire-affected rates threshold functions of tree cover, such as those shown in Fig. S7a, while assuming that vegetation remains spatially uncorrelated [for derivation, see 10].Therefore the mean-field analyses presented here provide the fairest points of comparison.
To estimate the loss due to fire, we will show two alternative approaches.The first is equivalent to our approach in the main text, using grassland-weighted forest perimeter to estimate exposed forest.The second estimates exposed forest via standard results from percolation theory, which are valid here due to the assumption of uniform random placement of forest.
FIG. S8.Steady states and bifurcations of forest cover [F] in the two-timescale mean field of the FGBA process via Eq.S18 (blue: steady state manifold, green/orange contours at fixed axis values, red: saddle-node bifurcations) as a function of (scaled) fire ignition rate ϕN and: (a) forest spreading rate α, (b) spontaneous forest growth rate β, (c) spontaneous forest mortality rate γ.Parameters others than the ones on the axes are the same as those chosen in simulations -see Table M1.
1. Using the grassland-weighted forest perimeter ⟨ [ FG ] ⟩ cg (see Eq. 7).According to this approach, fires spread perfectly to the forest perimeter, where a fraction of the forest is burnt.The difference with the main text is that the landscapes in which fire spreads have uniform random placement of forest.We indicate this difference below by the superscript u in ⟨ [ FG ] ⟩ u cg .The resulting loss per fire is such that the loss function is where subscript pu refers to percolation on a square lattice with uniform random placement of forest.The final mean-field model is where we made explicit that ⟨ [ FG ] ⟩ u cg is a function of [ F ] .This mean-field model shows clear bistability for the same parameter ranges as in simulations (Fig. S8).For the exact same parameters, this mean-field is qualitatively most comparable to simulations, but the saddle is much flatter (Fig. S9, solid blue line versus black dots).
2. Using percolation theory.Alternatively, we can estimate forest loss using mean burning probabilities from percolation due to a single fire.The forest perimeter exposed to fire for a given realisation is the interface of burnt grass with forest at the end of the fire [ FA ] * .When we take the expectation (for given total grass cover) and forest cells are assumed to be uniformly randomly placed, we have where ⟨Q pg ⟩ is the mean proportion of grass that burns for given p g and [ G ] (see Eq. S4; shown in Fig. S5b).
The loss per fire in this case is then such that the loss function is (multiplying by ϕN [ G ] ) The final mean-field model is then This mean-field model has very similar steady states as Eq.S18 but the saddle is slightly lower (dotted blue line in Fig. S9).In the limit of large domain size, ⟨Q pg ⟩ may be replaced by the percolation probability P ∞ , which is defined as the probability that a grass cell belongs to the giant component [26] (see also Fig. S11a).
Both of the estimates above assume that p f = 0 for fire spread and that p f is small for loss of (uniformly randomly placed) forest due to fire.The first further assumes that fire spreads perfectly on grass (p g = 1).Therefore the two estimates are equivalent when the spreading process is pure site percolation, for which (equating Eq.S17 and Eq.S21), which is confirmed by Fig. S7b ('×' and '+' symbols).Comparing the estimates to recorded forest loss in fire simulations where only the assumption of random placement is taken ('•' in Fig. S7b), one sees that the second estimate ('•' in Fig. S7b) is more accurate than the first estimate ('×' in Fig. S7b), despite that it carries more assumptions.Hence, the error due to the assumption of grass perfectly spreading compensates the error by the assumption of forest perfectly blocking fire.We expect that the difference between the two approaches will be smaller for landscapes with spatial aggregation of forest, where fires spread in pockets of high grass cover, for which ⟨Q 1 ⟩ − ⟨Q 0.9 ⟩ is smaller (Fig. S7a).

C. Detailed comparison against simulations
Here, we compare the time-separated mean-field models to the simulations of the main text and also to other simulations with different spreading ranges for forest and/or fire.In Fig. S9, dots with error bars are simulations and lines are mean-field approximations.The black dots are the simulations from the main text, with nearest-neighbour spreading for both fire and forest.Purple dots are simulations where both fire and forest can spread to any other cell.Blue dots are simulations where forest can spread to any other cell, but fire spreads along nearest neighbours.And finally, orange dots denote simulations where forest spread occurs in a Gaussian neighbourhood with standard deviation 60m (two cells).
All simulations and approximations agree fairly well on the parameter value where the lower saddle-node bifurcation occurs.All except the fully well-mixed case agree on the stable steady states.The disagreement occurs particularly for the unstable steady states, where the effect of spatial structure is hence most pronounced and process-dependent.There is little or no bistability in the well-mixed two-timescale mean-field (purple line) but its good agreement with uniformly mixed simulations (purple dots) further indicates the validity of the assumption of time-scale separation with well-separated fire events.The mean field where fire is a site percolation process on landscapes with uniform random placement of forest (blue line), is qualitatively more correct but it has a much flatter saddle than the simulations (black dots).Hence, even the mean-field model that takes into account the effects of percolation while keeping forest cells spatially uncorrelated remains strongly biased due to the importance of spatial aggregation of forest cells.Taking larger neighbourhood sizes in the simulations does not change this (orange dots).Even compared to simulations with uniform forest dispersal (blue dots), the mean field with site percolation shows some bias, indicating that the fire spreading alone already induces some spatial structure.This is most likely caused by lower survival rates of solitary compared to aggregated forest cells.
The bias of the mean field is even more apparent in the dynamics.Figure S3 shows the same figure as Fig. 4, but with the corresponding values of the mean field with site percolation on top.Panels a-c show large differences between [ FG ] and ⟨ [ FG ] ⟩ cg of simulations (scattered dots) versus those from the mean field (curves).In particular, while for the mean field, the perimeter is the parabola [ ), the perimeter of simulations lies below this parabola for any [ F ] .That the simulated perimeter is lower for given forest area means that forest is more spatially aggregated in simulations.This results in lower forest growth rate at any cover value (panels d-f) because fewer forest cells can expand into grass.Fire-induced damage is lower below [ F ] ≈ 0.4 and higher above (panels d-f).This is so because damage per fire (at given cover) is determined by two effects: exposure of forest and clustering of grass.Below [ F ] ≈ 0.4, there is no clustering, such that only decreased exposure due to aggregation can decrease forest loss.Above [ F ] ≈ 0.4, aggregation decreases clustering, such that grassland stays fully connected at higher forest cover than in the case with uniform random placement, with larger fires as a consequence.This further leads to an upward shift of the unstable forest state compared to the mean field (panels (g-i), see also Fig. S9).The effect of forest aggregation on fire spread has an equivalent in disease spread: in the SIR process, aggregation of immune individuals lowers the epidemic threshold, such that it elevates the population immunisation threshold to eradicate the epidemic [24].Note though, that, as argued above, the equivalent epidemic process to tropical fire spread in forest-grassland landscapes is not the regular SIR process, but one that has a mix of two populations: susceptibles (grass) and imperfectly immunised individuals (forest).

S5. EVOLUTION OF FRONTS -HETEROGENEOUS STATES
Here, we illustrate the case where grass and forest are initially separated into two contiguous areas with their interface extending along a straight line.Because for this type of initial conditions, the single-cluster approximation (Eq.11) is valid, we can focus on the evolution of the interface.As spontaneous conversion between forest and grass (with rates β and γ) increases independence between cells and promotes homogeneity at large scale, we expect the effects of heterogeneous initial conditions to be most persistent when the spontaneous conversion rates β and γ are small.Therefore, we will set β = γ = 0, for which Eq. 9 becomes Hence, the precise shape of the interface [ FG ] does not affect the location of the steady states, only the rate at which they are approached or receded from.The trivial steady states of Eq.S23 are [ F ] = 0 and [ F ] = 1 (where [ FG ] = 0), which are stable, and between them, there is the saddle As seen in Fig. S12, this analytical prediction (solid black) matches the controlled simulations with p g = 0.9999 (shaded blue).For p g = 0.9 (shaded red), which we used before, there is a small bias.In the limit of N → ∞, Eq.S24 converges to [ F ] * ∞ = 1, implying that in an infinite domain, any positive fire rate leads to extinction of forest below [ F ] * = 1.When ignoring the effect of ash, this would also occur for heterogeneous initial conditions.That is, considering an infinite domain with many grass clusters of which the size is a random variable (with support [0, ∞)), there will be initial grass clusters of arbitrarily large size, which will expand and eventually drive forest extinct.However, such determinism does not occur in the simulations because at high fire rates, patches with ash start to block fires, and the rate of exposure of the forest interface to fire becomes limited by the rate at which ash is converted back into grass.As a first correction for this, one can multiply p f with the average proportion of grass sites that are in the ash state after the expected waiting time between fires 1 − exp(−λ/ϕN ), such that [ F ] * ∞ = 1 − α/λp f .Keeping in mind that we are focusing on heterogeneous states, the analysis here implies that for γ = β = 0, there is a critical patch size above which the forest patch expands and below which it contracts.The intuition is that above the critical forest patch size, there is not enough grass area to reach the minimum number of ignitions required to erode the forest.

S6. A NOTE ON FINITE SIZES AND BISTABILITY
In simple bistable chemical systems, it is known that bistability converges to an abrupt phase transition in the thermodynamic limit (N →∞) [57,58], at a value known as the Maxwell point [e.g.59], making the macroscopic state of the system deterministically dependent on the parameters (e.g.pressure or temperature).With only forest and grass or provided that savanna and forest components are sufficiently decoupled, the same behaviour occurs in spatial mean-field models of tropical forest, with a front between forest and nonforest that is depends deterministically on environmental drivers [18].We do not expect such determinism as N →∞ to arise in the FGBA cellular automaton.The infinite FGBA cellular automaton possesses grass clusters of arbitrary size, such that, even when assuming that fire spreads instantly on grass, there will always be some parts of the forest perimeter shielded from intruding fires by adjacent ash cells.Were it not for this shielding effect, then there would be a deterministic dependence of the dynamics on fire ignition rate away from the absorbing states: ϕ=0 would lead to forest spread while ϕ>0 would lead to forest extinction (see Section S5, for β=γ=0).In reality, finite fire spreading rates and, in particular, the effects of heterogeneity in space or time impose stronger limits on the reach of fires.
In finite domains, both the cellular automaton and scalar reaction-diffusion equations (with bistable reaction term) show bistability due to critical patch sizes or domain shapes, and dependence on interfacial characteristics [e.g.19,60,61], but this correspondence requires further scrutiny.In realistic scenarios, we then suspect that the amount of bistability depends (besides the parameters) on the ratio of the characteristic interaction scale and the scale of observation.The range of interaction in turn depends on e.g.fire spreading and/or plant dispersal ranges.E.g., if we take as interaction scale the typical size of a savanna fire (assuming that it exists) O(10 1...2 km 2 ) [62], this corresponds to an area in the cellular automaton of 100×100 to 300×300 cells.This also corresponds to our observation scale (domain size) in the main text.Hence, it may be that the bistability observed in our work is a finite-size effect, and that a larger observation scale leads to more gradual transitions due to existence of multiple stable patterns [28,40].

S7. HOW TO INCLUDE FIRE PARAMETERS IN EXISTING MEAN-FIELD MODELS
Previously derived mean-field models of tropical tree cover bistability did not include parameters that relate to fire.Here we give suggestions on how to include fire ignition rate and the appropriate percolation quantities, focusing on the Staver-Levin mean-field model [10,12] of tropical savanna and/or forest bistability.We assume the reader is familiar with [10,12].
By running an infection process on clusters obtained by standard site percolation, [10] used a mixed site-uncorrelated bond-correlated percolation process for fire (although not using this terminology).The site percolation is due to uniformly randomly distributed tree cells perfectly blocking fires and the bond percolation due to flammable cells (grass and savanna saplings in [10]) spreading fires with a given probability.The correlation in bond occupation probability occurs due to the infection dynamics, as explained in Section S2C.One can use the complement of the burning probability of this mixed percolation process, i.e. 1 − ⟨Q p b (p s )⟩, as survival probability instead of that used in [10] (see Fig. S11b for a plot of ⟨Q p b (p s )⟩ as a function of the infection probability between flammable cells p b and the probability of a cell being flammable p s ).If one does so, one can write the mean-field recruitment rate of savanna saplings during a small time interval as ω(p b , p s ) [ S ] , where ω is ω(p b , p s ) := max[ω 0 − ϕN p s ⟨Q p b (p s )⟩∆ω, 0], (S25) with ∆ω>0 the per fire decrease in recruitment rate due to burning and ϕ the fire ignition rate in flammable cells.Note that, according to [10], the total flammable area is p s = [ S ] + [ G ] (i.e., the area of grass and savanna saplings).The reasoning is that there are on average ϕN p s ignitions, each causing a fire that on average affects a proportion ⟨Q p b (p s )⟩ of flammable cells, thereby lowering the recruitment rate by an amount ∆ω.If [ S ] and [ T ] cells are uniformly randomly placed in the area affected by fire, then it follows that the recruitment rate is ω(p b , p s ) [ S ] .
For the approximate effect on forest trees [as included in 12], one needs to take into account that forest trees are assumed [in 12] to block fires perfectly.Therefore, they are not in the flammable part of the landscape, but instead share an interface with it.The mean-field rate of forest loss due to fire during a small time interval is then ζ(p b , p s ) [  with p s = [ S ] + [ G ] also.The reasoning here is as follows.There are on average ϕN p s ignitions, each causing a fire that affects a proportion ⟨Q p b (p s )⟩ of the landscape.Assuming that occurrences of burnt and forest cells are uncorrelated, one can write the interface between them as the number of forest-burnt pairs in a lattice: 4(⟨Q p b (p s )⟩p s ) [ F ] .For each such pair, there is a probability p f of spreading into forest, such that (when using the approximation of small p f as in the main text), the resulting forest loss is 4ϕN p 2 s p f ⟨Q p b (p s )⟩ [ F ] .For large domains, one may replace ⟨Q p b (p s )⟩ by the percolation probability P ∞ (p b , p s ) (shown in Fig. S11a), which is the probability that a flammable cell is part of the giant connected component.Percolation probability is defined as the probability that any grass cell belongs to the giant component [26].Susceptibility is defined here as in [56], using Q as order parameter.The dash-dotted blue line indicates the location of the infinite-size percolation threshold for uncorrelated mixed site-bond percolation (taken from [54]).For a domain size of 100x100 cells and at the shown resolution, the percolation threshold of mixed site-bond percolation matches that of the infinite size system (f).The GBA process' percolation threshold lies at higher values (c) due to spatial correlation of pg as explained in the text.The top row is more noisy than the bottom row because for standard mixed percolation, we were able to obtain the whole distribution of cluster sizes for a each realisation and computed their statistics using percolation theory [26], whereas for the GBA process, each simulation only resulted in one sample.The colour scale was taken from [55].

FIG. 1 .FIG. 2 .
FIG. 1.The FGBA stochastic cellular automaton: (a) state transition diagram (coloured rates: spread to neighbouring cell, black rates: spontaneous conversion within cell), (b) example time series of a simulation starting at zero tree cover, (c,d) 10 2 x and 10 7 x zoom of (b), (e-h) snapshots of a simulation at indicated times for low fire ignition rate (ϕN =0.075).The fire in (f) spreads throughout grassland in the whole domain whereas that in (g) went extinct locally because forest splits grassland in clusters (notice the area of ash near the top).Remaining parameters are shown in Table M1.Domain size: 200x200 cells.

FIG. 5 .
FIG. 5. Resilience of forest-grass landscape to perturbations.(a-c) Spatially resolved contributions to forest gain (green) and loss (red) rates at the forest perimeter of the largest grass cluster: (b) for the saddle solution (where ∆ gain F ≈∆ loss F ), (a) for a perturbation of the saddle with more forest at the perimeter (resulting in ∆ gain F >∆ loss F ), and (c) for a perturbation of the saddle with less forest at the perimeter (resulting in ∆ loss F >∆ gain F FIG. S3.Emergent relations between key quantities and forest area[F] compared to the mean field with site percolation (dots: simulations, lines: corresponding mean field quantities from Section S4B): (a-c) forest perimeter[FG] (green) and grasslandweighted forest perimeter ⟨[FG]⟩cg (orange), where green line equals 4[F][G] and orange line was called ⟨[FG]⟩ u cg in the main article, (d-f) forest gain terms and loss terms in Eqs. 5 and 6, (g-i) forest area rate of change (d/dt)[F] from Eq. 9. Columns correspond to vertical dashed lines in Fig. 2 (ϕN =0.257, ϕN =0.38, ϕN =1.32).Simulation results are identical to Fig. 4. Domain size: N =100x100 cells.See Section S4B for details of derivation for ⟨[FG]⟩ u cg and [FG] mf .
FIG. S5.GBA process with ϕ=λ=0 (equivalent to square lattice SIR spreading): ⟨Q⟩, χ versus bond occupation probability p b := pg (Eq. 3) and site occupation probability ps=[G] 0 .(a) Mean cumulative probability of being burnt ⟨Q⟩ (expectation of Eq.S4).(b) Susceptibility χ (Eq.S5).(c) Susceptibility χ compared to the mean-field percolation threshold p s,mf , given in Eq.S6 (dashed black).The dash-dotted blue line indicates the location of the infinite-size percolation threshold for uncorrelated mixed site-bond percolation (taken from[54]).The GBA model's percolation threshold lies at higher values (b) due to spatial correlation of pg as explained in the text.The colour scale was taken from[55].See Fig.S11for more detail.
FIG. S9.Comparison of steady state forest as a function of ignition rate in the time-separated mean-field and in simulations (dots with error bars: simulations, lines: mean field).See legend for details.Note, 'uniform random' refers to the placement of forest cells in the domain.
F ] , where ζ is ζ(p b , p s ) := 4ϕN p 2 s p f ⟨Q p b (p s )⟩,(S26) FIG.S11.GBA process with ϕ=λ=0 (equivalent to SIR spreading on the square lattice) in terms of bond occupation probability p b := pg (Eq. 3) and site occupation probability ps :=[G] 0 (a-c) and compared to standard mixed site-bond percolation (df).Shown quantities as a function of bond and site occupation probability (based on 1024 realisations for (a-c) and on 512 realisations for (d-f)): (a,d) percolation probability P∞, (b,e) mean proportion of burnt grass cells ⟨Q⟩ (see Eq. S4), and, (c,f) susceptibility χ=⟨Q 2 ⟩/⟨Q⟩.Percolation probability is defined as the probability that any grass cell belongs to the giant component[26].Susceptibility is defined here as in[56], using Q as order parameter.The dash-dotted blue line indicates the location of the infinite-size percolation threshold for uncorrelated mixed site-bond percolation (taken from[54]).For a domain size of 100x100 cells and at the shown resolution, the percolation threshold of mixed site-bond percolation matches that of the infinite size system (f).The GBA process' percolation threshold lies at higher values (c) due to spatial correlation of pg as explained in the text.The top row is more noisy than the bottom row because for standard mixed percolation, we were able to obtain the whole distribution of cluster sizes for a each realisation and computed their statistics using percolation theory[26], whereas for the GBA process, each simulation only resulted in one sample.The colour scale was taken from[55].
[ FG ]⟩ cg lie on a narrow band around some steady state functions[FG ] * and ⟨ [ FG ] ⟩ * cg of [ F ] (and ⟨

TABLE M1 .
Reactions and rates (y −1 ) in the CA. β