
==== Front
J R Soc Interface
J R Soc Interface
RSIF
royinterface
Journal of the Royal Society Interface
1742-5689
1742-5662
The Royal Society

38321922
10.1098/rsif.2023.0585
rsif20230585
10041816179Life Sciences–Earth Science interface
Research Articles
A synthetic microbial Daisyworld: planetary regulation in the test tube
A synthetic microbial Daisyworld: planetary regulation on the test tube
Maull Victor Conceptualization Formal analysis Investigation Methodology Visualization Writing – review & editing 1 2
Pla Mauri Jordi 1 2
Conde Pueyo Nuria 2 3
http://orcid.org/0000-0001-6974-1008
Solé Ricard Conceptualization Formal analysis Funding acquisition Investigation Methodology Resources Supervision Validation Writing – original draft Writing – review & editing ricard.sole@upf.edu
1 2 4
1 Institució Catalana de Recerca i Estudis Avançats, Psg Lluis Companys, Barcelona, Spain
2 Complex Systems Lab, Universitat Pompeu Fabra, Barcelona 08003, Spain
3 EMBL Barcelona, European Molecular Biology Laboratory (EMBL), Barcelona 08003, Spain
4 Santa Fe Institute, 1399 Hyde Park Road, Santa Fe, NM 87501, USA
7 2 2024 Feburary 7, 2024
2 2024
21 211 202305857 10 2023 October 7, 2023
12 1 2024 January 12, 2024
© 2024 The Author(s)
2024
https://royalsociety.org/-/media/journals/author/Licence-to-Publish-20062019-final.pdf https://royalsociety.org/journals/ethics-policies/data-sharing-mining/ Published by the Royal Society. All rights reserved.

The idea that the Earth system self-regulates in a habitable state was proposed in the 1970s by James Lovelock, who conjectured that life plays a self-regulatory role on a planetary-level scale. A formal approach to such hypothesis was presented afterwards under a toy model known as the Daisyworld. The model showed how such life-geosphere homeostasis was an emergent property of the system, where two species with different properties adjusted their populations to the changing external environment. So far, this ideal world exists only as a mathematical or computational construct, but it would be desirable to have a real, biological implementation of Lovelock’s picture beyond our one biosphere. Inspired by the exploration of synthetic ecosystems using genetic engineering and recent cell factory designs, here we propose a possible implementation for a microbial Daisyworld. This includes: (i) an explicit proposal for an engineered design of a two-strain consortia, using pH as the external, abiotic control parameter and (ii) several theoretical and computational case studies including two, three and multiple species assemblies. The special alternative implementations and their implications in other synthetic biology scenarios, including ecosystem engineering, are outlined.

Daisyworld
, homeostasis
, Earth systems science
, synthetic biology
, terraformation
Agència de Gestió d'Ajuts Universitaris i de Recerca http://dx.doi.org/10.13039/501100003030 AGAUR 2021 SGR 00751JP
==== Body
pmc1. 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 [1,2]. 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 [3]. The deep connections between environment and life leave a mark in the geological past. How has life influenced climate and vice versa? As pointed out by Vernadsky [4], the emergence of life fundamentally transformed the geosphere. One specially interesting observation made by geochemist James Lovelock in the 1970s was the realization 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 an inhabitable planet [5]. In general, positive feedback loops can drive a planet to extreme steady states, from boiling temperatures to a snowball [3,6]. 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 [7–10]. To make this point more explicit, Watson & Lovelock proposed in 1983 a simple model of planetary regulation known as the Daisyworld model (DWM) [11,12]. 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, the DWM has become a canonical model for Earth system science [13] and has been used to explore a very diverse range of problems, from tipping points [14,15], to evolutionary dynamics [16,17] to climate bistability in exoplanets [18]. In contrast to other problems, the planetary scale makes it 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 [19,20], predator–prey systems [21] or even multispecies assemblies [22–24]. 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 defining an ecological model, namely (see [25] and references therein): 1.1 dxkdt=xk(βk(T)[1−∑ j=1nx j]−δk)

where xk (k = 1, …, n) is the population size of the kth species. The parameter βk(T) is a temperature-dependent growth rate, whereas δk is the corresponding decay (death) rate. A logistic saturation is introduced by the 1−∑ jx j term.

The basic feedback loops involved are represented in figure 1a for the standard two species (n = 2) model. Increasing luminosity can trigger the growth of both kinds of vegetation, at rates indicated by βk(T) (with k = 1, 2), here using k to refer to the different kinds of daisies, as a generalization, further on we specify which daisies are at play. In the best known (and simpler) version, two species are considered and the following model was explored, for a ‘world’ involving black (xb) and white (xw) daisies 1.2 dxwdt=xw(βw(T)[1−xw−xb]−δw)

and 1.3 dxbdt=xb(βb(T)[1−xw−xb]−δb)

where [1 − xw − xb] introduces the available bare soil that both species can occupy. In figure 1a, we sketch the basic model components using a spatial version [26–28], with dark, white and grey 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 βi(T) factors (with i = b, w), for black and white, defined by means of a single-humped function involving an optimum. The details of this model can be found elsewhere [11,12] 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 cool their environment locally. The area covered by each species contributes to the albedo and affect local temperatures in such a way that a global temperature regulation emerges, 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. Figure 1. The conceptual feedbacks in the Watson–Lovelock model, which is the canonical formulation of the DWM. Using a two-dimensional surface (a), increasing levels of solar luminosity L trigger the growth (after a threshold) of two populations of plants, indicated as B and W (white and black squares, grey squares stand for bare soil S, see inset on the left bottom, using a zoom on the area indicated). These are identified as black and white daisies, respectively, in the DWM. Opposite feedbacks emerge from the effects of albedo by the two kinds of daisies. As a consequence of these nonlinear couplings, as shown in (b), the planet temperature can be stabilized (instead of just simply growing with L, grey line) for a wide range of L values, thus indicating a homeostatic response due too the biosphere–climate system. Such stabilization is obtained by means of population arrangements between W and B states (c). In this paper, an equivalent system is proposed (d) using a bioreactor where an external input is also present (Hin+) that would increase the pH of the medium, unless feedback controls are present. In this case, two different strains of bacteria able to increase or decrease the pH would replace the daisies.

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 theoretical feasibility of constructing a synthetic Daisyworld.

Creating a living surrogate of the original DWM could help explore 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 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 also been controlled at the planetary level. One of them is acidity [9]: despite the tendency towards acidification associated with increasing oxidation of the atmosphere, the mean pH of the oceans has been remarkably stable over the Phanerozoic [29–32]. As pointed out in [3], '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 & Lovelock already pointed out that, along with temperature and atmosphere composition, ocean acidity has been under feedback control [9].

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 [33], we present a novel synthetic microbial consortium 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 Daisyworld: synthetic circuits

The success of implementing a pH-based synthetic DWM hinges on two key factors: (i) ensuring a strong alignment with the growth response assumptions originally formulated by Lovelock, which were based on local temperature and (ii) skilfully 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 behaviour. Notably, the analysis of the relative growth rate against pH reveals a distinctive inverted parabolic pattern [34], 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 Escherichia coli or Listeria sp., demonstrate symmetric functional responses within a pH range of [pHm, pHM], where pHm and pHM represent the zero-growth limits of the fitted curve [34]. Figure 2. The synthetic microbial Daisyworld. The logic of the zero-dimensional (non-spatial) DWM is summarized in (a) in terms of the interactions between planetary temperature and the distinct role played by the two kinds of daisies (B and W). Using our framework, where acidity would be the controlled variable. In this context, the acidity level within the medium influences cellular replication (self-replication loops). Simultaneously, both species reciprocally influence the acidity of the medium. A similar logic of feedback loops can be described (b) where now two microbial populations would also reduce or increase local acidity. The whole design, including the corresponding genetic constructs, is depicted in (c).

While different mathematical models, such as the Presser [35] or Lambert–Pearson [34,36] models, have been proposed, they all exhibit a near-parabolic behaviour. 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 designed strains 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 [37]. Unlike traditional methods that involve continuous monitoring and manual addition of sterile bases and acids as needed [33,38], 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 growing in a chemostat supplying a constant flow of nutrients and removing a constant proportion of cells. The acid-producing strain 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 external pH [33]. Conversely, the base-producing strain employs the kivd functional gene, encoding α-ketoisovalerate decarboxylase. This enzyme is involved in the decarboxylation of branched-chain α-keto acids derived from branched-chain amino acids transamination into aldehydes [39]. In simpler terms, it catalyses 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 catalyse ammonium ion production through alternative processes [33], using them might result in lower yields.

To ensure robustness in the experiments, we propose employing the E. 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 [40,41]. 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 changing one device to a different controllable promoter, such as pTet-ldhA or pTet-kivd. Our proposed model involves the production of strong acids/bases, but the generation of such highly reactive chemicals would be detrimental to the survival of the bacteria. Therefore, we suggest using an experimental system with production of weaker, autonomously controlled reactants. This approach aims to minimize differences between the model and the experiments.

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-pHluorin 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 eightfold increase in expression with increasing pH values, ranging from 5 to 9 [42,43].

This synthetic consortium would grow in a bioreactor environment (figure 1c) with a stable supply of fresh media ensuring nitrogen availability (required for a proper ammonia production).

2.2. Microbial Daisyworld: 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 in terms of relative cell populations. The second is a cell-level dynamics associated with 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: 2.1 dXadt=Xa[ϕ(X)β(pHa)−δ]

2.2 dXbdt=Xb[ϕ(X)β(pHb)−δ]

where Xa is the concentration of acid-producing cells, Xb is the concentration of base-producing cells, ϕ(X)=1−(Xa+Xb) describes the negative feedback associated with the finite amount of available space for cells to occupy and δ is the dilution rate. The functional form of β describing the change in growth rate depending on the perceived pH is here chosen as 2.3 β(pHx)=1−(pHopt−pHx)2(pHopt±pHlim)2

(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 pHx = pHopt, and positive output for an input in the range pHx ∈ (pHopt ± pHlim). It guarantees that all factors affecting cell growth that are pH-dependent, including the availability of nutrients, decreasing as pH strays from the optimum. The perceived pH depends on the level of acidity/alkalinity that each cell experiences on a local scale. Thus, we have: 2.4 pHa=pH f+ωaepHb=pH f−ωbe 

where both ae and be stand for the acid or base produced at the individual cell level and require specific dynamical assumptions (see below), and ω is the sensitivity to the produced substance. pHf is described as the free pH in the media, that depends on the external input pHin and the action of each population. We only take into account the external pHf instead of the usual ΔpH as it is assumed that the cell’s internal pH will remain almost constant [44]. 2.5 dpH fdt=pHin+Xaae−Xbbe−δpH f

Here, ae and be stand for the external level of acidity/alkalinity a cell is able to excrete to the external surrounding media (see below). Furthermore, if we (reasonable) assume that the pHf dynamics are fast, and thus we can use dpHf/dt ≈ 0 and in such a scenario we have: 2.6 pHf=pHin+Xaae−Xbbe

Additionally, we need to take in account the microscopic balances at the level of individual cells. Therefore, at 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. 2.7 ∅ →γ ai ⇌k2k1 ae →δ ∅∅ →γ bi ⇌k4k3 be →δ ∅

Here, as has been previously pointed, ae and be stand as the external concentration of acid or base in each type of cell, and ai and bi are the amount of acid or base that the cell produces inside the membrane. The production rate is γ. The interchange ratios between the inner cell and the exterior are k. The differential equations according to the reactions are 2.8 daidt=γ−k1ai+k2ae

2.9 daedt=k1ai−k2ae−δae

2.10 dbidt=γ−kebi+k4be

2.11 dbedt=k3bi−k4be−δbe

It can be assumed that those reactions happening at the molecular level, in fact, are much faster than the population dynamics, therefore, considering dΓ/dt=0 for each variable Γ∈{ai,bi,ae,be}. In this case, the equilibrium points can be described as constants. 2.12 ai∗=γ(k−δ)δkae∗=γδ

2.13 bi∗=γ(k−δ)δkbe∗=γδ

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 the assumptions made, it is possible now to study the expected dynamical behaviour of our proposed microbial Daisyworld. We first consider the two-species scenario (A) where a consortium of two microorganisms acts on pH and secondly the more general scenario where (B) a parasitic species is added and (C) multiple species are introduced. For low-dimensional (two to three species) cases, the system dynamics are computed numerically through an iterative Euler method, ensuring stability for both species and pHf under each external pH input. The multispecies equations were numerically solved using a Runge–Kutta fourth-order approximation. In table 1, all the parameters and functions used in our models are summarized.

3.1. Two-species consortium

Figure 3 illustrates the behaviour of the synthetic Daisyworld design, using the set of equations described in section §2.2 of the Methods. 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 & 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 (pHin), the actual pH in the media (pHf) and the γ parameter that gives age rates of the acid/base production. For very small γ values, the strains have no effect on the media pH and thus pHf = pHin. Figure 3. Steady state surface for the synthetic Daisyworld are depicted using equations (2.1) and (2.2). The flat surface displayed in (a) corresponds to the self-regulation illustrating dynamics of the synthetic microbial Daisyworld. The relationship between pHin and pHf is presented against γ. When γ = 0, neither base nor acid production occurs, resulting in the absence of regulation and a linear increase in pHf with the external input. Conversely, as γ increases, a diverse range of controlled pH values emerges, expanding the homeostatic region of optimal pH. In (b), one section to the surface at γ = 0.04 is made, showing the broad domain of self-regulation surrounded by the collapse domains. The corresponding abundances of each synthetic strain are displayed in (c). Xa stands for acid-producing species and Xb for base-producing. Here, we used: γ ∈ [0, 0.05], δ = 0.01, ω = 0.5 and both species use pHopt = 7 and pHlim = 9.

As γ 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 pHf, against pHin, bounded by two extreme scenarios that would lead to ecosystem collapse (grey 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 (Xb) grows rapidly at low pHin 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 ≤ pHin ≤ 12, allowing both species to grow within this span. These results are fully consistent with the original DWM and support our proposal that a 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 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 destabilize cooperative systems, sometimes pushing the population to the extinction threshold [45–47]. Such a scenario has been tested in the context of synthetic biology, using a mixotrophic consortium growing in a 2D surface along with a parasitic strain [20].

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 3.1 dXadt=Xa[ϕ(X)β(pHa)−δ]

3.2 dXbdt=Xb[ϕ(X)β(pHb)−δ]

3.3 dXcdt=Xc[ϕ(X)β(pHc)α−δ]

Here Xa and Xb are the two previously defined strains and Xc stands for new strain. The logistic term now reads ϕ(X)=1−(Xa+Xb+Xc). A growth rate parameter α ∈ [0, 2] has been introduced in dXc/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 β term for Xc reads as 3.4 β(pHc)=1−(pHopt−pH f)2(pHopt±pHlim)2

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 maintains 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. Figure 4. Effects of parasitism (using equations (3.1)–(3.3)) on self-regulation domains. In (a), the relationship between pHin and pHf is displayed against α as a control parameter, weighting the parasite advantage relative to the regulating species. In (b), the relationship between pHin and pHf is shown for parasites with pHopt ∈ [3, 11] and pHlim = pHopt + 2.5 are introduced, with α = 1. The sections obtained from the vertical planes Σ cutting through both surfaces are depicted in (c–e), planes Σ1, Σ2 and Σ3, respectively. The three cases correspond (for α = 1) to parasite optimal values of (c) pHopt = 7, (d) pHopt = 4 and (e) pHopt = 10. The corresponding population abundances are depicted in the right column. The remaining parameters are the same as those used in figure 3.

The section defined for the neutral case α = 1 (vertical, grey plane) is shown in figure 4c along with the corresponding populations (right column). The homeostasis region appears mildly perturbed, and remains constant with α (see surface figure 4a). The pHin versus Population graph shows that the parasite is able to thrive almost reaching carrying capacity when all share optimal pH. What if they do not share an optimal pH? Figure 4b shows the robustness of the homeostatic effect. Here the parasite pHopt ranges from 3 to 11 and α has been set to 1. Obviously, when pHopt = 7 we have the same plane Σ1. When the pHopt is acidic or alkaline (planes Σ2 and Σ3, respectively), the perturbation in the homeostasis happens at those specific ranges. Plane Σ2 (figure 4b,d) experiences a delay in the regulatory effect, diminishing the presence of Xb. The inverse effect is expected in Σ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 demonstrates that the properties of a pH-regulated DWM are equivalent to those that have already been described for more commonly studied thermostatic DWMs. Table 1. Functions and parameters defined throughout the paper.

variable	description	
Xa, Xb, Xc	acid-producing, base-producing and parasite cell concentration	
ϕ(Xi)	function representing negative feedback due to limited space	
β(pHi)	growth rate as a function of the perceived external pH value	
pHopt, pHlim	optimal pH value and pH limit, respectively, at which cells can grow	
pHa, pHb	external pH experienced by acid-producing and base-producing cells	
pHf	pH in the media	
pHin	external input pH	
ae, be	steady state external levels of acid and base produced by a cell	
ai, bi	steady state internal levels of acid and base produced by a cell	
γ	production rate of acid/base inside a cell	
k1, k2, k3, k4	interchange rates between inner cell and exterior for acid and base	
δ	dilution rate	
ω	cell sensitivity to the produced acid/base	
α	growth rate parameter for parasitic cells	
ϵ	immigration rate for each species in the multispecies model	
b	midpoint of the possible pH range in the multispecies model	

3.3. Multispecies synthetic Daisyworld: 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 [48]. Here too, the two-species example is our proof of concept that planetary regulation scenarios can be properly represented at small scales. What about a multispecies scenario? Previous work on the DW model has also considered the impact of biodiversity by introducing a range of competitor species displaying different parameters, also including different trophic levels [49]. Here, we aim to explore the capacity for self-organization within a community where members have different pHopt and pHlim 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 stabilizing the environmental pH. It is important to note that multi-species population dynamics under nonlinear regimes explored here are likely to develop oscillations and chaos [50]. 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: 3.5 dXidt=Xi[ϕ(X)β(pHi)−δ]+ϵ

Here, Xi represents the concentration of pH-modifying bacteria, where i = 1, 2, …, n. Here ϕ(X)=1−∑ jnX j is the negative feedback associated with finite resources, δ represents the dilution rate, and an immigration factor ϵ is introduced. The growth rate β has the same functional form as before, but now extended to multiple species: 3.6 β=1−1ΔH2[pHopt−pH f+ωδ∑i=1nγi]2

Here, ΔH = pHopt − pHmax, and γ 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 γ is given by γ ∈ ±[γmin, γmax]. The dynamics for pHf experience a slight change and can be described as 3.7 dpH fdt=(∑i=1nXiγiωδ)(pH fb(2b−pH f))

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

As it occurs with most high-dimensional, non-spatial dynamical systems, complex dynamical states might often involve high-amplitude oscillations, hindering species from establishing a suitable environment before collapse [52,53]. To address this limitation of well-mixed models, one effective approach is the introduction of a dispersal term or immigration factor [54]. Both theoretical and empirical studies have demonstrated that moderate levels of dispersal serve as a primary mechanism for survival, disrupting high-amplitude population dynamics [54–56]. Furthermore, in challenging environments, evolutionary rescue has been observed as a means of population recovery [57]. In a mean-field model such as ours, this effect can be effectively simulated through immigration, thus serving as a proxy for both space and evolution, while maintaining the model in its minimal state [58]. All species are set in an initial pH = 7 environment.

In order to characterize 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 us to see the frequency of different dynamical states as a function of diversity (figure 5a). Following previous methods developed elsewhere [59], 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: (i) stable fixed point attractors, (ii) oscillations and chaos and (iii) 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. In figure 5a, we employed a polynomial fitting of degree 4 for circles (raw data) and degree 3 for squares (raw data) to enhance the visualization of the dynamical states. Additionally, (figure 5a, right axis), we show the fraction of the original species that survive for each initial community size. Figure 5. Statistical analysis and specific examples are illustrated for the multispecies case study. Three different final outcomes are obtained, as summarized in (a): extinction due to collapse and failure of regulation, oscillatory, quasi-periodic or chaotic fluctuations (black) and self-regulating communities displaying constant values after a transient. Although chaotic, oscillatory and collapse states are frequent at low-species numbers, they become rare as diversity increases. Raw data in squares and circles. Two examples of stable (b) and chaotic (c) communities are also shown. In both cases, we show the time series of populations, the reconstructed attractor of one of them, the pH time series and its reconstructed attractor are displayed. The population and pH dynamics for case (c) define strange, chaotic attractors. In addition, the mean percentage of the original species that survive for each initial community size is displayed in (a), in the right y axis (blue). Standard deviation is also plotted with error bars. Here, we used N = 200 replicas for each community size (up to n = 30 species), where parameters have been randomly generated within the intervals pHopt ∈ [4, 9], pHlim = pHopt + [1, 4] and γ ∈ [0.05, 0.2], δ = 0.1, using ω = 1 and a small immigration factor of ϵ=10−3 has been used.

For each sampled model ecosystem, we discard a long transient to avoid non-stationary effects. Oscillatory or chaotic behaviours 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 pHf time series are displayed in figure 5b,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(t) on a delayed three-dimensional space using the vectors X(t), x(t + T), x(t + 2T) using a time delay T [60,61].

Interestingly, the observed fractions of dynamical states as a function of diversity allow us to formulate a strong prediction: despite the dominant presence of fluctuating populations for high-dimensional ecosystems, increasing biodiversity has a suppressor effect, favouring point attractors. In other words, species diversity enhances the desired self-regulation effects. It is also worth mentioning that earlier higher diversity seems to increase the chances of a stable sub-community surviving an early upheaval period and going on to regulate their environment due to their complementary environmental impacts, as can be seen in figure 5a. This is in line with previous models on planetary self-regulation that have shown how community restructuring can lead to stability [62].

Is this homeostasis preserved when exogenous fluctuations of pH are present? Will fluctuating populations (such as those in the ‘oscillations’ phase of figure 5) cope with a time-dependent, noisy input? We have explored this problem by introducing a stochastic component in the time evolution of pHin. The multispecies synthetic Daisyworld not only self-organizes to configure its environment optimally, but also exhibits the capacity to buffer external pH perturbations, sustaining oscillatory homeostasis in the face of a fluctuating environment. The external perturbation takes the form of an additive (discrete) uniform white noise sampled at every time step, characterized by a zero mean, and magnitudes between [−5, 5] over an interval. In figure 6, we show several examples of these responses for four different community sizes associated with either chaotic or stable point attractors, respectively. In all these cases, it is found that the community responses display pH buffering and thus robust self-regulation. Further work should address the statistical patterns displayed by stochastic versions of our previous models. Figure 6. Homeostatic behaviour for multispecies models under stochastic variations in the pH input. Here, we have chosen for different conditions associated with both stable (point attractors) and oscillating communities (as defined from the phase space in figure 5a). Specifically, the cummulative input pHin fluctuates in time (blue lines) starting from pHin(0) = 7. The time series for populations (left) and pH values (right) are shown. All communities respond by adjusting pHf (black lines) which typically fluctuate within abounded range of values around the deterministic regulated value.

4. Discussion

The DWM 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 [13]. Can new synthetic ecosystems help further develop the field by allowing experimental testing in 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 [63] or the potential for collapse (ecological suicide) in microbial communities where sustained modification of environmental pH by the microbes can end in their extinction [64]. 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 stabilizing a global environmental driver. We propose a specific design for the genetic constructs required to self-tune the system. The one-humped nature of growth responses against pH displayed by microorganisms makes our candidate designs perfectly fit to match the original DWM assumptions. Moreover, under a multispecies context, biodiversity is a firewall to prevent the system from becoming de-regulated. This, in our view, represents a novel concept.

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 [65]. Another one concerns those physiological systems (such as glucose regulation by glucagon and insulin) that can be described in similar terms [66,67]. 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. Previous work on evolutionary dynamics in DWMs [68–71]. The role played by evolution has been a matter of intense discussion within the context of the DWM [3,72] 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.

Acknowledgements

We thank Daniel Amor for insight into synthetic multispecies communities. V.M. thanks Dimitrije Ivancić and Nastassia Knödlseder for useful discussions. R.S. thanks the late Patrick Marcos-Nikolaus for introducing him to Vladimir Vernadsky’s work. Special thanks to the Santa Fe Institute, where this research was done and to David Krakauer, who allowed us to work in the monastery on the night shifts.

Ethics

This work did not require ethical approval from a human subject or animal welfare committee.

Data accessibility

The code used for the simulations is available from the Zenodo repository: https://doi.org/10.5281/zenodo.10475383 [73].

Declaration of AI use

We have not used AI-assisted technologies in creating this article.

Authors' contributions

V.M.: conceptualization, formal analysis, investigation, methodology, visualization, writing—review and editing; J.P.M.: conceptualization, investigation, writing—review and editing; N.C.P.: conceptualization, supervision, writing—review and editing; R.S.: conceptualization, formal analysis, funding acquisition, investigation, methodology, resources, supervision, validation, writing—original draft, writing—review and editing.

All authors gave final approval for publication and agreed to be held accountable for the work performed therein.

Conflict of interest declaration

We declare we have no competing interests.

Funding

This work was funded by grant nos. FIS2016-77447-R MINECO/AEI/FEDER, and AGAUR 2021 SGR 00751. J.P.M. has been funded by the PRE2020-091968 grant from the Spanish government.
==== Refs
References

1. Knoll AH. 2015 Life on a young planet. Princeton, NJ: Princeton University Press.
2. Knoll AH, Nowak MA. 2017 The timetable of evolution. Sci. Adv. 3 , e1603076. (10.1126/sciadv.1603076)28560344
3. Lenton T, Watson A. 2013 Revolutions that made the Earth. Oxford, UK: Oxford University Press.
4. Vernadsky VI. 1998 The biosphere. New York, NY: Springer Science.
5. Mackwell SJ, Simon-Miller AA, Harder JW, Bullock MA, eds., 2014 Comparative climatology of terrestrial planets. Tucson, AZ: University of Arizona Press.
6. Corsetti FA, Olcott AN, Bakermans C. 2006 The biotic response to Neoproterozoic snowball Earth. Palaeogeogr. Palaeoclimatol. Palaeoecol. 232 , 114-130. (10.1016/j.palaeo.2005.10.030)
7. Lovelock JE. 1972 Gaia as seen through the atmosphere. Atmos. Environ. 6 , 579-580. (10.1016/0004-6981(72)90076-5)
8. Lovelock JE, Margulis L. 1974 Atmospheric homeostasis by and for the biosphere: the Gaia hypothesis. Tellus 26 , 2-10. (10.3402/tellusa.v26i1-2.9731)
9. Margulis L, Lovelock JE. 1974 Biological modulation of the Earth’s atmosphere. Icarus 21 , 471-489. (10.1016/0019-1035(74)90150-X)
10. Lovelock J. 2016 Gaia: a new look at life on earth. Oxford, UK: Oxford University Press.
11. Watson AJ, Lovelock JE. 1983 Biological homeostasis of the global environment: the parable of Daisyworld. Tellus B 35 , 284-289. (10.3402/tellusb.v35i4.14616)
12. Lovelock JE, Watson AJ. 1982 The regulation of carbon dioxide and climate: gaia or geochemistry. Planet. Space Sci. 30 , 795-802. (10.1016/0032-0633(82)90112-X)
13. Steffen W, Richardson K, Rockström J, Schellnhuber HJ, Dube OP, Dutreuil S, Lenton TM, Lubchenco J. 2020 The emergence and evolution of Earth system science. Nat. Rev. Earth Environ. 1 , 54-63. (10.1038/s43017-019-0005-6)
14. Ackland GJ, Clark MA, Lenton TM. 2003 Catastrophic desert formation in Daisyworld. J. Theor. Biol. 223 , 39-44. (10.1016/S0022-5193(03)00069-9)12782115
15. Wilkinson DM. 2003 Catastrophes on Daisyworld. Trends Ecol. Evol. 18 , 266-268. (10.1016/S0169-5347(03)00097-1)
16. Lenton TM, Lovelock JE. 2001 Daisyworld revisited: quantifying biological effects on planetary self-regulation. Tellus B 53 , 288-305. (10.1034/j.1600-0889.2001.01191.x)
17. Lenton TM, Daines SJ, Dyke JG, Nicholson AE, Wilkinson DM, Williams HT. 2018 Selection for Gaia across multiple scales. Trends Ecol. Evol. 33 , 633-645. (10.1016/j.tree.2018.05.006)30041995
18. Murante G et al. 2020 Climate bistability of Earth-like exoplanets. Mon. Not. R. Astron. Soc. 492 , 2638-2650. (10.1093/mnras/stz3529)
19. Shou W, Ram S, Vilar JM. 2007 Synthetic cooperation in engineered yeast populations. Proc. Natl Acad. Sci. USA 104 , 1877-1882. (10.1073/pnas.0610575104)17267602
20. Amor DR, Montanez R, Duran-Nebreda S, Solé R. 2017 Spatial dynamics of synthetic microbial mutualists and their parasites. PLoS Comput. Biol. 13 , e1005689. (10.1371/journal.pcbi.1005689)28827802
21. Balagadde FK, Song H, Ozaki J, Collins CH, Barnet M, Arnold FH, Quake SR, You L. 2008 A synthetic Escherichia coli predator–prey ecosystem. Mol. Syst. Biol. 4 , 87.
22. Mee MT, Wang HH. 2012 Engineering ecosystems and synthetic ecologies. Mol. Biosyst. 8 , 2470-2483. (10.1039/c2mb25133g)22722235
23. Johns NI, Blazejewski T, Gomes AL, Wang HH. 2016 Principles for designing synthetic microbial communities. Curr. Opin. Microbiol. 31 , 146-153. (10.1016/j.mib.2016.03.010)27084981
24. Hu J, Amor DR, Barbier M, Bunin G, Gore J. 2022 Emergent phases of ecological diversity and dynamics mapped in microcosms. Science 378 , 85-89. (10.1126/science.abm7841)36201585
25. Wood AJ, Ackland GJ, Dyke JG, Williams HT, Lenton TM. 2008 Daisyworld: a review. Rev. Geophys. 46 , RG1001. (10.1029/2006RG000217)
26. Adams B, Carr J, Lenton TM, White A. 2003 One-dimensional Daisyworld: spatial interactions and pattern formation. J. Theor. Biol. 223 , 505-513. (10.1016/S0022-5193(03)00139-5)12875827
27. Ackland GJ, Wood AJ. 2010 Emergent patterns in space and time from Daisyworld: a simple evolving coupled biosphere–climate model. Phil. Trans. R. Soc. A 368 , 161-179. (10.1098/rsta.2009.0203)19948549
28. Punithan D, Kim DK, McKay RB. 2012 Spatio-temporal dynamics and quantification of daisyworld in two-dimensional coupled map lattices. Ecol. Complex. 12 , 43-57. (10.1016/j.ecocom.2012.09.004)
29. Ridgwell A, Zeebe RE. 2005 The role of the global carbonate cycle in the regulation and evolution of the Earth system. Earth Planet. Sci. Lett. 234 , 299-315. (10.1016/j.epsl.2005.03.006)
30. Kump LR, Bralower TJ, Ridgwell A. 2009 Ocean acidification in deep time. Oceanography 22 , 94-107. (10.5670/oceanog.2009.100)
31. Pelejero C, Calvo E, Hoegh-Guldberg O. 2010 Paleo-perspectives on ocean acidification. Trends Ecol. Evol. 25 , 332-344. (10.1016/j.tree.2010.02.002)20356649
32. Honisch B et al. 2012 The geological record of ocean acidification. Science 335 , 1058-1063. (10.1126/science.1208277)22383840
33. Li C et al. 2020 Intelligent microbial cell factory with genetic pH shooting (GPS) for cell self-responsive base/acid regulation. Microb. Cell Fact. 19 , 1-13. (10.1186/s12934-020-01457-3)31898497
34. Lambert RJ. 2011 A new model for the effect of pH on microbial growth: an extension of the Gamma hypothesis. J. Appl. Microbiol. 110 , 61-68. (10.1111/j.1365-2672.2010.04858.x)20880208
35. Presser KA, Ratkowsky DA, Ross T. 1997 Modelling the growth rate of Escherichia coli as a function of pH and lactic acid concentration. Appl. Environ. Microbiol. 63 , 2355-2360. (10.1128/aem.63.6.2355-2360.1997)9172355
36. Lambert RJW, Pearson J. 2000 Susceptibility testing: accurate and reproducible minimum inhibitory concentration (MIC) and non-inhibitory concentration (NIC) values. J. Appl. Microbiol. 88 , 784-790. (10.1046/j.1365-2672.2000.01017.x)10792538
37. Lv X, Hueso-Gil A, Bi X, Wu Y, Liu Y, Liu L, Ledesma-Amaro R. 2022 New synthetic biology tools for metabolic control. Curr. Opin. Biotechnol. 76 , 102724. (10.1016/j.copbio.2022.102724)35489308
38. Kambale SD, George S, Zope RG. 2015 Controllers used in pH neutralization process: a review. Int. Res. J. Eng. Technol. 2 , 354-361.
39. Mikami Y, Yoneda H, Aoki W, Ueda M. 2017 Ammonia production from amino acid-based biomass-like sources by engineered Escherichia coli. AMB Express 7 , 83. (10.1186/s13568-017-0385-2)28429328
40. Núñez MF, Kwon O, Wilson TH, Aguilar J, Baldoma L, Lin ECC. 2017 Transport of L-Lactate, D-Lactate, and glycolate by the LldP and GlcA membrane carriers of Escherichia coli. Biochem. Biophys. Res. Commun. 290 , 824-829.
41. de la Plaza M, Peláez C, Requena T. 2009 Regulation of α-ketoisovalerate decarboxylase expression in Lactococcus lactis IFPL730. J. Mol. Microbiol. Biotechnol. 17 , 96-100. (10.1159/000178018)19033676
42. Reifenrath M, Boles E. 2018 A superfolder variant of pH-sensitive pHluorin for in vivo pH measurements in the endoplasmic reticulum. Sci. Rep. 8 , 11985. (10.1038/s41598-018-30367-z)30097598
43. Hartmann FSF, Weiss T, Shen J, Smahajcsik D, Savickas S, Seibold GM. 2022 Visualizing the pH in Escherichia coli colonies via the sensor protein mCherryEA allows high-throughput screening of mutant libraries. Msystems 7 , e00219-22. (10.1128/msystems.00219-22)35430898
44. Padan E, Zilberstein D, Schuldiner S. 1981 pH homeostasis in bacteria. Biochim. Biophys. Acta 650 , 151-166. (10.1016/0304-4157(81)90004-6)6277371
45. Smith JM. 1979 Hypercycles and the origin of life. Nature 280 , 445-446. (10.1038/280445a0)460422
46. Boerlijst MC, Hogeweg P. 1991 Spiral wave structure in pre-biotic evolution: hypercycles stable against parasites. Physica D 48 , 17-28. (10.1016/0167-2789(91)90049-F)
47. Sardanyes J, Solé R. 2007 Spatio-temporal dynamics in simple asymmetric hypercycles under weak parasitic coupling. Physica D 231 , 116-129. (10.1016/j.physd.2007.04.009)
48. Sneppen K. 2014 Models of life. Cambridge, UK: Cambridge University Press.
49. Lovelock JE. 1992 A numerical model for biodiversity. Phil. Trans. R. Soc. Lond. B 338 , 383-391. (10.1098/rstb.1992.0156)
50. Munch SB, Rogers TL, Johnson BJ, Bhat U, Tsai CH. 2022 Rethinking the prevalence and relevance of chaos in ecology. Ann. Rev. Ecol. Evol. Syst. 53 , 227-249. (10.1146/annurev-ecolsys-111320-052920)
51. Ratzke C, Gore J. 2018 Modifying and reacting to the environmental pH can drive bacterial interactions. PLoS Biol. 16 , e2004248. (10.1371/journal.pbio.2004248)29538378
52. Harrison S. 1991 Local extinction in metapopulation context: an empirical evaluation. Biol. J. Linn. Soc. 42 , 73-88. (10.1111/j.1095-8312.1991.tb00552.x)
53. Alldredge AI, Silver MW. 1988 Characteristics, dynamics and significance of marine snow. Prog. Oceanogr. 20 , 41-82. (10.1016/0079-6611(88)90053-5)
54. Gokhale S, Conwill A, Ranjan T, Gore J. 2018 Migration alters oscillatory dynamics and promotes survival in connected bacterial populations. Nat. Commun. 9 , 5273. (10.1038/s41467-018-07703-y)30531951
55. Solé R, Bascompte J, Valls J. 1992 Nonequilibrium dynamics in lattice ecosystems: chaotic stability and dissipative structures. Chaos 2 , 387-395. (10.1063/1.165881)12779988
56. Solé R, Gamarra JG. 1998 Chaos, dispersal and extinction in coupled ecosystems. J. Theor. Biol. 193 , 539-541. (10.1006/jtbi.1998.0716)9735280
57. Bell G, Gonzalez A. 2011 Adaptation and evolutionary rescue in metapopulations experiencing environmental deterioration. Science 332 , 1327-1330. (10.1126/science.1203105)21659606
58. Morozov A, Poggiale JC. 2012 From spatially explicit ecological models to mean-field dynamics: the state of the art and perspectives. Ecol. Complex. 10 , 1-11. (10.1016/j.ecocom.2012.04.001)
59. Maull V, Solé R. 2022 Network-level containment of single-species bioengineering. Phil. Trans. R. Soc. B 377 , 20210396. (10.1098/rstb.2021.0396)35757875
60. Crutchfield JP, Packard NH, Shaw RS, Farmer JD. 1980 Geometry from a time series. Phys. Rev. Lett. 45 , 712-716. (10.1103/PhysRevLett.45.712)
61. Schaffer WM. 1985 Order and chaos in ecological systems. Ecology 66 , 93-106. (10.2307/1941309)
62. Worden L. 2009 Notes from the Greenhouse World: a study in coevolution, planetary sustainability, and community structure. Ecol. Econ. 69 , 762-769. (10.1016/j.ecolecon.2009.06.017)
63. Williams HT, Lenton TM. 2007 The Flask model: emergence of nutrient-recycling microbial ecosystems and their disruption by environment altering ‘rebel’ organisms. Oikos 116 , 1087-1105. (10.1111/j.0030-1299.2007.15721.x)
64. Ratzke C, Denk J, Gore J. 2018 Ecological suicide in microbes. Nat. Ecol. Evol. 2 , 867-872. (10.1038/s41559-018-0535-1)29662223
65. Conde-Pueyo N, Vidiella B, Sardanyés J, Berdugo M, Maestre FT, De Lorenzo V, Solé R. 2020 Synthetic biology for terraformation lessons from mars, earth, and the microbiome. Life 10 , 14. (10.3390/life10020014)32050455
66. Saunders PT, Koeslag JH, Wessels JA. 1998 Integral rein control in physiology. J. Theor. Biol. 194 , 163-173. (10.1006/jtbi.1998.0746)9778431
67. Saunders PT., Koeslag JH, Wessels JA. 2000 Integral rein control in physiology II: a general model. J. Theor. Biol. 206 , 211-220. (10.1006/jtbi.2000.2118)10966758
68. Wood AJ, Ackland GJ, Lenton TM. 2006 Mutation of albedo and growth response produces oscillations in a spatial Daisyworld. J. Theor. Biol. 242 , 188-198. (10.1016/j.jtbi.2006.02.013)16581088
69. Wood AJ, Coe JB. 2007 A fitness based analysis of Daisyworld. J. Theor. Biol. 249 , 190-197. (10.1016/j.jtbi.2007.07.021)17854837
70. Nuño JC, De Vicente J, Olarrea J, López P, Lahoz-Beltra R. 2010 Evolutionary Daisyworld models: a new approach to studying complex adaptive systems. Ecol. Inform. 5 , 231-240. (10.1016/j.ecoinf.2010.03.003)
71. Punithan D, McKay RB. 2012 Evolutionary dynamics and ecosystems feedback in two dimensional daisyworld. Artif. Life 13 , 91-98.
72. Lenton TM. 1998 Gaia and natural selection. Nature 394 , 439-447. (10.1038/28792)9697767
73. Maull V, Pla Mauri J, Conde Pueyo N, Solé R. 2024 Code from: A synthetic microbial Daisyworld: planetary regulation in the test tube. Zenodo. (10.5281/zenodo.10475383)
