A Synthetic Microbial Daisyworld: Planetary Regulation on the Test Tube

1. Introduction

Our biosphere is the result of a long-term evolutionary experiment where living life forms and their environments have been interacting closely and on multiple scales through millions of years. The composition of our biotas has been changing, sometimes in dramatic ways, as shown by the fossil record of life [2,3]. This evolutionary process has been taking place on a planet that has also experienced profound changes. Some were caused (or driven) by astronomical phenomena, from deterministic orbital cycles to fateful asteroid impacts. However, major changes took place as a consequence of the entangled nature between climate and the biosphere [4]. The deep connections between environment and life leave a mark in the geological past. How has life influenced climate and viceversa? As early pointed out by Vladimir Vernadsky [5], the emergence of life fundamentally transformed the geosphere. One specially interesting observation made by geochemist James Lovelock in the 1970s was the realisation that our planet should have been driven by our Sun into higher temperature regimes.

We know from our two closest planetary neighbours, Mars and Venus, that steady changes in physical parameters can trigger runaway effects leading to a planet devoted of life [14]. In general, positive feedback loops can drive a planet to extreme steady states, from boiling temperatures to a snowball [4]. And yet, life on Earth has emerged and diversified, somehow dealing with the impacts of external drivers. To explain such stability, Lovelock suggested that an active coupling between life and physical systems make up the planet [6,7,8,9]. To make this point more explicit, Watson and Lovelock proposed in 1983 a simple model of planetary regulation known as the Daisy World Model (DWM) [22,25]. In a nutshell, the DWM considered an ideal planet where two kinds of agents, namely black and white daisies, had a distinct impact on planet albedo thus changing the local temperature in a way that could allow for climate stability over a wide range of solar luminosity values. Since its formulation, DWM has become a canonical model for Earth System Science [26] and has been used to explore a very diverse range of problems, from tipping points [11,27], to evolutionary dynamics [23,24] to climate bistability in exoplanets [1,15]. In contrast with other problems, the planetary scale makes rather unlikely to approach this coupling between environment and ecology in terms of controlled experimental conditions. Within synthetic biology, successful modelling and implementation of engineered ecosystems has been achieved, including cooperative consortia [16,17], predator-prey systems [18] or even multispecies assemblies [19,20,21]. Could such a synthetic counterpart be found for the DWM?

A general formulation of the Watson-Lovelock model is described by means of a set of coupled differential equations, namely (see [28] and references therein):

$$\frac{d{x}_k}{dt}=x_k\left({\beta }_k\left(T\right)\left[1-\sum _{j=1}^n{x}_j\right]-{\delta }_k\right)$$

(1)

which correspond to an ecological model of species interactions. The basic feedback loops involved are represented in Figure 1a for the standard two-species WLM. Increasing luminosity can trigger the growth of both kinds of vegetation, at rates indicated by ${\beta }_k\left(T\right)$ (with $k=1,2$). In the best known (and simpler) version, two species are considered and the following model was explored, for a "world" involving black and white daisies

$$\begin{array}{c}\hfill \frac{d{x}_w}{dt}=x_w\left({\beta }_w\left(T\right)[1-x_w-x_b]-{\delta }_w\right)\end{array}$$

(2)

$$\begin{array}{c}\hfill \frac{d{x}_b}{dt}=x_b\left({\beta }_b\left(T\right)[1-x_w-x_b]-{\delta }_b\right)\end{array}$$

(3)

where $[1-x_w-x_b]$ introduces the available bare soil that both species can occupy. In Figure 1a we sketch the basic model components using a spatial version (see [10,12,13]), with dark, white and gray sites indicating black and white daisies and bare soil, respectively (lower inset, Figure 1a). These equations remind us of a competition model, but the relevant complexity is captured within the ${\beta }_i\left(T\right)$ factors (with $i=b,w$), defined by means of a single-humped function involving an optimum. The details of this model can be found elsewhere [22,25] but the crucial intuition is that the growth rates of each kind of daisy is affected by temperature in opposite ways. At low temperature, black ones will thrive since they will warm their local environment and spread. As T increases, white daisies are favoured because they reduce warming. The area covered by each species contributes to the albedo and affect local temperatures in such a way that a self-adjustment occurs, leading to a stable regime of local temperatures, as shown in Figure 1b. Here the grey, increasing trend indicates how temperatures will rise in the absence of biological control. The underlying populations of each type are also shown in Figure 1c.

The DWM, despite its specific traits and potential challenges in translating it to the real biosphere, offers a valuable principle of living homeostasis that warrants consideration under broader assumptions. The concept of a living ecosystem that can respond to external forces and readjust itself to preserve diversity or some global property is inherently intriguing. It opens the door to creating mathematical and computational models representing such homeostatic ecologies, but an even more relevant prospect lies in constructing a real living ecosystem that echoes Lovelock’s vision. Synthetic biology emerges as an ideal candidate for achieving this ambitious goal. Through genetic engineering techniques, scientists have successfully designed cells with novel functionalities and orchestrated their interactions in intricate ways. Consequently, the main objective of this paper is to demonstrate the feasibility of constructing a synthetic daisy world..

Creating a living surrogate of the original DWM could help exploring the general problem in novel ways, as well as providing a rich context to explore the role of lower-scale features (such as molecular regulation) on the global regulation processes. The challenge is not minor: an experimental surrogate of the planetary coupling between life and environment involving a cooperative feedback is far from obvious. However, although a potential choice could involve using temperature as the driving parameter, there are other no less important properties that have been also controlled at the planetary level. One of them is acidity [8]: despite of the tendency towards acidification associated to an increasing oxidation of the atmosphere, the mean pH of the oceans has been remarkably stable over millions of years. As pointed out in [4], "the daisies are a surrogate for any kind of life that can affect the global temperature – and temperature could equally well be any other environmental variable that life cares about, for instance the oxygen concentrations, or pH." And indeed in their early analysis of the problem, Margulis and Lovelock already pointed out that, along with temperature and atmosphere composition, ocean acidity has been under feedback control [8].

In this work, we choose pH as a relevant environmental parameter that can be used as an external, tuneable input. Thus, we instead consider an alternative that allows a straightforward approach that captures all the relevant feedbacks and allows for a microcosm/mesocosm implementation (Figure 1d) based on a synthetic microbial ecosystem where pH is tuned by two populations that will play the role of our daisies.

2. Methods

This paper delves into two crucial aspects concerning the definition of a synthetic microbial Daisyworld. Firstly, we explore a collection of genetic circuits linked to a two-strain consortia specially engineered to regulate the environmental $pH$. Drawing inspiration from recent research on the controlled management of industrial fermentation [29], we present a synthetic approach aligned with the regulatory feedback principles of the DWM. Secondly, we demonstrate how the well-established impacts of acidity on microbial growth can be mapped onto ecological network models, akin to those of the DWM based on temperature and albedo.

2.1. Microbial Daisy World: synthetic circuits

The success of implementing a $pH$-based synthetic DWM hinges on two key factors: (a) ensuring a strong alignment with the growth response assumptions originally formulated by Lovelock, which were based on local temperature, and (b) skillfully engineering microbial-environment interactions to effectively regulate local $pH$, simulating the feedback mechanisms depicted in Figure 2a. Extensive research has been conducted on the growth responses of microorganisms to $pH$, leading to the development of various mathematical models to characterize their behavior. Notably, the analysis of the relative growth rate against $pH$ reveals a distinctive inverted parabolic pattern [30], which harmonizes with the foundational assumptions of the DWM. These assumptions entailed a smooth curve with a single optimum and well-defined limits. Empirical data from diverse species, such as E. coli or Listeria sp., demonstrate symmetric functional responses within a $pH$ range of $[p{H}_m,p{H}_M]$, where $p{H}_m$ and $p{H}_M$ represent the zero-growth limits of the fitted curve.

While different mathematical models, such as the Presser [31] or Lambert-Pearson [30,32] models, have been proposed, they all exhibit a near-parabolic behavior. Consequently, the effects of $pH$ closely align, in mathematical terms, with the control space assumptions of the DWM. This congruence underscores the potential feasibility of establishing a synthetic microbial Daisyworld based on $pH$ regulation.

Is it possible to design a synthetic consortium comprising two strains that can effectively control the acidity of the environment? Our aim is to create a pair of synthetic organisms capable of responding to changes in $pH$ in a manner that mirrors Lovelock’s concept, as illustrated in Figure 2a,b. In this scenario, the two engineered strains would function in opposing ways: although they share a common optimum $pH$," they would either increase or decrease the environmental acidity/alkalinity levels. This complementary behaviour would lead to a mutual self-regulation, ideally maintaining a constant environmental $pH$ value. In a broader context, this approach represents a specific application of metabolic engineering facilitated by synthetic biology [33]. Unlike traditional methods that involve continuous monitoring and manual addition of sterile bases and acids as needed [29,34], our case study focuses on developing a system capable of self-correcting deviations from the optimal $pH$. By doing so, we aim to eliminate the need for constant fine-tuning, making the process more efficient.

We propose the utilization of two engineered cell types (Figure 2c): the acid-producing and the base-producing strains. The first is the acid-producing strain, which incorporates the ldhA functional gene responsible for expressing lactate dehydrogenase. This enzyme facilitates the conversion of pyruvate to lactic acid, resulting in a decrease in pH [29]. Conversely, the base-producing strain employs the kivd functional gene, encoding $\alpha$-ketoisovalerate decarboxylase. This enzyme is involved in the decarboxylation of branched-chain $\alpha$-keto acids derived from branched-chain amino acids transamination into aldehydes [35]. In simpler terms, it catalyzes the conversion of 2-keto-acid to aldehyde, leading to an overproduction of ammonia, subsequently reducing the medium’s acidity by forming ammonium ions. While other genes like glsA or gadA could have been considered to catalyze ammonium ion production through alternative processes [29], using them might result in lower yields.

To ensure robustness in our experiments, we propose employing the Escherichia coli knockout strain with a deletion in the glutamine synthetase gene (Keio collection JW3841-1), yielding the ΔglnA genotype. By blocking ammonia re-uptake through this deletion, we can enhance the performance of the base-producing strain. Additionally, we recommend a deletion in ΔlldP, which encodes an inner membrane permease involved in lactate uptake [36,37]. This modification further supports the function of the acid-producing strain, ensuring its efficiency in the system.

Regarding promoter usage, we propose employing pTac in both genetic devices, as it is commonly used to control and overexpress recombinant proteins. However, any controllable promoter without leakiness and with high fold change expression could also be suitable. If an imbalance is detected in the production of both ldhA and kivd genes, the genetic devices could be easily modified to address the issue by using different controllable promoters, such as pTet-ldhA or pTet-kivd.

Furthermore, to facilitate real-time dynamics of the two strains and track the experimental process, fluorescent reporter genes would be used. We suggest using pH-sensitive fluorophores such as green and red fluorescent proteins. The green fluorescent protein GFP-pHlorine and the red fluorescent protein mCherryEA have been identified as ratiometric pH sensors, where the protonation of the chromophore is pH-dependent. They exhibit an approximately 8-fold increase in expression with increasing pH values, ranging from 5 to 9 [38,39].

2.2. Microbial Daisy World: two-species model

Using the known parabolic shape of relative growth responses to acidity, a simple mathematical model can be build under the assumption that the two synthetic strains described above can be engineered. The model now requires taking into account the fact that changes can be observed on two, connected scales. One is described (as in the WLM) in terms of relative cell populations. The second is a cell-level dynamics associated to the pH balances between the intracellular concentration and the one they perceive from their local environment.

The equations that are used to describe the cell dynamics are a coupled set of equations describing the growth of each cell type, similar to standard replicator equations from population genetics:

$$\begin{array}{c}\hfill \frac{d{X}_a}{dt}=X_a\left[\varphi \left(\overrightarrow{X}\right)\beta \left(p{H}_a\right)-\delta \right]\end{array}$$

(4)

$$\begin{array}{c}\hfill \frac{d{X}_b}{dt}=X_b\left[\varphi \left(\overrightarrow{X}\right)\beta \left(p{H}_b\right)-\delta \right]\end{array}$$

(5)

where $X_a$ is the concentration of acid-producing cells, $X_b$ is the concentration of base-producing cells, $\varphi \left(\overrightarrow{X}\right)=1-(X_a+X_b)$ describes the negative feedback associated to the finite amount of available space for cells to occupy, and $\delta$ is the dilution rate. The functional form of $\beta$ describing the change in growth rate depending on the perceived $pH$ is here chosen as

$$\beta \left(p{H}_x\right)=1-\frac{{(p{H}_{opt}-p{H}_x)}^2}{{(p{H}_{opt}\pm p{H}_{lim})}^2}$$

(6)

(other choices gave similar results, provided that the dependence is one-humped). It describes a symmetric, single-peaked function with its maximum output located at $p{H}_x=p{H}_{opt}$, and positive output for an input in the range $p{H}_x\in \left(p{H}_{opt}\pm p{H}_{lim}\right)$. The perceived $pH$ depends on the internal level of acidity/alkalinity that each cell experiences on a local scale. Thus we have:

$$\begin{array}{cccc}\hfill p{H}_a& =p{H}_f+a_i,\hfill & \hfill p{H}_b& =p{H}_f-b_i,\hfill \end{array}$$

(7)

Here, both $a_i$ and $b_i$ stand for the acid or base at the individual cell level and require specific dynamical assumptions (see below). $p{H}_f$ is described as the free $pH$ in the media, that depends on the external input $p{H}_{in}$ and the action of each population.

$$\frac{dp{H}_f}{dt}=p{H}_{in}+X_a{a}_e-X_b{b}_e-\delta p{H}_f$$

(8)

Here, $a_e$ and $b_e$ stand for the external level of acidity/alkalinity is able to excrete to the external surrounding media (see below). Furthermore, if we (reasonable) assume that the dynamics of protocol concentration is fast, then we can use $dp{H}_f/dt\approx 0$ and in such scenario we have:

$$p{H}_f=p{H}_{in}+X_a{a}_e-X_b{b}_e$$

(9)

Additionally, we need to take in account the microscopic balances at the level of individual cells. Therefore in the molecular level we can assume stoichiometry ruling the internal levels of acid or base production of each cell and the extracelular transport rates. We have the following reactions.

(10)

Here, as has been previously pointed, $a_e$ and $b_e$ stand as the external concentration of acid or base in each type of cell, and $a_i$ and $b_i$ are the amount of acid or base that the cell produces inside the membrane. The production rate is $\gamma$. The interchange ratios between the inner cell and the exterior are k. The differential equations according to the reactions are:

$$\begin{array}{cc}\hfill \frac{d{a}_i}{dt}& =\gamma -k_1{a}_i+k_2{a}_e\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill \frac{d{a}_e}{dt}& =k_1{a}_i-k_2{a}_e-\delta a_e\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill \frac{d{b}_i}{dt}& =\gamma -k_e{b}_i+k_4{b}_e\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill \frac{d{b}_e}{dt}& =k_3{b}_i-k_4{b}_e-\delta b_e\hfill \end{array}$$

(14)

It can be assumed that those reactions happening at the molecular level, in fact, are much faster than the population dynamics, therefore considering $d\mathsf{\Gamma }/dt=0$ for each variable $\mathsf{\Gamma }\in \{a_i,b_i,a_e,b_e\}$. In this case the equilibrium points can be described as constants. If we further assume all k as equal, we can pack $k-\delta /k$ as $\omega$.

$$\begin{array}{cccc}\hfill a_i^{*}& =\frac{\gamma (k-\delta )}{\delta k}=\omega \frac{\gamma }{\delta }\hfill & \hfill a_e^{*}& =\frac{\gamma }{\delta }\hfill \end{array}$$

(15)

$$\begin{array}{cccc}\hfill b_i^{*}& =\frac{\gamma (k-\delta )}{\delta k}=\omega \frac{\gamma }{\delta }\hfill & \hfill b_e^{*}& =\frac{\gamma }{\delta }\hfill \end{array}$$

(16)

Using the fast-relaxation assumption, we can now proceed to the analysis of the system dynamical patterns.

3. Results

Using the previous mathematical model approach (and under the assumptions made above) it is possible now to study the expected dynamical behaviour of our proposed microbial Daisyworld. We first consider the two-species scenarios where a consortium of two microorganisms acts on pH and secondly the more general scenario where (a) a parasitic species is added and (b) multiple species are introduced.

3.1. Two-species consortium

Figure 3 illustrates the behaviour of the synthetic DW design, using the set of equations described above. Here, two species with the same physiological preferences in $pH$ but opposite effects on the surrounding pH establish a wide zone of homeostasis. Similar to the classic work of Watson and Lovelock with temperature, the key parameter ($pH$) also increases in the environment where both species coexist. The surface in Figure 3a fully captures the presence and range of the homeostatic self-regulation achieved by our system. The three axes correspond to the external $pH$ input ($p{H}_{in}$), the actual $pH$ in the media ($p{H}_f$) and the $\gamma$ parameter that gives age rates of the acid/base production. For very small $\gamma$ values, the strains have no effect on the media pH and thus $p{H}_f=p{H}_{in}$.

As $\gamma$ increases, self-regulation emerges leading to a bell-shaped curve. This phenomenon can be observed more clearly in Figure 3b, which represents a two-dimensional graph showing $p{H}_f$, against $p{H}_{in}$, bounded by two extreme scenarios that would lead to ecosystem collapse (gray areas). If species have no effect on the surrounding $pH$ or there are simply no species alive, the result is a straight line representing both pH measures, therefore producing an uncontrolled environment. When a sufficiently high production rate is reached, homeostasis around the optimum appears.

In Figure 3c, we plot the corresponding equilibrium populations against input $pH$. The involved strains undergo a switch, therefore inducing homeostasis. The base production species ($X_b$) grows rapidly at low $p{H}_{in}$ since it is able to thrive in acidic environments by alkalizing its surroundings. When the $pH$ is high enough, the acidifying strain experiences a boost and out-competes the first strain until the pH becomes too high, resulting in the collapse of the system. Due to their presence, the $pH$ is regulated around seven within an effective range $3\le p{H}_{in}\le 12$, allowing both species to grow within this span. These results are fully consistent with the original DWM and support our proposal that an synthetic, small-scale implementation of planetary regulation can be formulated. How robust is this result? In the next two sections we answer this question under two distinct and relevant scenarios.

3.2. Two-species plus parasites

The conceptual framework of the DWM is grounded in the presence of a cooperative control of the environmental fluctuations. Such control is operated, as described above, by a two-species consortium. The effective outcome of the interactions between the members of the consortium is the stabilization of both populations, although submitted to a marked bias related to the pH input value. How robust is this mechanism of environmental stabilization? One way of answering this question is to consider the introduction of a parasitic component. Parasites are known to destabilise cooperative systems, sometimes pushing the population to the extinction threshold [40,41,42]. Such scenario has been tested in the context of synthetic biology, using a mixotrophic consortium growing in a 2D surface along with a parasitic strain (Amor et al., 2019).

What is the effect of an added parasitic strain? To answer this question, we consider an extended, three-species model where the new set of equations reads as

$$\begin{array}{c}\hfill \frac{d{X}_a}{dt}=X_a\left[\varphi \left(\overrightarrow{X}\right)\beta \left(p{H}_a\right)-\delta \right]\end{array}$$

(17)

$$\begin{array}{c}\hfill \frac{d{X}_b}{dt}=X_b\left[\varphi \left(\overrightarrow{X}\right)\beta \left(p{H}_b\right)-\delta \right]\end{array}$$

(18)

$$\begin{array}{c}\hfill \frac{d{X}_c}{dt}=X_c\left[\varphi \left(\overrightarrow{X}\right)\beta \left(p{H}_c\right)\alpha -\delta \right]\end{array}$$

(19)

Here $X_a$ and $X_b$ are the two previously defined strains and $X_c$ stands for new strain. The logistic term now reads $\varphi \left(\overrightarrow{X}\right)=1-(X_a+X_b+X_c)$. A growth rate parameter $\alpha \in [0,2]$ has been introduced in $d{X}_c/dt$ in order to tune the relative advantage of the extra species in relation to the pair of regulating species. The rest of parameters remain the same. In order to define a parasitic interaction, the functional form of the $\beta$ term for $X_c$ reads as:

$$\beta \left(p{H}_c\right)=1-\frac{{(p{H}_{opt}-p{H}_f)}^2}{{(p{H}_{opt}\pm p{H}_{lim})}^2}$$

(20)

Therefore, as defined, the added species takes advantage of the two other strains, but (as it must be for a parasite) has no effect on the free $pH$ and thus makes no contribution to regulation.

In Figure 4, we summarize the results of our theoretical model. The surface shown in Figure 4a dramatically illustrates the fact that the parasite does not undermine the homeostatic effect; instead, it enhances it, enabling the parasite to thrive in an environment with increasing pH levels. In Figure 4a we clearly appreciate the presence of a very broad domain of regulation.

The section defined for the neutral case $\alpha =1$ (vertical, grey plane) is shown in (c) along with the corresponding populations (right column). The homeostasis region appears mildly perturbed, and increases with $\alpha$ (see surface 4a). The $p{H}_{in}$ vs $Population$ graph shows that the parasite is able to thrive almost reaching carrying capacity when all share optimal $pH$. What if they don’t share an optimal $pH$? Figure 4b shows the robustness of the homeostatic effect. Here the parasite $p{H}_{op}$ ranges from 3 to 11 and $\alpha$ has been set to one. Obviously, when $p{H}_{op}=7$ we have the same plane ${\sum }_1$. When the $p{H}_{op}$ is acidic or alkaline (planes ${\sum }_2$ and ${\sum }_3$ respectively), the perturbation in the homeostasis happens at those specific ranges. Plane ${\sum }_2$ (Figure 4b,d) experiences a delay of the regulatory effect diminishing the presence of $X_b$. The inverse effect is expected in ${\sum }_3$ (Figure 4b,e). In summary, the introduction of a parasitic component does not (at all) have a detrimental effect on self-regulation. Our model predict instead that such extra population would favour the homeostatic patterns shown by the DWM.

3.3. Multispecies SDW: order and chaos

Low-dimensional ecosystems made of two or three species provide the simplest limit cases where general principles can be derived. As it occurs with genetic switches or simple oscillatory systems, a few components coupled through a nonlinear set of interactions can display robust behaviour [43]. Here too, the two-species example is our proof of concept that planetary regulation scenarios can be properly represented at small scales. Let us consider a general case study involving a multispecies consortium. In third example, we aim to explore the capacity for self-organization within a community where members have different $p{H}_{op}$ and $p{H}_{lim}$ values that we generate at random, as well as their different effects on pH. The question we seek to answer is whether the community can persist while stabilising the environmental pH. It is important to notice that multi species population dynamics under t nonlinear regimes explored here are likely to develop oscillations and chaos [44]. Will homeostasis also emerge in these cases? How is it affected by species diversity?

To address these questions, a new set of equations can be used. Those describing cell populations are already familiar to us:

$$\frac{d{X}_i}{dt}=X_i\left[\varphi \left(\overrightarrow{X}\right)\beta \left(p{H}_i\right)-\delta \right]+ϵ$$

(21)

Here, $X_i$ represents the concentration of pH-modifying bacteria, where $i=1,2,...,n$. Here $\varphi \left(\overrightarrow{X}\right)=1-{\sum }_j^n{X}_j$ is the negative feedback associated with finite resources, $\delta$ represents the dilution rate, and an immigration factor $ϵ$ is introduced. The growth rate $\beta$ has the same functional form as before, but now extended to multiple species:

$$\beta =1-\frac{1}{{\Delta }_H^2}{\left[p{H}_{opt}-p{H}_f+\frac{\omega }{\delta }\sum _{i=1}^n{\gamma }_i\right]}^2$$

(22)

Here, ${\Delta }_H=p{H}_{opt}-p{H}_{max}$, and $\gamma$ represents the production rate of acid or base, which can be positive or negative depending on whether it is an acid or base producer. The range of $\gamma$ is given by $\gamma \in \pm [{\gamma }_{min},{\gamma }_{max}]$. The dynamics for $p{H}_f$ experience a slight change and can be described as:

$$\frac{dp{H}_f}{dt}=\left(\sum _{i=1}^n{X}_i\frac{{\gamma }_i\omega }{\delta }\right)\left(\frac{p{H}_f}{b}(2b-p{H}_f)\right)$$

(23)

Here, the expression differs from our previous models First, there is no $p{H}_{in}$, indicating that the environment does not change due to external factors; rather, the community itself is the only changing factor. Additionally, following [45], since multiple factors (species) influence the pH, a saturation term has been introduced that ensures that $p{H}_f$ remains within the range of $[0,2b]$ (here $b=7$).

In order to characterise the general patterns of organization emerging from this high-dimensional set of equations, a statistical analysis of their dynamical behaviours has been performed, using different randomly generated communities. The analysis allows to see the frequency of different dynamical states as a function of diversity (Figure 4a). Following previous methods developed elsewhere [46], we use the number of species as a parameter to determine the likelihood to observe our final community in one of these three dynamical states: (a) stable point attractor, (b) oscillations and chaos and (c) collapse. These correspond to communities achieving populations that remain constant over time (after some transient), populations that exhibit deterministic fluctuations (either periodic, quasi periodic or chaotic) and finally those that experience extinction after a short initial instability.

Specifically, we used $N=50$ different consortia are simulated for each community size (from $n=1$ to $n=30$ species), where each species parameters has been randomly generated within the intervals $p{H}_{op}\in [4,9]$, $p{H}_{lim}=p{H}_{op}+[1,4]$ and $\gamma \in [0.05,0.2]$, $\delta =0.1$, using $\omega =1$ and a small immigration factor of $0.001$. All species are set in an initial $pH=7$ environment.

For each sampled model ecosystem, we discard a long transient to avoid non-stationary effects. Oscillatory or chaotic behaviour are found by looking at the amplitude variations of the time series of any of the species involved, and we do not distinguish among the different kinds of dynamics. Examples of both point attractors and fluctuating populations, along with the corresponding $p{H}_f$ time series are displayed in figures (b) and (c), respectively. We also display three-dimensional reconstructions of the underlying attractors. These are obtained by plotting the dynamics of one chosen species with a population $X\left(t\right)$ on a delayed 3-d space using the vectors $\left(X\right(t),x(t+T),x(t+2T)$ using a time delay T [47,48].

The observed fractions of dynamical states as a function of diversity allow us to formulated a strong prediction: despite the dominant presence of communities experiencing collapse or displaying oscillatory behavior, increasing numbers of species have a suppressor effect, favouring point attractors. In other words: species diversity enhances the desired self-regulation effects associated with our original SDWM scenario. Further theoretical model should explore the origins of such diversity-stabilizing pattern.

4. Discussion

The Daisyworld model provides a rationale for a stable self-organized biosphere resulting from a feedback between living beings and their environment, which they modify in predictable ways. Despite (or perhaps because of) its simplicity, the DWM has been instrumental in developing the field of Earth Systems Science [26]. Can new synthetic ecosystems help further develop the field by allowing experimental testing on the test tube? Previous work on synthetic ecosystems, both in vivo and in silico, have considered communities that illustrate the success or failure of feedback control in microbial ecosystems. Examples include ecosystem-level nutrient recycling in the Flask model [49] or the potential for collapse (ecological suicide) in microbial communities where sustained modification of environmental pH by the microbes can end in their extinction [45]. Our proposal instead is that of an engineered biological system capable of self-adjusting itself under a given range of external parametric conditions, actively acting on stabilising a global environmental driver. This, in our view, represents a novel concept.

Synthetic biology has been developing as a unique way to interrogate nature while exploring alternative designs not present in the biological world. One particularly important path involves the construction of synthetic counterparts of well known ecological systems. Over the last decade such developments have been extended into some of the classic problems of theoretical ecology, from tipping points to multispecies stability. In this paper, we propose to move a step forward by mapping the modelling planetary self-regulation problem into our test-tube analogy. Our proposed system will allow experimental testing of the DWM proposal. Along with the standard two-component species regulation, the generality of the design is supported by the robustness of the system under the presence of a parasite, which typically acts as a destabiliser force in cooperative systems. Instead, in the SDW system parasites do not change the potential for pH self-regulation. Moreover, an important, specific prediction from our work is that diverse communities should act in a self-organized way by supressing wide fluctuations and guaranteeing the convergence to a global controlled pH.

There are several potential extensions of our work that are worth exploring. On one hand, the regulatory nature of the system described here can be extended to other systems beyond the ecological context considered here. One avenue is the potential of terraformation scenarios based on synthetic biology, where cooperative consortia might be required for a successful spread over large scales [51]. Another one concerns those physiological systems (such as glucose regulation by glucagon and insulin) that can be described in similar terms [52,53]. These "rein control" systems could be constructed synthetically and used to target given regulation goals within a model organism. Finally, we have not included a major player in Lovelock’s picture: evolutionary dynamics. The role played by evolution has been a matter of intense discussion within the context of the DWM [4,54] and we may wonder if some of these debates could be settle by evolving synthetic communities. Extensions of our model approach could guide future developments towards this goal.