Population dynamics of synthetic terraformation motifs

Ecosystems are complex systems, currently experiencing several threats associated with global warming, intensive exploitation and human-driven habitat degradation. Because of a general presence of multiple stable states, including states involving population extinction, and due to the intrinsic nonlinearities associated with feedback loops, collapse in ecosystems could occur in a catastrophic manner. It has been recently suggested that a potential path to prevent or modify the outcome of these transitions would involve designing synthetic organisms and synthetic ecological interactions that could push these endangered systems out of the critical boundaries. In this paper, we investigate the dynamics of the simplest mathematical models associated with four classes of ecological engineering designs, named Terraformation motifs (TMs). These TMs put in a nutshell different ecological strategies. In this context, some fundamental types of bifurcations pervade the systems’ dynamics. Mutualistic interactions can enhance persistence of the systems by means of saddle-node bifurcations. The models without cooperative interactions show that ecosystems achieve restoration through transcritical bifurcations. Thus, our analysis of the models allows us to define the stability conditions and parameter domains where these TMs must work.

Ecosystems are complex systems, currently experiencing several threats associated with global warming, intensive exploitation and human-driven habitat degradation. Because of a general presence of multiple stable states, including states involving population extinction, and due to the intrinsic nonlinearities associated with feedback loops, collapse in ecosystems could occur in a catastrophic manner. It has been recently suggested that a potential path to prevent or modify the outcome of these transitions would involve designing synthetic organisms and synthetic ecological interactions that could push these endangered systems out of the critical boundaries. In this paper, we investigate the dynamics of the simplest mathematical models associated with four classes of ecological engineering designs, named Terraformation motifs (TMs). These TMs put in a nutshell different ecological strategies. In this context, some fundamental types of bifurcations pervade the systems' dynamics. Mutualistic interactions can enhance persistence of the systems by means of saddle-node bifurcations. The models without cooperative interactions show that ecosystems achieve restoration through transcritical bifurcations. Thus, our analysis of the models allows us to define the stability conditions and parameter domains where these TMs must work. Here the different scales are shown at the left and potential mathematical or formal descriptions at the right. Several layers of complexity need to be considered (in principle) including (a) whole population dynamics of ecosystems, with the population levels of each species described in terms of N j variables. The time evolution of each of these variables would follow a deterministic or stochastic Lotka-Volterra formulation. A smaller scale level (b) considers the dynamics of a few given species, which can include synthetic candidates (here indicated as S) derived from a wild-type strain (here indicated as S w ) and a plant host H. This small subset can be described, as a first approximation, by means of a small number of coupled equations, defining a three-species subgraph. At the species, cell level, the mathematical description of molecular interactions typically involves many coupled equations with nonlinear responses among genes and signalling molecules (here we just indicate a typical form of these cooperative interactions terms Γ 1 , Γ 2 ). Finally, at the gene-sequence level, designed constructs (d) must be engineered in order to operate under predictable circumstances.
its new ecological function, but also the impact and the evolutionary potential. In this context, a rational design of either ecological interactions and habitat constraints should guarantee its efficiency limiting the impact and the evolvability of synthetic organisms. It was argued that some special classes of engineered ecological motifs might well provide such strategy [23,26]. We labelled these basic designs Terraformation motifs (hereafter TM).
Our approach requires dealing with multiple scales, as summarized in figure 1. Here we display four levels of complexity: the ecological networks and the flow of resources within them (figure 1a); the specific nature of interactions among pairs of species including both wild-type and synthetic strains with their hosts (H), figure 1b; the required regulatory components, considering the cellular networks that will operate within cells (and how to engineer them, figure 1c); as well as the bottom level description where genetic sequences and available genetic toolkits are considered (figure 1d). In this paper, we approach these TMs from the point of view of their underlying population dynamics (figure 1b). We will consider the minimal mathematical models associated with each TM and determine the conditions for the survival or extinction of the synthetic organisms.
The aim of this paper is to make a few initial steps towards a general ecological theory for TMs seeking to understand the way these ecosystems might behave and how their populations will achieve different equilibria. By general theory, we mean the development of theoretical models able to explain the expected behaviours for TMs under different initial conditions and parameter values. Following the scheme outlined in figure 1, two important levels of complexity are presented when studying whole communities (figure 1a), and also when we consider a detailed description of species as cellular networks (figure 1c). As occurs with standard population dynamics, we need to start with the simple, few-species models indicated in figure 1b, somewhat averaging the details defining each particular partner at the smaller scale and also ignoring the multispecies nature of interaction defining the upper scale. As will be discussed below, such mesoscopic approach makes sense in the contexts discussed here.
All of the motifs, which are displayed in electronic supplementary material, Fig. S1, are explicit instances of the four basic classes of TMs presented in Solé et al. [26]. Electronic supplementary material, Fig. S1, also shows possible model organisms that could be used as candidates for their bioengineering. Specifically, in this article, we explore cooperative motifs with: direct ( §2.2.1; electronic supplementary material, Fig. S1a) and indirect ( §2.2.4; electronic supplementary material, Fig. S1b) cooperation. We also investigate two different restoration systems: function-and-die motif ( §2.3.1; electronic supplementary material, Fig. S1c) and sewage and landfill motif ( §2.3.2; electronic supplementary material, Fig. S1d). Before studying these particular motifs, we first introduce mean field models for the evolutionary and ecological process that are common to the four TMs ( §2.1; electronic supplementary material, S1).

Population dynamics of wild-type and synthetic strains under mutation
Let us first consider the simplest scenario tied to the TMs. It does not define a TM in itself, but it does contribute to an important piece of the underlying population dynamics. Here we consider the presence of two populations: the wild-type (hereafter WT) and the synthetic (hereafter SYN) strains. We assume that the SYN has been obtained by engineering the WT strain and that the SYN can revert to the WT by mutation or by the loss of the introduced genetic construct during the process of replication or during the life of the SYN.

Wild-type strain-synthetic strain system with mutation
First, we introduce a system with a population of SYN (S) and WT (S w ) strains. Figure 2a displays schematically a manipulated organism (microbe) where a given gene has been modified, perhaps adding some additional synthetic constructs either inserted within the genome (such as S 1 ) or within a plasmid π (indicated as S 2 ). Either case, we simply assume that these constructs can be lost at a rate μ, giving place to the WT strain. The previous processes can be summarized using a transition diagram as shown in figure 2b, where the replication rates for the WT and the SYN strains are given by ρ w and ρ, respectively. If nothing else is included, the dynamics for these two populations can be described by the following couple of linear differential equations: The previous model is rather simple and the time solutions can be obtained explicitly (see electronic supplementary material, Section 1.1). A key process that needs to be considered is competition, which introduces a first nonlinearity constraining the persistence of the SYN strain. The dynamics with competition is explored in the following section.

Model with mutation and competition
Let us consider a situation where we simply consider the two species, the WT and the SYN, competing for available resources. In figure 3a, we summarize the structure of their interactions. They compete (as shown by the mutual negative feedback) and, following the model introduced in the previous section, also replicate at rates ρ and ρ w , respectively. The system with competition can be described by  Figure 2. (a) Transitions from a synthetic (SYN) towards a wild-type (WT) strain can occur, for example, once the introduced constructs (either within the genome as S 1 or as an external plasmid element π , indicated as S 2 ) are lost. Such a loss of genetic material will typically occur at some rate μ. The basic scheme associated with this process is shown in (b). S w and S are the population sizes of the WT and the SYN strains, with replication rates ρ w and ρ, respectively. (b) Transcritical bifurcation for three different values of μ at decreasing ρ, with ρ w = 1 involving the extinction of the SYN strain at ρ = ρ c (see electronic supplementary material, Section S1.2). (c) Phase diagram (ρ, μ) displaying persistence (grey) and extinction (white) of the SYN species. The dynamics in one-dimensional phase space S are illustrated as the parameters (from bottom-right to upper-left) approach and cross the transcritical bifurcation value (black and white circles denote stable and unstable fixed points, respectively).
Competition is introduced assuming a constant population (CP) constraint by means of the outflow term Φ(S) = Φ(S, S w ) (this term will be used in all of the following models). The CP assumption captures the competition of the WT and the SYN in the same ecological niche. As mentioned, when a given microbe replicates, daughter cells might lose the gene constructs introduced in the engineering process. Such scenario should be expected (and occurs often in experimental conditions) if the fitness advantage of the synthetic organism does not compensate for the metabolic burden associated with the maintenance of additional genetic information. Assuming a CP, the previous equations can be collapsed to one by considering S + S w = 1 andṠ + S w = 0, having dS dt which is formally equivalent to a SIS-like model of epidemic spreading [32]. The dynamics of this system is summarized in figure 3   . Terraformation motifs involving cooperation among synthetic (SYN) engineered micro-organisms and multicellular hosts (H). The SYN can revert to the WT either by mutation or by losing the engineered construct at a rate μ. Two motifs involving direct (a) and indirect (b) positive interactions among both partners defining a mutual dependency are explored. One potential scenario for this class is provided by dryland ecosystems, with plants being the hosts. In these habitats, the soil crust (c) provides a spatially, well-organized community of microbial species that can help engineering cooperative interactions. An engineered microbe capable of improving moisture retention can have a very strong effect on the underlying plant species. In soil crusts, a whole range of species adapted to water-poor conditions exist (drawing after Belnap & Lange [35]). Here we display: (1) mosses, (2,3) lichens, (4,5,7,9) cyanobacteria, has two fixed points: S * 0 = 0 (extinction of the SYN strain) and S * 1 = (1 − (μρ)/(ρ − ρ w )). This system, which behaves like the logistic equation with density-independent degradation, suffers a transcritical bifurcation, as shown in figure 3b as a function of parameter ρ, or by means of a two-parameter bifurcation diagram (figure 3c; see also the associated phase space diagrams). The mutation bifurcation value is given by μ c = 1 − (ρ w /ρ). Note that replication rates for the WT and SYN can also act as bifurcating parameters.
In the model above, we considered that the processes of mutation or gene loss were tied to the process of replication. However, the genetic construct could be modified or lost during the lifetime of the SYN strain, not necessary during cells' duplication. A model taking into account this process also undergoes a transcritical bifurcation similar to the model previously discussed (see electronic supplementary material, Section 1.2, for details on fixed points and bifurcation values for this model).

Mutualistic terraformation motifs
A candidate organism to be engineered for modifying ecological systems should not be capable of decoupling itself from other species in such a way that becomes an expanding invader [33,34]. One especially appealing scenario is given by engineered mutualistic interactions ( figure 4a,b).
In this context, cooperation has been successfully engineered experimentally in yeast [36] as well as in spatially extended populations of Escherichia coli [37]. Mutualism requires a double positive feedback where the synthetic species S benefits-and is benefited by-its host H. Ideally, design failure should end in the disappearance of the modified species. Because mutualism deals with two partners, the synthetic species will be constrained by the population of its mutualist partner and such a tight bond is especially convenient, as shown below.
Several targets can be conjectured. One particularly relevant class is given by the bacteria-root dependencies exhibited by plants [38] and particularly plant crops [39] with their surrounding microbiome. The main case study where this motif applies is provided by semi-arid ecosystems, already discussed in the Introduction section. In these ecosystems, a patchy vegetation cover is usually present, with species adapted to low moisture, extreme temperatures and high UV radiation.
A crucial component of these ecosystems is the biological soil crust (figure 4c-e), defining a complex living skin enclosed within a few centimetres of the topsoil [40,41]. Interactions between the soil crust and the vegetation have been modelled mathematically [42,43]. These are remarkable communities hosting a wide variety of species and largely mediating the energy and matter flows through the soil surface. They are known to help to preserve biodiversity and provide a reliable monitorization system for ecosystems' health. In general, the more arid the environment the less diverse is the community, and, since plants and the biocrust are strongly related to each other, increased aridity leads to a smaller vegetation cover, less organic carbon reaching the soil, decreased micro-organism diversity, reduced plant productivity and loss of multifunctionality [44][45][46].
The basic schemes representing the interactions between the different components of the motif are shown in figure 4a,b. These are of course oversimplified pictures, since we ignore the multispecies composition of the biocrust. This simplification is done with the goal of understanding the behaviour displayed by minimal models, in the spirit of fundamental population dynamics models [47][48][49]. The first case (figure 4a) involves a direct impact through some tight relationship with the host plant, which can be, for example, an engineered symbiosis [50]. The second case (figure 4b) relies on an indirect cooperation mediated by the influence of the S w species on, for example, moisture. Let us consider and analyse the two scenarios separately.

Direct cooperation
The first type of cooperation motif deals with a synthetic strain that enhances the replication rate of the target (host) species (figure 2a). Here the best example would be to start from a free-living species and engineer it in order to build a new strain that becomes an obligate mutualist. Such transition has been shown to be possible and has been created by artificially forcing a strong metabolic dependence [51][52][53][54]. These systems have shown that the final product can be a physically tight interaction between the two partners [37].
This case study can be approached by a system of coupled differential equations as follows: In this model, we make the assumption that the two strains (SYN with state variable S; and WT with state variable S w ) compete for available space and/or resources while the engineered strain is involved in a cooperative interaction with the host. Here Γ (H, S) is a growth function for the host (see below), K being its carrying capacity. Parameter is the density-independent death rate of the host. Constant η is the growth of strain S due to the mutualistic interaction with the host. The other parameters ρ, ρ w and μ have the same meaning as in the previous sections. The function Φ(S) also stands for the outflow of the system, introducing competition. As we previously did, under the CP assumption, we can collapse the dynamical equations for the microbial strains into one, here with Φ(S) = ηHS + ρS + ρ w S w , and thus the equation for the synthetic population reads now Here we consider Γ (H, S) = (r + γ S)H, assuming that the host is capable of growing (at a rate r) in the absence of the microbial strains, whereas the term γ HS stands for the cooperative interaction. Two

Strict cooperation
The first scenario considers strict cooperation, which requires that the host can only grow via cooperation (i.e. r = 0). For this particular case, the system has four biologically meaningful equilibrium points, given by ) and the pair P * ± = (H * ± , S * ± ) (see electronic supplementary material, Section S2). The fixed point P * 1 will be outside the positive (biologically meaningful) phase space when μ > ρ − ρ w . Under the condition μ = ρ − ρ w , the fixed points P * 0 and P * 1 will collide since P * 1 | μ=ρ−ρ w = P * 0 = (0, 0). As we will show below, such condition involves a transcritical bifurcation between equilibria P * 0 and P * 1 . Let us now study the local stability of the fixed points P * 0 and P * 1 by means of linear stability analysis. Since the expression of the fixed points P * ± is cumbersome, the analytic derivation of the eigenvalues for these fixed points is rather difficult, and their stability character will be determined numerically by means of phase portraits representation (all of the numerical results presented in this article are obtained by means of the fourth-order Runge-Kutta method with a time step δt = 0.1).
The stability of the fixed point P * 0 is computed from the characteristic equation |J(P * 0 ) − λ (0) I| = 0, J being the Jacobian matrix of the system and I the identity matrix. The eigenvalues of this fixed point are given by λ The first eigenvalue is always negative, and thus the stability of this fixed point is given by λ 2 . Hence, the fixed point P * 0 will be stable when μ > ρ − ρ w , a scenario under which the host and the SYN strain will become extinct.
Let us now characterize the stability of the second fixed point, P * 1 . This equilibrium point, if stable, involves the extinction of the host and the survival of the SYN strain. The eigenvalues computed from The equilibrium point P * 1 will be stable if The previous results on the different fixed points and their stability nature are displayed in figure 5. First, we display the dynamics in the parameter spaces (ρ, μ) (figure 5a) and (γ , μ) (figure 5b). Here, for each pair of parameters we solved the system numerically plotting those parameter combinations where H and SYN persist at equilibrium (grey region). Here, there exists a frontier separating the grey and white zones that is given by a saddle-node bifurcation, which creates the pair of fixed points P * + and P * − , which are a stable node and a saddle. These two equilibria are interior fixed points, and the stable node governs the survival of the host and the SYN strain.
The dashed line in the grey region of figure 5a separates two scenarios where the bifurcation between fixed points P * 0 and P * 1 takes place. In the grey region, the dynamics can be bistable, and the system can achieve persistence or extinction of the H and the SYN depending on the initial conditions and whenever the origin is stable. Figure 5c,d shows two bifurcation diagrams by tuning μ and ρ. The phase portraits of figure 5 display all possible dynamical scenarios with H-SYN extinction (figure 5e), H-SYN coexistence under bistability (figure 5f,g) and the H-SYN persistence without bistability, since after the bifurcation the origin becomes unstable and the node P * + is a global attractor (figure 5h).

Facultative reproduction and cooperation
The cooperative system considering facultative reproduction of the host (r > 0) has five fixed points. Two of them are also given by P * 0 = (0, 0), and P * 1 = (0, S * 1 = 1 − μ/(ρ − ρ w )) and P * ± = (H * ± , S * ± ) (with r > 0; see electronic supplementary material, Section S2), as we found in the previous section. For this system a new fixed point is found, named P * 2 = (K(r − )/r, 0). This new fixed point, if stable, will involve the persistence of the H and the extinction of the SYN.
Linear stability analysis reveals that the eigenvalues for the fixed point P * 0 are λ Hence, this equilibrium point will be stable provided r < and μ > ρ − ρ w . The stability of   the fixed point P * 1 is given by the eigenvalues The first eigenvalue will be negative when while the second eigenvalue remains the same as for the strict cooperation model analysed above. Also, the stability of the fixed point P * 2 is determined from the sign of the eigenvalues, given by λ This fixed point will be stable provided r > and Notice that the stability of this equilibrium also depends on parameters ρ, ρ w , μ, K, and .
The stability of the fixed points P * ± is also characterized solving the equations numerically, as we did in the previous section. For the numerical study we will use (if not otherwise specified) a value of r = 0.5 > = 0.05. By doing so we ensure that the fixed point P * 0 (which, if stable, involves the extinction of both H and SYN) is unstable. Figure 6 summarizes all the dynamical outcomes of this system. First, we display the equilibrium states of the system in the parameter spaces (r, μ) (i) and (ρ, μ) (ii) computed numerically. The space (r, μ) contains three different phases. For those values of > r, the outcome of the system is the extinction of both the H and the SYN (the black region in figure 6a), since the fixed point P * 0 is stable. The transition to the scenario with H > 0 and SYN = 0 is governed by a transcritical bifurcation between the fixed points   P * 0 and P * 2 (see phase portraits in figure 6d,e). After this bifurcation, the fixed point P * 0 becomes unstable (white circle in the phase portraits of figure 6e-h) and P * 2 stable (black circle). Similar to the case for strict cooperation studied in the previous section, the boundary of the region where both H and SYN persist defines a saddle-node bifurcation responsible for the creation of the fixed points P * + (node) and P * − (saddle), involving bistability. Then, further increase of r makes the saddle point to collide with the fixed point P * 2 in another transcritical bifurcation. Such a collision involves the interchange of stability between points P * − and P * 2 . After the collision, the fixed point P * 2 becomes unstable, and the saddle leaves the positive (biologically meaningful) phase space. Here, the stable node becomes asymptotically globally stable, bistability being replaced by monostability. Finally, figure 6c displays a bifurcation diagram using r as a control parameter. Note that this bifurcation diagram corresponds to the values of parameter r displayed with the horizontal dashed line in figure 6a (setting μ = 0.6). Here one can follow the series of bifurcations discussed above.

Indirect cooperation
As previously discussed, one of the most obvious candidates to apply the approach taken here is provided by semi-arid ecosystems and other water-controlled habitats where soil water interacts with a diverse range of soil and community properties, including carbon assimilation, transpiration rates or biomass production [55]. The biological soil crust is composed by a network of mutualists [56] and provides the ecological context suitable for vascular plants (figure 5e). It strongly influences key ecosystem processes and its diverse composition offers multiple opportunities for engineering cooperative loops. Specifically, we aim at describing the impact of an engineered strain capable of improving water retention in the biocrust. This is illustrated by the production of extracellular polysaccharides by cyanobacteria [57][58][59], which have been shown to affect hydrological soil properties, as well as other important features such as soil carbon and maintenance of structural soil integrity [60]. These and other molecules (from vitamins to phytohormones) have been recognized to play a key role in helping plant growth and development [61].
The potentially beneficial role of increased extracellular molecules has been exploited in field experiments in desert habitats and sustainable agriculture under adverse ecological and edaphic conditions [62] through the direct addition of polysaccharides [63], cultivating outdoors combinations of cyanobacteria and plants [64], and massive inoculation of selected microbial strains [65]. Several positive results have been reported from these studies, suggesting a major role played by ecological interactions connecting molecular, cellular and population responses.
A minimal model that encapsulates the indirect cooperative interactions (see electronic supplementary material, Fig. S1b) involving water (state variable W) is provided by the following set of equations: The function Ξ (W, H) = βψWH is the growth rate of the host, which depends on the availability of water.
Here, the constant β is the growth rate of the H depending on the water, while ψ is the fraction of water used by the H to grow. The population of the host has a logistic growth restriction with carrying capacity K, and a density-independent death rate, . Concerning the dynamical equation for the SYN strain, constant η is the growth rate of the S tied to the cooperation with H, ρ is the replication rate of the S, ρ w is the replication rate of the WT and μ is the rate of gene loss or mutation. Now, the dilution flow is given by Φ(S) = (ηH + ρ)S + ρ w S w . Following Klausmeier [66], the last equation includes three terms in the RHS: (i) a constant water input, a (rate of precipitation); (ii) water loss; and (iii) a term of water consumption by the vegetation. The function of water loss reads L being the maximum evaporation rate and S c the rate of inhibition of the evaporation due to the presence of S. This function introduces a specific modulation by means of an inhibition function which includes both the term LW as well as a nonlinear decay associated with the presence of the synthetic population, which is capable of reducing water loss. As we did for the previous models, we can reduce the system by using the linear relation S w = 1 − S, now having For the sake of simplicity, we will useρ = ρ − ρ w . This system has five different fixed points, P * 1...5 , with and two more fixed points P * 4 and P * 5 (see electronic supplementary material, Section S3, for the values of these fixed points and the Jacobian matrix). The eigenvalues for the fixed point P * 1 are given by while the eigenvalues for P * 2 are  Relevant parameters that could be engineered are the replication efficiency of the SYN (ρ), the benefit the SYN obtains from the host (η), and the evaporation inhibition due to the action of the SYN microbia (S c ). Figure 7 displays two-parameter phase diagrams, where the different dynamical scenarios can be visualized. Synthetic organisms shall have a higher expression load due to the SYN construct. This load would make the SYN grow slower than the WT (ρ < ρ w ). In order to counterbalance this effect and make the SYN strain able to survive, the SYN can take advantage of the vegetation. If there is no symbiosis (η = 0) the SYN organism only survives providedρ > μ (figure 7A). For η > 0, a bistable region exists and it becomes larger as the strength of symbiosis increases. The SYN survives even for ρ = 0 if the symbiosis strength is 1 ( figure 7B). For higher η values the vegetation survives for large evaporation rates (L), even with a low replication rate (figure 7C).
For a semi-arid ecosystem with a given fraction of vegetated area, the increase in the evaporation rate (e.g. increase in temperature) involves the full desertification of the system (see the changes between the phase portraits displayed in figure 7a,b). This process and the associated changes in the topology of the phase space can be visualized in electronic supplementary material, video S1. If evaporation increases in the engineered ecosystem without changing the replication rate, there is a region of bistability, and a saddle-node bifurcation is achieved (see electronic supplementary material, video S2). Ifρ is high enough the fixed point P * 4 collides with P * 2 in a bifurcation (this collision can be seen in electronic supplementary material, video S3). Once the system is optimally engineered, the vegetation can survive even if ρ < ρ w and the evaporation is much higher (figure 7d). Ifρ < μ, the SYN will not be able to survive alone (figure 7c; electronic supplementary material, Fig. S2a), otherwise the engineered organism will survive even without vegetation (figure 8e; electronic supplementary material, Fig. S2f).
The engineering of the system can change the dynamics and the vegetation-desert transition from a transcritical bifurcation (involving P * 3 and P * 1 ) to a saddle-node bifurcation (where both P * 4 and P * This scheme applies to a diverse range of possible targets, such as marine plastic debris (a) where plastic is the resource, entering the system at a rate α and spontaneously degraded but also actively degraded by the microbial strains S and S w which appear to be supported by the substrate. The formal motif diagram is shown in (b). In (c), we display a simplified motif where the synthetic strain has not been derived from a wild-type variant.
collide). This process can be seen in the bifurcation diagrams of electronic supplementary material, Fig. S3. The saddle-node bifurcation involves the emergence of bistability. If ρ is high enough the saddle-node bifurcation takes place after the transcritical bifurcation between the points P * 3 and P * 1 . In this scenario, the re-vegetation will take place if the evaporation decreases (electronic supplementary material, Fig. S3e). The two stable states are given by: the coexistence of the H, S and W; and, depending onρ, the desert alone (figure 7d) or the coexistence of both SYN and W (electronic supplementary material, Fig. S2g).
Another relevant parameter to be modified by the engineering of the SYN strain is S c . This parameter is related to the effect in the evaporation of water depending on the amount of SYN. When the threshold S c is low, the amount of SYN population needed to have a positive effect in the environment is very low. This means that with a small fraction of the SYN the system could support a vegetated state. The SYN strain could be able to produce, for example, a large amount of polysaccharides capable of retaining water. The rate of production will be closely related to the inverse of the threshold. Depending on the difference between the replication efficiency (ρ) the region of bistability changes. Ifρ < μ the bistability region is broader (electronic supplementary material, Fig. S2A), and ifρ > μ (electronic supplementary material, Fig. S2B) the coexistence state is broader. However, the limit where the saddle-node bifurcation takes place does not change (electronic supplementary material, Fig. S2).

Restoration motifs
One of the outputs of the Anthropocene is the accumulation of waste. In order to face this problem, we propose and investigate two restoration motifs given by the: function-and-die and sewage and landfill motifs.

Function-and-die motifs
The TMs do not necessarily need to act within natural ecosystems. An alternative scenario, to be considered here and in the next section, is to take advantage of extensive wasteland habitats that have been created by humans and in some sense are already 'synthetic'. These synthetic habitats offer opportunities for using bioremediation (including removal of undesired molecules), but can also act as novel substrates that can host useful synthetic micro-organisms.
A synthetic strain could use this substrate as a physical surface allowing it to grow and perhaps disperse. Oceanic plastic debris is an example of this situation (figure 8a). Plastics have been entering oceans and concentrating into large garbage gyres [67] since the 1950s at a rapid and global scale, leading to the generation of the so-called plastisphere [68,69]. This widespread class of anthropogenic waste is made of non-natural macromolecular structures that were not present in nature and thus there were no biological mechanisms expected at that time to degrade it. However, mounting evidence indicates that some species of micro-organisms have been adapted to this special class of substrate, effectively degrading it (see [70] and references therein) or at least contributing to its fragmentation and decay.
The TM analysed here takes its name from an intrinsically relevant property that defines an ecological firewall. The key idea is illustrated by a specific example that has been developed with the goal of repairing (self-healing) concrete cracks, which are a major challenge for the maintenance of infrastructures. The alkaline environment makes it difficult for most species to thrive but some species can be used for this purpose. Since repair requires filling a given volume (something that living things can do by reproducing themselves) and do it by means of a suitable material (and bacteria can do that too) synthetic biology appears as a potentially useful approximation here [71].
A microbe can be designed to grow and replenish cracks with calcium carbonate along with a secreted macromolecule that merges into a strong material [72,73]. A major advantage of this problem is that anaerobic bacteria are not going to survive outside the crack and thus selection immediately acts once the task is finished: the function (repair) is done and afterwards the synthetic strain is unable to survive. The right combination of genetic design and ecological constraints create a powerful safeguard.
More generally, we consider here the potential conditions for survival of a synthetic strain living on a given substrate that enters the system and is degraded. The synthetic strain can just degrade the resource or can additionally perform some given functionality. Since removal of plastic debris might actually be part of the goal, it might be unnecessary to use existing species associated with this substrate. Instead, it could be more efficient to simply design or evolve a highly efficient species capable of attaching to the plastic surface, being also able to outcompete other present species.
The mathematical model associated with the function-and-die motif (see electronic supplementary material, Fig. S1c) is given by The state variables for this model are the resource (R) and both the SYN (S) and the WT (S w ) strains. Here α is the constant income rate of the resource (e.g. plastic), δ is the spontaneous degradation rate of the reource and σ is the elimination rate of the resource due to the action of the SYN species. Additionally, η is the growth rate of the SYN strain associated with the degradation of the resource (μ and ρ w have already been defined).
For this model, the outflow term is given by Φ(S) = ησ SR + ρ w S w . Assuming again a CP constraint, S + S w = 1, the system can be reduced to the following model: This system has three fixed points, given by the only-plastic system) and the pair P * 1,2 = (R * 1,2 , S * 1,2 ). The coordinates of the fixed point P * 1 = (R * 1 , S * 1 ) are given by with Ψ = −4α(δ + σ )ηη w + (αη + δη w + σ (ρ w + μ)) 2 . The coordinates of the fixed point P * 2 = (R * 2 , S * 2 ) read like the coordinates above but with a change of sign, with Numerical results obtained for this model suggested that the coordinate S * 1 > 1 within the range 0 ≤ η ≤ 1. For this case, the fixed point P * 1 is outside the simplex and it is not biologically meaningful (recall that the CP constraint assumes that S w + S = 1, and thus S * 1 cannot be higher than 1). Under this scenario, the dynamics in the interior of the simplex is governed by the fixed point P * 2 . In order to check whether the fixed point P * 1 lives inside the simplex, we performed a simple numerical test. We computed the value of the coordinate S * 1 for 10 10 combinations of random parameters (with uniform distribution) within the ranges: α ∈ [0, 50], and δ, σ , η, ρ w , μ ∈ [0, 1]. For all these combinations, we obtained values of S * 1 > 1. The Jacobian matrix for this model reads From det |J (P * 0 ) − λI| = 0, we compute the associated eigenvalues for P 0 , given by We note that λ 1 is always negative and thus the stability of P * 0 will entirely depend on λ 2 . The change of stability of this point can be computed from λ 2 = 0. The critical values of the parameters in λ 2 that involve a change of sign of this eigenvalue are Following the previous critical conditions, the fixed point P * 0 will be unstable (i.e. saddle-point with λ 2 > 0, meaning that the SYN strain will survive) when μ < μ c , η > η c , η w < η c w , α > α c or δ < δ c . For example, at μ = μ c , both fixed points P * 0 and P * 2 collide since At the bifurcation value, these fixed points also interchange stability. Hence, a transcritical bifurcation is found for this motif. The same behaviour is found at η = η c , ρ w = ρ c w , α = α c and δ = δ c . Some examples of the bifurcation diagrams associated with this model are shown in figure 9a,b. Here we represent the equilibrium populations S * (computed numerically) of the SYN strain against the efficiency parameter η. A continuous transition given by the transcritical bifurcation takes place for η = η c , when the SYN strain overcomes the competitive advantage of S w . Given the definition of η c , for a fixed input and degradation of the resource and mutation rate, the condition η > η c is achieved once the advantages of the engineered strain overcome the growth rate of the WT.
An interesting feature of this diagram is that, even for large values of the reversal parameter μ we obtain high population values provided that η is large enough. The changes in the phase space at increasing η are displayed in figure 9c-f. In figure 9c the fixed point P * 0 is globally asymptotically stable (here η < η c ). Once η > η c (figure 9d-f ), the fixed point P * 2 enters into the phase plane having exchanged the stability with P * 0 at η = η c via the transcritical bifurcation, thus becoming globally stable. As mentioned, the increase of η involves the motion of P * 2 towards higher population values, meaning that the SYN strain populations dominate over the WT ones. The results of the analysis reveal that, provided that the resource is not scarce, we just need a slight advantage of the engineered strain to make it successful and reaching a high population. Moreover, if the resource declines over time, S w will remain high.
The potential relevance of this scenario is illustrated by the observation that pathogenic strains of Vibrio sp. might be a major player in the marine plastisphere (particularly plastic microplastic) as revealed by sequencing methods [74]. If this is the original (WT) strain, an engineered strain with no toxin genes and improved attachment to plastic substrates could be designed to replace the WT. Moreover, one could also consider a rather orthogonal scenario where the plastic garbage constitutes an opportunity for synthetic communities to thrive. Since plastic debris is known to be used by a large number of species as a stable substrate, it can be argued that it can be seen as an artificial niche that provides new opportunities for developing complex communities.
It should be noted that plastic degradation by micro-organisms is not necessarily good news: degradation (or accelerated fragmentation) of large plastic items leads to a faster transfer towards smaller plastic sizes, particularly microplastic that can be easily transferred to food webs [75]. Is removal of the human-generated waste a necessary condition for these designed motifs? An alternative terraformation approach could be using synthetic species that attach to the substrate without actively degrading it. The synthetic micro-organism could carry some beneficial function, such as providing useful molecules enhancing the growth or establishment of other species, thus again acting as ecosystem engineers.
Mounting evidence indicates that a rich community of species adapted to these substrates has been developing over the years. Metagenomic analyses indicate an enrichment of genes associated with surface-associated lifestyles [76]. Within a surrounding environment that is oligotrophic and speciespoor, the plastic garbage defines a novel niche that has been fairly well colonized by a wide variety of species attached to the plastic substrate. In many cases, the resulting microbial community provides the scaffold for other species to thrive. Some proposals on using synthetic biology to address the problem of plastic garbage included a project aimed to facilitate the stable adhesion of plastic pieces with the goal of creating plastic islands (see [23] and references therein). In such a context, we could consider a TM where the colonization by a given species performing other functionalities could be designed, perhaps taking advantage of the niche as an opportunity to build a synthetic ecosystem.

Sewage and landfill motifs
The last TM analysed is connected to a major class of waste generated by farming as well as by urban and, especially, by mega-urban areas associated with domestic, municipal and industrial sources. Urban centres incorporate massive infrastructures associated with the treatment of waste as an end part of the city metabolism [77]. Sewage systems offer a specially interesting opportunity to apply our approach. They contain large amounts of organic matter, along with a wide repertoire of molecules of different origins, from drugs to toxic chemicals. Because of the potential damage caused by organic matter-rich   . (a, b) A physical container (grey boxes) allows one to regulate input and output flows for these systems. In (a) the graph shows the resource-consumer structure of this motif, where both strains are supported by the same resource, while they compete and are connected through gene loss. A simpler alternative (b) does not require engineering of extant species since it is a completely artefactual ecosystem and its preservation is not required. A typical scenario would be sewage-related infrastructures (c, image obtained from Wikimedia Commons) where a rich microbial community is known to exist.
waters (which can promote blooms of heterotrophic organisms leading to oxygen depletion in rivers) sewage treatment deals with a combination of organic particles along with diverse filters and a treatment of the resulting sludge from anaerobic micro-organisms [78]. Similarly, landfills have been widely used as a cheap solution of storing waste, despite the environmental consequences involving pollution on a local scale associated with leaching as well as contributing to global warming due to methane emissions. Some other problems are related to the treatment of heavy metals, organic pollutants (particularly aromatic hydrocarbons) and other hazardous waste. In order to address these environmental problems, strains of micro-organisms could be used to target these molecules. This approach is known as bioremediation and has been used with different degrees of success [79,80]. In recent years, it has been suggested that the use of genomic search, along with systems and synthetic biology approximations should be considered as a really effective approach to this problem [27,81].
The basic designs, summarized in figure 10a,b (see also electronic supplementary material, Fig. S1d), incorporate some kind of 'container' indicating the presence of a physical boundary. This represents, for example, the sewage system of an urban centre or the spatial domain defining the limits of a landfill. In this way, we can define inputs and outputs associated with the inflow of water, organic matter or chemicals on the one hand and the outflow carrying other classes of molecules as well as micro-organisms on the other.
Here too it might be less relevant to preserve the existing species of microbes, given the less relevant motivation of preserving WT strains, thus making unnecessary to engineering from S w (figure 10b). Of course, the real situation is more complex in terms of species and chemical diversity, and the single box indicating the resource R encapsulates a whole universe of chemical reactions. But we can also consider R as a very specific target for the SYN strain.
The mathematical model for the sewage and landfill motif (electronic supplementary material, Fig. S1d; and figure 10a) is given by the following system of nonlinear dynamical equations: Here, as in the previous model, α and δ R are, respectively, the rate of waste income and the spontaneous degradation rate of this resource. The degradation of the resource by the SYN and the WT strains is parametrized by σ S and σ w , respectively. Some fraction of the degraded resource can be invested for growth and reproduction for both the SYN and the WT strains. The constants η and η w parametrize these processes. Assuming the CP constraint, again with S + S w = 1 andṠ +Ṡ w = 0, the outflow term is given by Φ(S) = (η w σ w + ρ w − δ w ) + S(ησ − η w σ w +ρ −δ). Using the previous conditions, the three-variable system is reduced to a two-dimensional dynamical system describing the dynamics of the resource (R) and the SYN strain (S), according to Here we also plot the population equilibria along the black dashed lines from the panels tuning μ andσ . In all of the analysis, we set α = 1 and δ R = 0.05.
We note that here, for simplicity, we setδ = δ S − δ w ,ρ = ρ − ρ w ,σ = σ S − σ w andη = η − η w . The fixed points for equations (2.4)-(2.5) are given by and the pair of fixed points P * ± (see electronic supplementary material, Section S4, for their values). The Jacobian matrix for the system above reads where Θ = [η(σ + σ w ) + η wσ ]. The eigenvalues of the first fixed point, obtained from det |J (P * 1 ) − λI| = 0, are given by Note that λ 1 is always negative, and thus the stability of P * 1 entirely depends on λ 2 . From λ 2 , we can define the critical μ value: When μ > μ c , λ 2 < 0 and thus P 1 is stable. Under this stability condition the SYN strain will become extinct.
In order to focus on the most interesting parameters from the engineering point of view (i.e.σ andη), we will hereafter take into account that both S and S w strains reproduce at the same rates in the absence of R (i.e.ρ = 0), also assuming that both strains have the same death rates (δ w = δ S , i.e.δ = 0). Under these assumptions, the equations now read and (2.8) For the system above, the fixed point P * 1 remains the same, but now there exists a single interior fixed point, given by P * The eigenvalues of the first fixed point, obtained from det |J (P * 1 ) − λI| = 0, fixingρ =δ = 0, are given by As mentioned above, the stability of this fixed point will depend on λ 2 , and now the critical μ value involving a change in the stability of P * 1 is given by Also, note that all the rest of the model parameters apart from μ are in λ 2 . Hence, the bifurcation can be also achieved tuning these parameters. For the case of μ, it can be shown that meaning that P * 1 and P * 2 collide at this bifurcation point, also exchanging their stability (see figure 11; electronic supplementary material, Fig. S4). Hence, we can conclude that this system suffers a transcritical bifurcation.
The impact of parametersσ ,η and μ on the equilibrium concentrations of the resource and the SYN strain is displayed in figure 11a. We note that the resource in the parameter space (σ ,η, μ) is always present. However, the boundaries causing extinction for the SYN strain are clearly seen (figure 11a (right)). These transitions are given by the transcritical bifurcation previously discussed. Similar to the previous models, the increase of μ involves the extinction of S, all the population being formed by the WT strain.
As expected, the increase of both σ w and η w also causes the extinction of the SYN strain. This actually means that if the WT has a fitness advantage in terms of resource degradation or metabolization, this population will outcompete the SYN strain. Note that this effect takes place when bothσ andη decrease (i.e. σ w and η w increase, respectively), and the WT strain is able to degrade the resource faster or to better metabolize this resource. Under these conditions, and due to the competitive exclusion principle, the SYN strain will be outcompeted by the WT one. This process is accelerated at increasing parameter μ. The dynamics tied to different and decreasing values ofσ are represented in the phase portraits of figure 11a(i) and (iv).
Two specific cases can be considered here from the model given by equations (2.7)-(2.8) involving two different ecological scenarios that could be achieved by means of different engineering strategies. First, we will consider that both strains degrade the resource at the same rate (symmetric degradation) while their reproduction due to the consumption of such resource can be different. In the second scenario, we will consider that both SYN and WT strains degrade the resource at different rates but the metabolic efficiency of the resource is equal for both strains. Consider first the case where both the WT and the SYN strains degrade the resource at the same ratesσ = 0, with σ w = 0 and σ S = 0. However, their ability to metabolize the resource and use it for reproduction can be asymmetric (i.e.η = 0). The fixed points under these assumptions are given by P * 1 (which remains the same as described above) and Note that the transcritical bifurcation will take place when S * 2 = 0 and when λ 2 = 0. In particular, it can be shown that the SYN strain will survive when μ < μ c , with The diagrams in figure 11b indicate such a smooth transition involving the extinction of the SYN strain. Together with μ, the other parameters involved in the transcritical bifurcation are α,η, σ w and δ R . For instance, the bifurcation diagram tuningη in figure 11b displays the transition at the valueη =η c , withη As previously mentioned, the resource is always present, and its population is constant for this case since its equilibrium value does not depend on μ andη. Note that the coordinate R * 2 is the same as R * 1 , which depend on parameters α, δ R and σ w . This is the reason why the equilibrium of the resource does not change at the bifurcation. It is important to highlight that the concentration of the resources will decrease at increasing σ w . Hence, under symmetric degradation of the resources, the degradation efficiency of the WT will determine the equilibrium concentration of the resource. The corresponding bifurcation diagram of figure 11b shows the transitions of the SYN strain at μ = μ c (left) andη =η c (right).
As shown in figure 11c the transition takes place atσ =σ c , with σ c = μ(δ R + σ w ) αη w . Figure 11c shows that the equilibrium concentration of the resource depends on parameters μ andσ , although the concentration of the resource for the parameters analysed remains large. Increasing σ w or μ involves the extinction of the SYN strain.

Conclusion
As our population rapidly increases our ecological demands also grow, pushing ecosystems to their limits. We probably move, thereby, to an uncertain future, where both ecosystems and societies could face the threat of catastrophic responses [82,83]. Humans have been highly successful as engineers, particularly as large-scale ecosystem engineers [84,85] since our activities have deeply modified the flows of matter and energy. To counterbalance runaway effects, we might need to develop new forms of ecological engineering. In this context, technological solutions grounded in geoengineering [86] have been suggested over the last decades [87][88][89] to counterbalance the damaging effects of rising carbon dioxide levels. The bioengineering alternative [23] offers a rather different (but perhaps complementary) path, where the technological solution has the potential for self-replication and thus for scaling up to the ecosystem-level requirements. However, a proper understanding of the implications and likelihood of these methods asks for a new synthesis of ideas and fields, from molecular biology to community ecology. In this paper, we have mathematically explored four large classes of TMs that encapsulate different strategies. All of them involve engineering a given (microbial) species that is already present in the given ecological context. The proposed design goes beyond a standard bioremediation strategy, aiming at taking advantage of nonlinear responses, such as those displayed by (but not limited to) ecosystem engineers [90,91]. In this context, the proposed framework aims to extend the bioremediation picture [27] to larger scales and under a new set of ecological-grounded rules. Our models represent a first step in defining a population dynamics theory of synthetic ecosystems (considering the so-called TMs), that incorporate some classical modelling approaches with a constraint imposed by the nature of the synthetic component (which can revert to a WT strain). Our analysis provides a first glimpse of the potential relationships between key parameters allowing the SYN strain (and the implemented functionalities) to persist.
How reliable are our predictions? Several important factors have been left aside in our study, from diverse types of functional responses to stochasticity. However, as occurs with other nonlinear systems, we conjecture that some key results will not be strongly affected by the incorporation of additional nonlinearities or stochastic fluctuations. In this context, some fundamental types of bifurcations pervade the qualitative types of behaviours and parameter domain shapes reported above. The TMs including mutualistic interactions (studied in §2.2) revealed the presence of saddle-node bifurcations responsible for the extinction of both hosts and SYN strains. Interestingly, the replication rate of the SYN strain becomes a key parameter in allowing the persistence of the entire system. This phenomenon has been identified in all of the models including cooperation ( §2.2). Our results indicate that restoration motifs operate gradually under transcritical bifurcations.
Because of the inevitable simplification imposed by our low-dimensional models, it can be argued that many potential biases might arise from diversity-related factors and spatial effects. Community dynamics might limit or even prevent the spread of the engineered strain, but we also need to consider how the changes derived from the engineering can propagate through the system (e.g. by horizontal gene transference). However, indirect evidence from manipulation experiments suggests that inoculation of micro-organisms can successfully change the organization and functionality of a given ecosystem in predictable ways, particularly in relation to soil crust ecosystems [32,59,[92][93][94][95]. Future work should experimentally validate the predictions made here and further explore the limitations and potential extensions of our formalism.
Data accessibility. All the needed methods and results to replicate this article are provided as part of the paper and the electronic supplementary material.
Authors' contributions. R.V.S. conceived the original idea and wrote the manuscript with J.S. and B.V. in close collaboration with the rest of the authors. All authors built and analysed the mathematical models. All authors gave final approval for publication.
Competing interests. We declare we have no competing interests. Funding. This study was supported by a European Research Council Advanced Grant (SYNCOM), by the Botin Foundation, by Banco Santander through its Santander Universities Global Division, by grant FIS2015-67616-P, by the PR01018-EC-H2020-FET-Open MADONNA project and by the Santa Fe Institute. This work has also counted with the support of Secretaria d'Universitats i Recerca del Departament d'Economia i Coneixement de la Generalitat de Catalunya. J.S. has been also partially funded by a 'Ramón y Cajal' Fellowship (RYC-2017-22243) and by the CERCA Programme of the Generalitat de Catalunya. The research leading to these results has received funding from 'la Caixa' Foundation.