Synthetic soil crusts against green-desert transitions: a spatial model

Semiarid ecosystems are threatened by global warming due to longer dehydration times and increasing soil degradation. Mounting evidence indicates that, given the current trends, drylands are likely to expand and possibly experience catastrophic shifts from vegetated to desert states. Here, we explore a recent suggestion based on the concept of ecosystem terraformation, where a synthetic organism is used to counterbalance some of the nonlinear effects causing the presence of such tipping points. Using an explicit spatial model incorporating facilitation and considering a simplification of states found in semiarid ecosystems including vegetation, fertile and desert soil, we investigate how engineered microorganisms can shape the fate of these ecosystems. Specifically, two different, but complementary, terraformation strategies are proposed: Cooperation-based: C-terraformation; and Dispersion-based: D-terraformation. The first strategy involves the use of soil synthetic microorganisms to introduce cooperative loops (facilitation) with the vegetation. The second one involves the introduction of engineered microorganisms improving their dispersal capacity, thus facilitating the transition from desert to fertile soil. We show that small modifications enhancing cooperative loops can effectively modify the aridity level of the critical transition found at increasing soil degradation rates, also identifying a stronger protection against soil degradation by using the D-terraformation strategy. The same results are found in a mean-field model providing insights into the transitions and dynamics tied to these terraformation strategies. The potential consequences and extensions of these models are discussed.

Semiarid ecosystems are threatened by global warming due to longer dehydration times and increasing soil degradation. Mounting evidence indicates that, given the current trends, drylands are likely to expand and possibly experience catastrophic shifts from vegetated to desert states. Here, we explore a recent suggestion based on the concept of ecosystem terraformation, where a synthetic organism is used to counterbalance some of the nonlinear effects causing the presence of such tipping points. Using an explicit spatial model incorporating facilitation and considering a simplification of states found in semiarid ecosystems including vegetation, fertile and desert soil, we investigate how engineered microorganisms can shape the fate of these ecosystems. Specifically, two different, but complementary, terraformation strategies are proposed: Cooperation-based: C-terraformation; and Dispersion-based: D-terraformation. The first strategy involves the use of soil synthetic microorganisms to introduce cooperative loops (facilitation) with the vegetation. The second one involves the introduction of engineered microorganisms improving their dispersal capacity, thus facilitating the transition from desert to fertile soil. We show that small modifications enhancing cooperative loops can effectively modify the aridity level of the critical transition found at increasing soil degradation rates, also identifying a stronger protection against soil degradation by using the D-terraformation strategy. The same results are found in a mean-field model providing insights into the transitions and dynamics tied to these terraformation strategies. The potential consequences and extensions of these models are discussed.

Introduction
Global warming is changing the dynamics and resilience of ecosystems, damaging many of them and creating the conditions for widespread diversity loss [1][2][3]. Because of the presence of nonlinearities, many ecosystems could suffer the so-called tipping points [4][5][6], responsible for community collapses and for transitions towards more degraded states [7,8] (e.g. coral bleaching [9] or lake eutrophication [10]). Among these systems, drylands (which comprise arid, semi-arid and dry-subhumid ecosystems) are a specially fragile subset of major importance: they include more than 40% of terrestrial ecosystems and host a similar percentage of the current human population [11,12]. Increasing aridity is pushing these ecosystems towards serious declines in microbial diversity, land degradation and loss of multifunctionality as desert states are approached [11,13,14]. Dedicated efforts have been addressing several avenues to both understanding how transitions can be anticipated by means of warning signals [15][16][17] and even prevented [11,18,19]. Catastrophic transitions i.e. the transition that occurs at the tipping point, from vegetated to desert ecosystems can be really fast [19], even in large ecosystems [20]. This type of sudden transitions have been identified in the evolution of the vegetated Sahara to the current desert state, which occurred less than 6000 years ago [21][22][23].
Drylands are characterized by the presence of organisms that have adapted to low moisture availability, damaging UV radiation, and high temperatures [24][25][26]. A rich community structure and the maintenance of physical soil coherence are essential to prevent drylands from degradation [27][28][29]. In this context, some universal types of interactions that occur among species in arid habitats inevitably lead to breakpoints associated with the existence of multiple alternative states, identified in field data [14,30]. A well-known class of these interactions takes place among vascular plants and is known as facilitation [31][32][33]. Facilitation consists on non-trophic interactions between individuals mediated through changes in the abiotic environment favouring individual growth and reproduction [34][35][36]. An example of this process is the fertile islands in the savannah ecosystems, where the presence of trees increase the soil organic carbon thus increasing soil moisture and nutrients [33]. This process facilitates the implantation and growth of new vegetation (grasses and trees) [37,38].
A common outcome of spatially extended ecosystems is the emergence of spatial patchiness e.g. peatlands (see [39] and references therein) or the Kalahari desert [38], to cite a few. These spatial patterns have been extensively studied using computational models [39,40], for instance, models only considering competition [41]. Other models including both competition and facilitation processes also typically display self-organized patchiness [38,[42][43][44][45]. Such patterns are often remarkably organized in space [45]. Transitions between different ecological states, e.g. vegetated towards bare soil, can be sometimes observed in the same landscape under the presence of environmental gradients [31]. In this context, it has been conjectured that spatial correlations can be used as indicators of forthcoming green-desert shifts [15,46,47]. The presence of catastrophic transitions deeply modifies our perception of risks associated with climate change and land degradation. Once a tipping point is reached, unstoppable runaway processes are unleashed [8]. Despite much progress having been made in modelling drylands [47][48][49] as well as in identifying warning signals [15][16][17]46], feasible strategies to prevent green-desert transitions are currently scarce.
In a recent paper [42], it has been shown that, by tuning some particular features of models exhibiting catastrophic shifts, one can modify their nature and location in parameter spaces. Based on a theoretical approach, the authors suggested that the amount of stochasticity (both demographic and extrinsic) or population dispersal range can play that role. Indeed, the impact of plant dispersal ranges on the nature of transitions has been recently investigated in a simple metapopulation model with facilitation [43], showing that short-range dispersal might involve continuous transitions. Examples of gradual regime shifts in spatial patterns have been also described theoretically [44]. What type of microscopic mechanisms could actually avoid green-desert transitions? In this paper, we are extending previous models of dryland dynamics [19,47] to predict the impact of bioengineering strategies on these fragile ecosystems. In this context, several bioremediation approaches have been developed in the last two decades to enhance and stabilize soil biocrusts (for an overview, see [26] and references therein). These include a diverse range of approaches addressed to improve moisture and soil texture, from enrichment of key nutrients [50] or enrichment by cyanobacteria [51][52][53], to large-scale straw chequerboard barriers used to stabilize sand dunes [54].
More recently, ecosystem terraformation has been suggested as a novel approach against green-desert shifts [55][56][57][58][59]. In a nutshell, some existing microorganisms in a degraded ecosystem, such as cyanobacteria of the soil crust, could be minimally modified by means of synthetic biology techniques to help improve soil moisture and create a cooperative feedback between vegetation or moss cover (see [60]) and the soil microbiome [61]. This engineering scenario aims at building synthetic soil composed of synthetic strains mixed with natural microorganisms embedded, for example, in the soil crust. By doing so, potential tipping points could be made much more difficult to be achieved [58]. What might be the impact of these strategies on the large-scale, long-term dynamics of semiarid ecosystems? How can new or enhanced ecological interactions (cooperation and microbes dispersal capacities) modify the presence of tipping points? As shown below by means of both mathematical and computational models, the engineering of microorganisms present in arid and semiarid soils may easily expand the potential of ecosystems' persistence.

Modelling terraformation of semiarid ecosystems
In order to predict the impact of synthetic bioengineering strategies in arid and semiarid ecosystems, a spatial model given by a stochastic cellular automaton (CA) is employed. Let us first consider the microscopic rules associated with the dynamics, as described by a set of probabilistic transitions among the three defined ecological states. After this general description of the transitions, we will investigate a mean-field model without stochasticity and the stochastic CA. The mean-field model will allow us to identify the nature of the tipping points and the expected equilibria (which may be useful for low stochasticity levels i.e. large system's size).
The basic transition diagram is displayed in figure 1c, which includes the three potential states characterizing a given patch, namely These states, as mentioned, are defined by desert (D) patches, by fertile soil (S) composed by patches occupied by engineered microorganism(s) with their natural community (e.g. soil crust), and by vegetation (V), respectively. This is, of course, an oversimplification that ignores most of the complexity and diversity involved, but allows for an analysis of the dynamics considering key ecological interactions such as transitions between states and both competition and facilitation processes. The model described here is derived from the one proposed by Kéfi et al. [47], now considering that unoccupied sites (non-degraded and non-vegetated) contain soil crust plus (engineered) microorganisms. As discussed in a previous work [58], we aim at describing how an appropriate synthetic modification of microorganisms could help to maintain the stability and resilience of semiarid ecosystems.
The first rule, which considers the transition from degraded to fertile soil containing synthetic strains, will take place with a probability Here, r is the probability of spontaneous recovery of fertile soil due to external factors such as increased humidity, accumulation of organic matter, etc. Since we are interested in identifying the impact of the engineered strains in the recovery of the desert soil to the fertile one we will set r = 0. This will allow us to provide clearer results concerning soil recovery derived from terraformation. Constant f denotes the facilitation due to the influence of the surrounding vegetation. The variable ρ v is the local density of vegetation (number of vegetated neighbours over the total number of neighbours). Finally, Γ(ρ s ) denotes the spreading capacity of the microorganisms associated with their synthetic engineered properties (see below).
The second transition considers the reverse situation, namely a fertile soil-to-desert transition where ɛ is the probability of loss of fertile soil due to increased aridity or to other soil degradation processes. In this paper, we consider (see below) a modification of the degradation rate mediated by a vegetation-dependent function Y(r v ). The third transition rule involves the colonization of available S patches by vegetation This term has the same form as the colonization used in [47]. The probability δ balances the influence of the local and global vegetation to produce the germination of new plants in a given site. Seeds are produced and germinate with probability b. However, they can also be degraded during the royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 7: 200161 dispersion (with probability c). Finally, an exponential decay of vegetated patches to fertile soil occurs following the transition rule where m is the death probability of plants, which can occur due to pathogens, grazing, or to the recollection by livestock activity.
Here two different, but complementary, terraformation strategies are considered and weighted by means of constants α and β. Their role within the transition diagram is displayed in figure 1c. They involve: 1. C-terraformation: Engineering cooperative loops between vegetation and soil organisms. Here, a decreasing function of the vegetation cover, namely introduces a reduction in the impact of desertification due to increased soil quality (by e.g. nutrient deposition) favoured by the synthetic microbial population. The efficiency of this term is weighted by the constant α: large values of α imply a lower soil degradation rate due to the action of the synthetic microbes which are cooperatively coupled with vegetation at a local scale. 2. D-terraformation: Engineering the capacity of microbial spreading (e.g. faster replicative rates of the microbes and/or increased formation of endospores able to successfully colonize local surroundings). The enlarged view shows the patchy structure of vegetation, surrounded by bare soil and soil crust (shown in (b)). Our computational model considers three well-differentiated ecological states: vegetation (V), functional soil (including soil crust, S) and bare soil (desert sites, D). The local dynamics follows a set of stochastic rules that allow transitions between the three different states, as indicated in (c) using black arrows. Moreover, the presence of vegetation cover has an impact on these transitions, as indicated by the links ending in black circles (e.g. vegetation provides organic matter to the soil which increases water retention and the microorganisms in there degrade this matter for vegetation uptake). These are the basic rules used in previous literature modelling semiarid ecosystems [19,46,47]. The proposed scenario of synthetic terraformation considered in this article is indicated using coloured links. The orange line denotes C-terraformation (a synthetic strain establishing cooperative feedbacks with vegetation), inhibiting the transition from S to D. The blue line indicates D-terraformation, allowing for increased dispersal of the engineered organisms (SYN) thus favouring the transition form D to S. (d) We use a stochastic cellular automaton including these patch types, here displayed in the lattice projected on a sphere. royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 7: 200161 The dispersal is here considered proportional to the term Here, ρ s is the local density of fertile soil. In other words, with probability β, engineered strains improve their dispersal capabilities, thus changing the properties of desert patches by retaining humidity, depositing organic matter, etc.
Setting α and β to zero, the CA model recovers the original model introduced by Kéfi et al. [47].

Models and results
In the next sections, we will explore the dynamics tied to the rules described above using both a well-mixed (mean-field) and a discrete spatial setting using a stochastic Cellular Automaton (CA). The first approach, developed in §3.1, employs differential equations to predict the presence of tipping points and changes in their location in the parameter space due to the effects of the addition of a synthetic strain. The second approach, studied in §3.2, explicitly deals with a spatially extended population where local interactions take place on a two-dimensional lattice and the set of rules are applied in a probabilistic manner, thus considering stochastic effects. Finally, we provide insights into the nature of the tipping points found in the spatial model, focusing on the early warning signals for the system without and with the synthetic engineered microorganisms.

Mean-field model
The differential equations model, which includes the interactions and processes tied to the bioengineered synthetic strains, is given by Equation (3.1) contains a logistic term that includes variable S as a multiplicative term, indicating that plants need viable soil to persist, and an exponential decay proportional to m. Equation (3.2) includes the positive effects triggered by vegetation and the existing soil cover, as well as negative terms that are in fact symmetric to those present in the previous equation for V. Here recall that 3) provides the dynamics of desert patches resulting from desertification, recruitment to viable soil patch, and an exponential decay from S to D introduced with Γ = β S.
Since the three classes of patches cover the entire lattice in the spatial model, it is possible to normalize and consider the states as population fractions in the mean field approach, i.e. V(t) + S(t) + D(t) = 1. This actually assumes a constant amount of available sites that can transition between them, also allowing to reduce the three-variables system to a two-variables one using the linear relation which needs to satisfy dV/dt + dS/dt + dD/dt = 0 (recall that r = 0) since the sum of all states remains constant through time. The reduced system is thus given by Note that now the fraction D is automatically obtained from equation (3.4) once the fractions of states V and S are determined. We must note that equations (3.5) and (3.6) have been recently studied in [19] for the non-engineered system (that is, taking α = β = 0).
We are specially interested in those scenarios involving a full dominance of the desert state, focusing in the transitions between states. We have identified four equilibria for equations (3.5) and (3.6), two of royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 7: 200161 them implying the extinction of vegetation and another one allowing for the stable coexistence of the three ecological states. The first equilibrium is at the origin, i.e. P Ã 1 ¼ (V ¼ 0, S ¼ 0), were the system becomes a desert (D = 1). The second equilibrium, keeps the fertile soil and the desert patches without vegetation, and is given by . This state is interesting from the point of view of the terraformation strategies, since its existence depends on the chosen engineering strategy. If β = 0, when no engineered organisms are found in the soil crust, this equilibrium is not biologically meaningful (S(t → ∞) = −∞).
Two more equilibrium points, labelled P Ã 3 and P Ã 4 , have been identified numerically (see electronic supplementary material, figures S1-S4 and S6-S8). The equilibrium P Ã 4 (numerical results suggest it is always stable within the simplex; see below) involves V(t → ∞) > 0, S(t → ∞) > 0 and V(t → ∞) + S(t → ∞) < 1, allowing the coexistence of the three ecological states when it is an interior equilibrium. P Ã 3 is a saddle point. See below and electronic supplementary material, figures S1-S4 for further details on the dynamics of the equilibria on the simplex (V, S).
The conditions defining the local stability of each equilibrium can be obtained by means of the eigenvalues of the Jacobian matrix J , given, for equations (3.5) and (3.6), by The eigenvalues of matrix J evaluated at an equilibrium provide its local (and linear) stability properties. The stability for equilibria P Ã 1,2 can be easily computed, while the stability of equilibria P Ã 3,4 will be characterized numerically. In the case of the desert equilibrium P Ã 1 ¼ (V ¼ 0, S ¼ 0), the eigenvalues are λ 1 = −m and λ 2 = β − ɛ. Since all the parameters are positive, this equilibrium will be stable when ɛ > β. For P Ã 2 ¼ (0, 1 À 1=b), its eigenvalues are λ 1 = ɛ − β and λ 2 = −m + b(1 − ɛ/β); see electronic supplementary material, figure S5 for a stability diagram for P Ã 2 in the parameter space (ɛ, β). Note that P Ã 1 and P Ã 2 suffer a transcritical bifurcation at β = ɛ, since the two conditions for this bifurcation are fulfilled: the two fixed points collide at the bifurcation value (at ɛ = β, P Ã 1 ¼ P Ã 2 ¼ (0, 0)), and they interchange the stability. The dynamics tied to increases of soil degradation rate are summarized in figure 2c (see also electronic supplementary material, figures S1-S4), where different phase portraits are displayed for the values of ɛ indicated in figure 2b panel (b.c) (from (i) to (iv)). The transition from the coexistence of vegetation and fertile soil to no vegetation is catastrophic, and is due to a saddle-node bifurcation between fixed points P Ã 3 and P Ã 4 (transition between phase portraits (ii) and (iii), see also electronic supplementary material, figures S1-S4). Once this bifurcation takes place, the only stable point is P Ã 2 (with fertile soil and desert areas). Then, further increase of ɛ involves the transcritical bifurcation previously mentioned, after which the origin becomes globally stable and the desert becomes the only possible state (transitions between phase portraits (iii) and (iv) of figure 2c), see also electronic supplementary material, figure S7 where the transcritical bifurcation is shown.
Finally, electronic supplementary material, figure S6 provides one-dimensional bifurcation diagrams tied to the increase of ɛ for the non-terraformed system and the three possible terraformation strategies royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 7: 200161 (C-, Dand (C, D)-terraformation). Here, we want to emphasize the presence of so-called delayed transitions (also called ghosts), which arise just after a saddle-node bifurcation takes place [62]. This is a dynamical phenomenon that involves extremely long transients once the bifurcations has occurred, and the time trajectories experience a long bottleneck before rapidly achieving, in equations (3.5) and (3.6), another attractor (the full desert state in electronic supplementary material, figure S6(a) and (b)). These long transients [63] are typically found in systems with strong feedbacks, such as cooperation, catalytic processes [64,65] or metapopulations with facilitation [43]. Also, this dynamical delay tied to saddle-node bifurcations has been recently described in both deterministic and stochastic well-mixed approaches for the non-terraformed system explored in this article (see [19] for more details).

Spatial stochastic model
The mean-field model studied in the previous section provides a first approximation to understand the qualitative dynamics arising from the nonlinear interactions of the studied system and the resulting tipping points, especially when including the proposed terraformation strategies. However, in order to test the robustness of these results in a more ecologically realistic setting, one has to take into account both local spatial correlations and stochastic effects [66]. To do so, we build a spatially explicit simulation model given by a stochastic cellular automation (CA). The CA incorporates the previously studied states (V, S, D) and transitions among them, which are now probabilistic. Spatial degrees of freedom are introduced by using a L Â L square lattice with periodic boundary conditions. At each time step, the following four transition probabilities defining the rate of change for each site k of the lattice (figure 1) are applied in an asynchronous manner: Here, M k is the neighbourhood of site k and M N ¼ 8 the number of neighbour (using a Moore neighbourhood). The nature of local (spatial) interactions is known to largely influence ecological dynamics [40,43,67]. This is particularly important in drylands, where carbon and water limitation deeply constrains the outcome of nonlinear exchanges [14,59], leading to spatial patterning [48,68]. Also, facilitation processes (involving strong nonlinearities) are known to introduce important changes in spatial systems, as compared with well-mixed ones. In this sense, recent research has found a shift from catastrophic tipping points to continuous ones, due to local spatial processes [42,43]. Similar analyses to those presented in figure 2 for the mean-field model are displayed in figure 3 for the CA simulations. For the sake of comparison, we have kept the same parameter values of the meanfield model, implemented as probabilities in the CA, also using as initial conditions a full vegetated lattice. The critical surface separating the parameter scenarios allowing the persistence of vegetation is displayed in figure 3a. Below the surface, V persists, while above the surface it is only possible to find states S and D. Here, similarly to the mean-field model, the yellow region indicates those values of α and β for which no critical degradation rate of the fertile soil is found (i.e. ɛ c = 1). Note that the surface for the spatial system differs from the one obtained with the mean-field model. This effect, taking into account that the lattice has 4 × 10 4 sites (large system size), is probably introduced by space more than by stochasticity. This change of the surface is especially visible for parameter α. This may be explained because soil degradation rate (ɛ) makes the vegetation decrease globally, but due to the local interactions of facilitation [33], plants remain in the ecosystem for larger rates of soil degradation.
The fraction of the three states computed at increasing ɛ is displayed in figure 3b. All diagrams (b.a-b.d) display the same phenomenon: the introduction of spatial correlations shifts the extinction of the vegetation to larger values of ɛ. Consistently with mean-field results, both the Cand D-terraformation displace the value of ɛ c to larger values, the (C, D)-terraformation strategy being much more efficient in doing so under the selected parameters. The spheres displayed in figure 3c show the spatial patterns at equilibrium for the regions labelled in panel (c) of figure 3b, with: (i-ii) coexistence of the three states; (iii) only fertile soil and desert patches; and (iv) the full desert scenario.
In order to identify the nature of the transitions tied to the increase of ɛ, we have computed the stationary values of the fraction of areas with vegetation and fertile soil. We recall that both mean field and well-mixed stochastic simulations for the system without terraformation revealed a catastrophic transition and delayed transitions due to a ghost [19]. Also, the mean-field model studied in §3.1 has revealed the presence of saddle-node bifurcations responsible for vegetation extinctions. Figure 4 shows these results together with the spatial patterns on the square lattice and time series. Specifically, figure 4a shows results increasing ɛ for the system without terraformation. For this case ɛ c ≈ 0.24. The same results are displayed by setting α = β = 1, for which the critical degradation moves to ɛ c ≈ 0.85. For both cases, the discontinuity of the transition is not so evident as in the mean-field model, even looking like a continuous one.
To identify the transition governing vegetation extinction, we have used an indirect method analysing the properties of the dynamics close to the transition value. The results displayed in electronic supplementary material, figure S9 indicate the presence of bottlenecking phenomena tied to the royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 7: 200161 extinction of the vegetation, meaning that the transition is governed by a saddle-node bifurcation. We must notice that we have used a lattice of side size L ¼ 50 for these analyses due to the huge computational cost of computing extinction transients in large spatial systems. Electronic supplementary material, figure S9(a) shows the transition for the non-engineered system. Specifically, panel (a.1) displays five realizations obtained with ɛ = 0.26, and the bottleneck region can be clearly seen (see electronic supplementary material, figure S6 for comparison with the mean field bottlenecks). The delaying effect of the ghost can also be seen in panels (a.2), (a.3) and (a.4), which indicate the values at which vegetation spends longer iterating before achieving the full desert state. The same results are shown for the system with (C, D)-terraformation. Here, the extinction dynamics of the vegetation (investigated setting ɛ = 0.86) is also given by a clear bottleneck. Previous research in catalytic hypercycles also revealed extinctions due to saddle-node bifurcations and the corresponding bottlenecking phenomena in a stochastic CA model [65].
As discussed above, the approach to tipping points is often associated with a qualitative change in the nature of the fluctuations exhibited by the system. How will our model behave when dealing with an extra phase associated with vegetation-free, synthetic soil crust? A divergence in a stochastic dynamical system can be detected by using standard deviation measures. We have used several measures that can detect the approach to criticality in both the non-manipulated and the terraformed scenarios (spatial variance and distribution of patch sizes). The spatial variance computed as standard royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 7: 200161 deviation σ ψ , with ψ referring to either V or S, is determined from where we have used the mean over the entire lattice, given by  The fluctuations diverge in a characteristic fashion, as figure 5a,b shows: whereas 〈σ V 〉 slowly grows and accelerates close to criticality, the same measure for 〈σ S 〉 displays a slight decline only growing very quickly close to ɛ c . Despite these trends, both variables exhibit high fluctuations close to the transition to the desert state.
An additional confirmation of the criticality associated with the transition to the desert state is obtained by determining the distribution of vegetation cluster sizes. In a nutshell, criticality is known to be linked to power law distributions of connected clusters. One particular instance of these scaling behaviours was found in [47]. Clusters of size S are composed by a set of S connected sites (i.e. sharing the same state and all elements in contact with another as nearest neighbours).
The number of clusters of a given size and their distributions have been obtained using a burning algorithm [69]. If P(S) indicates the probability of finding a cluster of size S, most measured patch size distributions found in arid and semiarid habitats follow a general form described as a truncated power law, namely P(S) S Àt exp (ÀS=S c ). This defines a power law decay with a scaling exponent τ that characterizes the statistical organization of systems close to percolation [69] along with an upper bound introduced by a cut-off S c in the exponential decay term. The cut-off defines the crossover from a regime of critical clusters ( power-law distributed) to that of non-critical clusters. As we get closer to percolation points, the system also displays divergent correlation lengths and S c ! 1 thus rapidly increasing the cut-off value and moving closer to pure power laws. Since the lattice is finite, we will always observe a rapid decay at the end of the distribution. These distributions can be systematically estimated from vegetation data [38,70,71]. Close to critical transitions we should expect a fat-tailed decay in P(S) following the power law form P(S) S Àt (when no characteristic scale is present, i.e. for S c ! 1). Instead, far below the transition point, the exponential term dominates and a characteristic scale is observed. In order to improve the smoothness of the distribution, the so-called cumulative distribution P . (S) will be used, and is given by where S 0 and S m indicate the smallest (the one-site scale) and largest (lattice) sizes. At criticality, the cumulative form scales as a power law P . (S) S Àtþ1 : (3:10) In our model, the shapes of the distributions change at increasing values of ɛ: from single-scaled to scalefree, as shown in figure 5c (original model) and figure 5d (terraformed system). In both cases, the same pattern is found, but the power-law behaviour is observed at much higher levels of soil degradation. Interestingly, the scaling exponent for both cases is the same: τ ≈ 1.75, in agreement with the results reported in Kéfi et al. [46]. The consideration of the extra processes tied to terraformation do not alter royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 7: 200161 a universal behaviour pattern, thus suggesting that, despite the parameter shift in ɛ c , the warning signals associated with both scenarios will behave in a similar way.

Discussion
The impact of anthropogenic-driven processes is rapidly increasing, leading to a degradation of extant ecosystems and, in many cases, pushing them closer to viability thresholds [1,2,72,73]. In the case of semiarid ecosystems, this degradation is produced mostly by global warming and intensive grazing [22]. The increasing pressure could end up in catastrophic transitions that are essentially irreversible. It is expected that the expansion of arid areas will increase in the next decades, with an associated increase in the likelihood of green-desert tipping points. Under this potential scenario, it is important to develop strategies of intervention aimed at avoiding these transitions. Theoretical results on the nature of these shifts and the parameter values at which they take place suggest that some generic factors (such as noise and dispersal [42,43]) could be crucial. Can realistic interventions help preventing green-desert transitions? As discussed above, several bioremediation strategies to face this problem have been proposed. Some recent proposals suggest how to exploit nonlinear features after tipping points based on periodic replanting [19]. More recently, an engineering approach based on synthetic biology has also been suggested [56][57][58]. In this article, we investigate this latest approach for arid and semiarid ecosystems employing both mathematical and computational models.
Beyond the standard vegetation-desert model, in our study we pay attention to the possibility that an engineered soil crust (a 'synthetic' soil) might have on the fate of fragile or endangered ecosystems. The original motivation is to describe, using a toy model, the impact of two different but complementary strategies: the so-called C-terraformation using a designed cooperative loop ( promoted by designed microbial strains); and the D-terraformation, consisting on engineered strains with increased dispersal capacities. Cooperative loops (included with the C-terraformation at the soil crust level in our modelling approach) have indeed been identified in natural ecosystems. For instance, termites activity may modify the nature of transitions to desert states by increasing water infiltration rates and soil moisture [74]. In this manuscript, we have explored the impact of these two synthetic terraformation strategies on the location and behaviour of tipping points.
We show that both terraformation strategies increase the resistance to soil degradation, pushing the tipping points to higher (or even much higher) critical degradation rates (ɛ c ). This is consistently predicted by both mean-field and spatial models. However, when space is explicitly included, the parameter space supporting this protection against catastrophic shifts is much larger, with a fluctuation dynamics (and underlying warning signals) that behave similarly between the nonengineered and the engineered systems, i.e. they show the same universality pattern.
Our theoretical predictions are limited by the simplifying assumptions that define our model. In particular, the description of soils ignores their development, spatial organization or diversity, as well as the network structure (e.g. trophic relations) of ecosystems [75]. These webs are known to display their own critical thresholds, and thus it remains open how such diverse communities may impact the results reported here. However, despite all these limitations, it is important to highlight that similar previous models of vegetation dynamics in drylands have been very successful in providing deep insight into their natural counterparts [15]. Some key universal properties might pervade the success of these models. In our context, this is likely to be related to the general patterns found close to phase transitions. Additionally, the terraformation scheme described here is likely to be effective when the cooperative interaction helps create the proper ecosystem engineering effect. By engineering ecosystem engineers, endangered arid and semiarid habitats might get protected (at least temporarily) from sudden collapse.
Finally, we may emphasize that ecosystems microbiome is a complex ecosystem itself within the full ecosystem. The introduction of a synthetic species could induce secondary effects at this level that could scale up towards higher ecological (e.g. trophic) levels. However, the synthetic biology approaches suggested here and in previous research [56][57][58] may re-use already existing microorganisms, which could be synthetically modified. The impact of a truly new species in the ecosystem may be investigated elsewhere by e.g. invasion (adaptive dynamics) or complex networks theory.
Data accessibility. The programming code for the cellular automaton model is provided as a electronic supplementary material.