Skip to main content

Stochastic simulations of a synthetic bacteria-yeast ecosystem

Abstract

Background

The field of synthetic biology has greatly evolved and numerous functions can now be implemented by artificially engineered cells carrying the appropriate genetic information. However, in order for the cells to robustly perform complex or multiple tasks, co-operation between them may be necessary. Therefore, various synthetic biological systems whose functionality requires cell-cell communication are being designed. These systems, microbial consortia, are composed of engineered cells and exhibit a wide range of behaviors. These include yeast cells whose growth is dependent on one another, or bacteria that kill or rescue each other, synchronize, behave as predator-prey ecosystems or invade cancer cells.

Results

In this paper, we study a synthetic ecosystem comprising of bacteria and yeast that communicate with and benefit from each other using small diffusible molecules. We explore the behavior of this heterogeneous microbial consortium, composed of Saccharomyces cerevisiae and Escherichia coli cells, using stochastic modeling. The stochastic model captures the relevant intra-cellular and inter-cellular interactions taking place in and between the eukaryotic and prokaryotic cells. Integration of well-characterized molecular regulatory elements into these two microbes allows for communication through quorum sensing. A gene controlling growth in yeast is induced by bacteria via chemical signals and vice versa. Interesting dynamics that are common in natural ecosystems, such as obligatory and facultative mutualism, extinction, commensalism and predator-prey like dynamics are observed. We investigate and report on the conditions under which the two species can successfully communicate and rescue each other.

Conclusions

This study explores the various behaviors exhibited by the cohabitation of engineered yeast and bacterial cells. The way that the model is built allows for studying the dynamics of any system consisting of two species communicating with one another via chemical signals. Therefore, key information acquired by our model may potentially drive the experimental design of various synthetic heterogeneous ecosystems.

Background

Advances in the field of synthetic biology have enabled the design of engineered cells performing human-defined functions at a single cell resolution[1, 2]. These functions include but are not limited to oscillators[3–5], bistable switches[6], bio-logical gates[7–9], riboregulators[10, 11] and molecular devices that control gene expression[12, 13]. Despite this progress, several limitations still exist. A major shortcoming is the decreased robustness and the limited potential complexity of single cell functions. Thus, attention has been shifted to synthetic systems based on communication between cells, rather than individual isolated cell functionality. Cooperation among cells is largely mediated by quorum sensing[14] and may be promising for the development of cell-systems that robustly perform complex tasks[15–17]. These tasks range from cells rescuing or killing one another[18–21] to cells synchronizing across a relatively long distance[22].

The potential advantage of microbial consortia compared to monocultures is two-fold. First, in contrast to monocultures, multicultures allow the different species to share the various required synthetic functions or the different steps of a synthetic function. This function sharing decreases the burden in the metabolism of the cells significantly. Second, the sharing of different functions, or steps, among different cells potentially renders microbial consortia more suited for fine-tuning of their artificial functionality[23].

It is now clear that mathematical models can accurately capture the behavior of synthetic systems comprising of either bacterial or yeast cell strains and allowing cell-to-cell communication[18–20, 22, 24–28]. You and his co-workers designed a synthetic bacterial ecosystem where cell-cell communication controls cell density by inducing a killer gene in the bacteria[19]. To mathematically investigate the dynamics of this system, they coupled their experiments with a simple deterministic model. Shou et al. designed a synthetic yeast system where cell growth was dependent on successful cell-cell communication[24]. To further explain their system behavior, they used a mathematical model comprised of algebraic equations. Basu and his co-workers designed a synthetic system, composed of bacteria, that forms different patterns of differentiation, such as rings and clovers, driven by cell-cell communication via N-Acyl homoserine lactone (AHL) signals[25]. In addition to experimentally designing this system, they used a deterministic mathematical model to explore the behavior of this system. Balagadde et al. designed a synthetic bacterial ecosystem where cell-cell communication enables cells to exhibit predator-prey dynamics by either killing or rescuing one another[18]. They initially developed a deterministic model to thoroughly study the dynamics of their synthetic ecosystem and then introduced a constant noise term to their model aiming to explore the influence of the stochasticity in their system.

Even though communication between different species using non-AHL signals has been demonstrated previously[29], no synthetic ecosystem has been developed that is composed of bacteria and yeast which communicate with and benefit from each other using AHL signals. Such a microbial consortium could exhibit interesting dynamics, such as oscillatory behavior, that stem from the substantial differences (e.g. different volume, growth rate, gene expression process) between prokaryotes and eukaryotes. Here, we investigate the behavior of such a synthetic heterogeneous community using stochastic modeling. To this end, we have modeled and simulated a synthetic consortium composed of Saccharomyces cerevisiae (S. cerevisiae) and Escherichia coli (E. coli) cells. This synthetic ecosystem was found to exhibit intriguing dynamic behavior that is commonly observed in natural ecosystems. Our model, capturing the behavior of this ecosystem, has been built in such a way that it can capture the dynamics of any system with two different species communicating with AHL signals. Thus, our model may drive the experimental design of artificial ecosystems with two different species (e.g. mammalian-yeast or mammalian-bacteria) which communicate with and regulate gene expression in one another.

Methods

Design of the synthetic ecosystem

In this study, we propose the design of a synthetic yeast-bacteria ecosystem that is based on diffusible chemical signals. Examples of these signals are the RhII/RhlR and LuxI/LuxR quorum sensing signals from Pseudomonas aeruiginosa and Vibrio fisheri quorum sensing systems, respectively, which are known for their sensitivity and the absence of signal cross-reactivity[30].

Each species exists in the presence of a molecule controlling growth, Gc. This molecule could be an antibiotic, such as Kanamycin, which is effective against both E. coli and S. cerevisiae[31]. Gc inhibits cell growth and therefore each species ultimately goes extinct. However, each species contains a resistance gene which counteracts the function of Gc and is controlled by the other species via diffusible molecules. Thus, when both species are present, they induce each other’s resistance gene through chemical signals, thereby rescuing one another. A schematic representation of the proposed ecosystem is illustrated in Figure1.

Figure 1
figure 1

Logic behind the synthetic yeast-bacteria ecosystem. S. cerevisiae cells produce AHL1 thereby activating resistance gene expression in E. coli and cell survival. Similarly, E. coli cells produce AHL2 that induces resistance gene expression in S. cerevisiae rescuing the latter.

More specifically, S. cerevisiae constitutively expresses a diffusible molecule, AHL1. AHL1 diffuses out of the S. cerevisiae cells, penetrating E. coli cells and binding to its cognate receptor, AHLR1. AHLR1 is constitutively produced in E. coli. The activated molecule in E. coli binds to the responsive element fused upstream of the Res promoter activating expression of Res. Subsequently, the resistance protein, Res, deactivates Gc in E. coli. Potential Res could be the Kanamycin resistance protein[31].

The second component of the feedback loop in E. coli functions in the same genetic fashion. It constitutively produces an autoinducer synthase, AHL2. Once AHL2 is produced, it diffuses out of the E. coli into S. cerevisiae, and is recognized by its cognate receptor, AHLR2, which is constitutively produced in S. cerevisiae as a fusion protein that allows it to be activated in eukaryotic cells. This activated molecule now binds to its responsive promoter and induces expression of the resistance gene, res. The resistance protein, in turn, represses the function of Gc in S. cerevisiae.

It is important to note that for the purposes of this study, we assume synthetic bacterial molecular components function in yeast. We hypothesize that their functionality may be retained when they are used in yeast. This is not an unreasonable hypothesis since the functionality of quorum sensing bacterial elements has been demonstrated experimentally in other higher organisms[32].

Here, we aim to computationally explore the behavior of a microbial consortium consisting of two different species, and how the differences of the two species affect its dynamics. The focus is therefore on the population dynamics. The functionality of such an ecosystem could in principle be achieved using any other molecular components with similar function.

Model description

As discussed in the previous section, numerous mathematical models that describe the behavior of synthetic ecosystems have been developed previously[18–20, 24–28]. The vast majority of these models are deterministic, ignoring the stochastic nature which is ubiquitous in biological systems[33–35]. Thus far, different methods have been described[36–42] and extensively applied to stochastically simulate the dynamics of biological systems in general and gene networks in particular[4, 8, 12, 13, 43–47].

In this study, we develop a stochastic model that accounts for the intrinsic and extrinsic noise and describes the dynamics of the synthetic bacteria-yeast ecosystem depicted in Figure1. The model takes into consideration the volume and the growth rate differences between E. coli and S. cerevisiae. In addition, it accounts for the gene expression dissimilarities between bacteria and yeast. Our model monitors the evolution of molecular species that usually exist in relatively high amounts allowing for the use of continuous stochastic models[38, 41]. Continuous stochastic computational approaches have also accurately described the experimental phenotype of synthetic cell communities[18]. We, therefore, employ chemical Langevin equations[41] to capture the evolution of the species participating in this synthetic ecosystem.

The model consists of 17 reactions whose dynamics are described using 9 Stochastic Differential Equations (see Additional file1). The equations were integrated in Matlab using the Euler Maruyama method[48]. The type of reactions as well as the kinetic parameters used were acquired from previously published studies involving experimental work. Our model is generic (i.e it may be used to capture the dynamics of various two-species ecosystems), but for the purposes of this study we assumed specific molecular components (and their associated kinetic parameters) that have been widely used in designing synthetic ecosystems. These components are presented in Table1. The current model may capture the behavior of any similar heterogeneous ecosystem by simply modifying the kinetic parameters according to the new system. The reaction network along with the kinetic parameters and the reaction rates capturing the behavior of our system is presented in Table2.

Table 1 Molecular components assumed in the model
Table 2 Reaction network capturing synthetic ecosystem’s behavior

The first two reactions describe the cell population growth. Consistent with previous mathematical models[18–20], and because the model refers to ecology, population growth follows logistic kinetics. Bacteria were considered to grow four times faster than yeast[49]; k1 was set four times smaller than k2. C max represents the carrying capacity of the bioreactor, i.e. the maximal population load that the bioreactor can sustain[51], and is set equal to 109 cells[18, 21]. Reactions 3 and 4 represent the cell death due to the presence of Gc (in our case Kanamycin). We assume a constant concentration of Gc (0.3 μM) as, according to the kinetic parameters used in our model, this concentration kills each single simulated cell colony when the two species are placed separately. Both the bacteria and yeast carry a resistance gene so the reaction rate is written such that the higher the amount of the resistance protein, the slower the cell death rate is. Similar reaction rates have been used previously to capture cell death due to killer proteins[18]. The correlation between Gc and the resistance protein is tuned through the parameter α. The parameter α was initially set equal to 5 · 104 molecules−1, due to the lack of literature values, and subsequently the sensitivity of the ecosystem’s behavior to changes in this parameter was investigated. Reactions 5 and 6 describe the production of the molecules responsible for the diffusible signals. AHL1 and AHLR2 are produced by S. cerevisiae whereas AHL2 and AHLR1 are produced by E. coli. The concentration of AHLR2 and AHLR1 is considered constant (0.5 μM) and equal to previously published values[26]. AHL1 and AHL2 production reactions are assumed to be first order, in accordance with previous studies[19, 20, 52]. The production rate of these diffusible molecules can vary significantly depending on the promoter strength of the associated genes. Using directed evolution, a wide range of quorum sensing production rates can be achieved[53]. The optimized behavior can be also achieved using computational approaches[54]. In our model, we initially adopted k4 and k5 from[18] and subsequently increased their values since our system required very long time to reach steady state under these conditions. Reaction 7 captures the binding of AHL2 to AHLR2 in S. cerevisiae. This reaction is considered a fourth order reaction (this reaction accounts for the volume of S. cerevisiae cells) since it has been demonstrated that a fourth order reaction can capture the experimental phenotypes well[26]. Resistance protein (Res) production is calculated using Hill type kinetics, in accordance with experimental observations[26], and is shown in reaction 8. This reaction also accounts for gene expression differences between eukaryotes and prokaryotes. In contrast to prokaryotes, eukaryotic transcription requires many transcription factors to be recruited before its initiation. Moreover, the translation process in prokaryotes is faster than in eukaryotes[55]. These two factors introduce a delay in eukaryotic gene expression rendering it slow compared to the prokaryotic gene expression. In our model, we represent this delay using reaction 8. In fact, we assume that a complex (preRes) must first be formed before Res production can take place. A similar approach has been used previously to capture transcription in yeast[56]. After performing a set of simulations, we set k8 equal to 5 h−1 since this value was found to cause a delay in our ecosystem compared to a model lacking this intermediate reaction (data not shown). The actual process of protein production is captured by reaction 9. Similarly to reaction 7, reaction 10 captures the binding of AHL1 to AHLR1 in E. coli (this reaction accounts for the E. coli cell volume). Reaction 11 is used to describe E. coli gene expression. Note that in this case there is no reaction describing a delay in gene expression. Finally, reactions 12-17 represent the degradation of the species participating in this network and they are all considered first-order.

Results and discussion

Testing synthetic ecosystem’s functionality

Initially, we explored whether S. cerevisiae cells can withstand Gc in the absence of E. coli cells and vice versa. We simulated the behavior of 50,000 yeast cells and 50,000 bacterial cells in the absence of Gc. Our simulations indicated that both S. cerevisiae and E. coli grow normally (data not shown). However, when each of the two populations is placed in a simulated bioreactor separately, in the presence of Gc, neither population is able to survive (data not shown).

Subsequently, we simulated the behavior of bacteria and yeast when both are present to test whether communication and cooperation between these two species can be successfully achieved. We simulated stochastically 100 different population colonies containing 50,000 E. coli and 50,000 S. cerevisiae cells. The results are depicted in Figures2A (yeast cells) and2B (bacterial cells). The different lines correspond to the cell population size in different trajectories. The evolution of all the species is provided in the Additional file1.

Figure 2
figure 2

Behavior of coexisting S. cerevisiae and E. coli cells. When the two species are placed together, obligatory mutualism is observed, i.e they benefit from each other and survive from Gc. The inset represents part of Figure2B and shows the fluctuations of E. coli population size.

In both cases, the two different species exploit communication with one another for successful survival in the presence of Gc. Our simulations demonstrate that yeast can successfully induce the expression of the resistance gene found in bacteria and vice versa. This is a common characteristic of ecosystems called obligatory mutualism. In other words, S. cerevisiae and E. coli cells are not able to survive separately but they are able to grow in concert. As expected, we observe that the number of E. coli cells is always higher than the number of S. cerevisiae cells. As discussed before, the reason for this is the high growth rate of bacteria relative to yeast.

Variation regarding the number of cells is observed when different colonies are simulated and is attributed to the stochasticity underlying biological functions. The average S. cerevisiae population (calculated over 100 trajectories) is around 4.90· 105 cells and the standard deviation is equal to 5.77· 104 cells. The mean E. coli population is about 9.99·108 cells and the standard deviation is approximately 6.68·104 cells. Note that the total number of cells cannot exceed 109. Both bacteria and yeast require approximately 16 hours to reach steady state. In every single bioreactor, neither bacteria nor yeast die from the presence of Gc. This demonstrates that communication can take place between S. cerevisiae and E. coli allowing for the survival of the two species.

Even though Figure2 establishes cell communication and obligatory mutualism between E. coli and S. cerevisiae cells, this refers only to the case described by this set of parameters. Thus, in order to investigate which parameters promote successful communication and cooperation between E. coli and S. cerevisiae cells, and to explore the dynamics of different parameter sets, a sensitivity analysis was performed. To implement this, we systematically modified different parameters within reasonable ranges and monitored the dynamics of the system. In what follows, we present the evolution of the average S. cerevisiae and E. coli population over 100 trajectories. In some cases, we further provide all the 100 trajectories with variation in the values of key parameters examined in our analysis.

Ecosystem’s sensitivity to parameter α

As discussed in the previous section, α represents a correlation between the molecule controlling growth and the resistance protein. More specifically, the larger the α, the smaller the amount of Res required for cells to survive from Gc (see Table2). Since this parameter is of high importance for our model, and because it was the only parameter not acquired from previously published models, we explored the influence of α on the system’s behavior. To this end, we performed multiple computational experiments modifying α and investigating our ecosystem’s dynamic behavior. Our simulation results showed that when α is larger than 100 nM−1, the total system’s behavior does not change appreciably (data not shown). For values of α smaller than 0.07 nM−1, the ecosystem is driven to extinction (data not shown). Importantly, our simulations’ data demonstrated that for values of α in the range of 0.07 nM−1 to 100 nM−1, the dynamics of the system, and specifically the time the system needs to reach steady state, becomes remarkably slow. Figures3A and3B show the mean values, along with the standard deviation, of 100 trajectories from the stochastic simulations for α equal to 25 nM−1 (red), 75 nM−1(green) and 5· 104 nM−1(blue). As observed in Figure3, when the value of α is lower than 100, the system cannot reach steady state even after 10,000 hours. To our knowledge, no synthetic ecosystem exists that reaches steady state after such a long time suggesting that in order for this system to be realistic, the value of α in our model should be higher than 100 nM−1. This is confirmed by the fact that our system reaches steady state approximately as fast as previously published synthetic ecosystems systems did[18–20].

Figure 3
figure 3

Average values and standard deviation of S. cerevisiae (A) and E. coli (B) population for different values of α. Mean values and standard deviation (grey shade) of 100 trajectories of S. cerevisiae (A) and E. coli (B) population size for different values of the parameter α.

Importance of Gc concentration

Previous studies describing similar synthetic ecosystems have demonstrated the importance of the concentration of the molecule controlling growth on the system’s dynamics[20]. Guided by this, we conducted a set of simulations where we modified Gc’s concentration. We monitored the dynamics of the system for three different Gc concentrations. The average population values (over 100 trajectories) for each concentration are shown in Figures4A and4B. 100 trajectories of the two species population for the different Gc concentrations are provided in Figures4C and4D.

Figure 4
figure 4

Average values and single trajectories of S. cerevisiae and E. coli population for different Gc concentrations. Average (over 100 trajectories) values (A,B) and 100 single trajectories (C,D) of S. cerevisiae and E. coli population size for Gc concentration equal to 60 nM (red), 250 nM (black) and 10 μM(blue). The synthetic ecosystem adopts different behaviors, that are commonly observed in natural ecosystems, in response to different Gc concentrations.

As Figure4 indicates, an increase on Gc’s concentration from 60 nM to 250 nM is followed by a decreased yeast population and an increased bacterial population. In other words, upon increasing Gc concentration in the bioreactor, E. coli cells benefit whereas S. cerevisiae cells are harmed. Based on the way our model was built, this is likely ascribed to the fact that yeast grow much slower than bacteria and can therefore resist only low Gc concentrations. As the antibiotic concentration increases, yeast die faster than bacteria and the latter, even though they grow slower than they would in the absence of Gc, take advantage of the higher nutrient levels in the bioreactor. This is an interesting characteristic and could be used as a means for controlling the bioreactor’s population, obviating the need of adding or removing cells. However, when Gc concentration is significantly high (e.g. 10 μ M), the average value (of the 100 trajectories) of both populations decreases dramatically as many single trajectories reach zero.

Further analysis of the system’s behavior indicated that changing the Gc concentration leads to an intriguing behavior commonly exhibited by natural ecosystems. More specifically, when Gc levels are low, each species can survive even in the absence of the other species. In particular, bacterial cells can withstand up to 250 nM Gc. On the other hand, yeast cells cannot survive even these Gc levels and they can only withstand Gc concentrations lower than 60 nM. Having said this, the behavior of the system for Gc levels up to 60 nM is analogous to facultative mutualism, i.e. both species benefit from but are not dependent on each other. However, when Gc’s level lies between 60 nM and 250 nM, the behavior of the system is similar to commensalism for bacteria, i.e. bacteria can survive without yeast but yeast are not able to survive without bacteria.

The lethal Gc concentration for cultures with both cell types present is 20 μ M. Based on this, we conclude that when Gc’s levels are between 250 nM and 20 μ M, the behavior of the system is homologous to obligatory mutualism as both species are completely dependent on each other and unable to survive individually. Finally, for Gc levels higher than 20 μ M, we observe ecosystem’s extinction. Such behaviors have been observed previously in similar synthetic bacterial ecosystems[20] and are shown in Figure4. The concentrations used in Figure4 represent the boundaries between different system’s behavior (note that instead of 20 μ M Gc, which is the boundary between obligatory mutualism and extinction, we considered 10 μ M Gc). The population dynamic behavior for Gc concentrations between the ones used here lies in the area between these lines.

Figures4C and4D demonstrate deviation among the different cell density trajectories. Note that this deviation could not be captured using deterministic simulations. The standard deviation (calculated over the 100 trajectories) of the population size at 50 hours and for 60 nM, 250 nM, and 20 μ M Gc is shown in Table3.

Table 3 Standard deviation of population size at steady state for different Gc concentrations

Ecosystem’s sensitivity to various carrying capacities and initial cell densities

We then explored the influence of c max and the initial cell population on the synthetic ecosystem’s behavior. As discussed above, the carrying capacity is the maximum number of (bacterial and yeast) cells that can exist in the bioreactor[51]. Here, we only show average values of the 100 trials since the single trajectories exhibit similar behavior as in the previous cases. Figures5A (S. cerevisiae) and5B (E. coli) show average population sizes for different c max values.

Figure 5
figure 5

Average (over 100 trajectories) S. cerevisiae (A) and E. coli (B) population size for different c max . Average (over 100 trajectories) values of S. cerevisiae (A) and E. coli (B) population size for different reactor capacities. The higher the reactor capacity the higher the steady state population density of the two species is.

As expected, an increase in c max causes an increase on both yeast and bacterial steady state populations as the nutrients in the culture suffice for more cells. Thus, both species grow faster and consequently survive in the presence of Gc. It is important to note that a minimum amount of nutrients must exist in the bioreactor for the cells to grow and survive. Thus, we ran simulations decreasing c max to find this minimum threshold under which the ecosystem goes extinct. According to our simulation results, the minimum c max in order for all the trajectories to end up in non-zero steady states (over a period of 3,000 hours) is equal to approximately 2 · 105 cells (data not shown). Thus, the model suggests that our synthetic ecosystem is fully functional only for reactor capacities equal to or higher than 2 · 105 cells.

As in the previous cases, deviation among the different population trajectories was observed. The steady state standard deviation of the population size for various reactor capacities is provided in Table4. Notably, the larger the reactor capacity, the higher is the deviation among the different population trajectories.

Table 4 Standard deviation of population size at steady state for different reactor capacities

We further explored the minimum initial number of total cells required in order for the two species population to cooperate favorably and survive. To do so, we performed different simulations starting with equal E. coli and S. cerevisiae populations and monitoring the system’s dynamics for 1,000 hours. Our results showed that for equal initial populations of the two species, the minimum number of S. cerevisiae and E. coli cells in the reactor should be approximately equal to 15 cells for the ecosystem to survive with Gc. Moreover, when the initial E. coli population is 50,000 cells, the minimum S. cerevisiae initial population required in order for the system to avoid extinction is 14 cells. On the other hand, when the initial S. cerevisiae population is 50,000 cells, the required E. coli initial population is 4 cells. This difference is ascribed to the fact that bacteria grow predominantly fast thereby quickly helping yeast to survive and therefore only 4 yeast cells are initially required to make the ecosystem functional. However, yeast grow and consequently rescue E. coli with a slower rate and therefore larger E. coli population is initially required for the ecosystem to function.

Effects of E. coli cell death rates on the ecosystem’s dynamics

It is clear from the aforementioned analysis that in most cases bacterial cell populations dominate yeast cell populations because of their high growth rate. We therefore introduced a bacteria degradation term in our network to enhance the competition between the population of the two species. We only considered E. coli degradation as bacteria grow significantly faster than yeast. This degradation could be achieved experimentally as bacteria can be engineered to stimulate their lysis in response to a human-defined signal. More specifically, introducing holin and lysozome genes that are activated via AHL signals, allows for controlling cell membrane destruction and consequently cell death[57].

Initially, we performed our analysis under the assumption that the deterministic term dominates the stochastic term, i.e. the intrinsic noise of the system is negligible. The results presented in what follows were therefore produced based only on the deterministic part of the equations 1-9. A similar approach has been used previously to explore the oscillatory behavior of a synthetic ecosystem[18].

As expected, high degradation rates cause bacterial cell death followed by yeast wash out due to obligatory mutualism (data not shown). In contrast, low degradation rates allow yeast domination, as bacterial populations quickly decreases due to both Gc and degradation, thereby allowing an increase in yeast population (data not shown).

Importantly, and as observed in other synthetic ecosystems composed of species with different growth rates[18], there is a range of bacterial degradation rate where S. cerevisiae and E. coli population exhibit sustained oscillations. These oscillations originate from the antagonism between the two species population and demonstrate a predator-prey like relationship between S. cerevisiae and E. coli cells. In particular, when the E. coli degradation rate, d, lies between 0.30 h−1 and 0.72 h−1, sustained oscillations are observed. For d smaller than 0.30 h−1, damped oscillations are exhibited. Finally, for d larger than 0.72 h−1, the two species population goes to zero. In other words, the ecosystem becomes extinct, since this high degradation rate results in bacterial death which in turn leads to yeast extinction because cell-cell communication cannot take place favorably anymore. Figure6 shows the behavior of the two species population for degradation parameters that lie in the aforementioned ranges. When d is equal to 0.25 h−1(Figures6A,6B), we observe damped oscillations that end up on a stable steady state. However, from d = 0.30 h−1 to d = 0.72 h−1, sustained oscillations whose amplitude scales with the degradation rate are observed. This trend is provided in Figures6C,6D where d = 0.50 h−1. Finally, when d is larger than 0.72 h−1 the system is driven to extinction, as depicted in Figures6E and6F.

Figure 6
figure 6

S. cerevisiae (A,C,E) and E. coli (B,D,F) population dynamics for different E. coli degradation rates. S. cerevisiae (A,C,E) and E. coli (B,D,F) population dynamics for E. coli degradation rate equal to 0.25, 0.50 and 0.75 h−1. For d = 0.25 h−1(A,B) the ecosystem exhibits damped oscillations. For d = 0.50 h−1(C,D) the population of the two species oscillates with sustained oscillations whereas for d = 0.75 h−1(E,F) goes to zero.

The bifurcation diagram describing our ecosystem’s oscillatory behavior is presented in Figure7. A Hopf point, where sustained oscillations of the two species population initiate, is observed for d approximately equal to 0.30 h−1(red). The Hopf point was further confirmed by eigenvalue analysis. The lines following the Hopf point correspond to the oscillation amplitude, as calculated from the transient analysis. Please note that for the sake of clarity, in Figure7B we show the upper limit of oscillations only for d up to 0.41 h−1. The complete bifurcation diagram is provided in the inset. The period of the oscillations was calculated using the FFT (Fast Fourier Transform) function in Matlab and is depicted as inset in Figure7A. As evident, the period of the oscillations scales with the E. coli degradation rate. This is an intriguing observation which suggests that the E. coli degradation rate could be used to control the period of our oscillatory ecosystem.

Figure 7
figure 7

Bifurcation diagram of the S. cerevisiae (A) and E. coli (B) population versus the E. coli degradation rate. Bifurcation diagram of the S. cerevisiae (A) and E. coli (B) population versus the degradation rate of E. coli cells. For the sake of clarity, Figure B shows only part of the bifurcation diagram whereas the complete bifurcation diagram is illustrated in the inset. The period of oscillation of S. cerevisiae and E. coli cells for different E. coli degradation rates is the same and presented as inset in Figure7A.

It should be stressed that including the stochastic terms in our simulations, leads to the ecosystem’s extinction. This has been observed before[18] and is caused by the fact that during the oscillations, the bacterial population reaches small values and therefore noise terms destroy the sustained oscillations by driving bacterial population to zero and consequently the ecosystem to extinction (since cooperation cannot occur). In fact, the smaller the noise amplitude, the higher the probability for the system to circumvent extinction and exhibit sustained oscillations. Figure8 shows a comparison between deterministic and stochastic simulations. For d = 0.50 h−1(Figures8A and8B), deterministic solution (black) provides sustained oscillations that end up in a steady state. Stochastic simulations (red) are consistent with the deterministic ones, i.e. demonstrate oscillations, but only for a small period of time and subsequently all the trajectories reach zero. Motivated by this observation, we performed several simulations (for d = 0.50 h−1) where we systematically decreased the noise terms amplitude. Our simulations demonstrated that when the noise terms are 1.25% or less of the current values, the stochastic behavior matches the deterministic one, i.e. the ecosystem population exhibits sustained oscillations. When the noise terms are between 1.30% and 100% of the current values, there are always stochastic trajectories that reach zero over a period of 1,000 hours. Figures8B and8C show three population density trajectories of the stochastic simulation (green, red, blue) compared with the deterministic simulation (black) when the noise terms are reduced to 1.25%. As evident, the ecosystem’s behavior provided by the two approaches is consistent and stochastic trajectories exhibit continuous oscillations.

Figure 8
figure 8

Population density of S. cerevisiae and E. coli for d = 0.5 h −1calculated using stochastic and deterministic simulations. A,B: Population size of S. cerevisiae (A) and E. coli (B) for d = 0.50 h−1calculated using stochastic (red) and deterministic simulations (black). C,D: Population size of S. cerevisiae (C) and E. coli (D) for d = 0.50 h−1calculated with deterministic (black) and stochastic (red, green, blue) simulations with 1.25% of the current intrinsic noise terms.

Overall, our simulations suggest that high amplitude intrinsic noise damages the ecosystem’s oscillatory behavior. On the other hand, less noisy environments stimulate the sustained oscillation of the two species population.

Conclusions

We presented the in silico design of the first synthetic bacterial-yeast ecosystem where communication between cells is achieved using AHL signals. The model, while developed to accurately depict these interactions, can be adapted to characterize any cell-to-cell communication and population dynamics mediated by diffusible chemical signaling.

We showed that when the two species coexist, they overcome Gc’s toxicity by inducing each other’s resistance gene via small molecule signalling and therefore survive. Our simulations suggest that the minimum reactor capacity required for this ecosystem to evolve is 2 · 105cells. By varying the Gc concentration, the ecosystem adopts different behaviors including obligatory and facultative mutualism, commensalism and extinction. Adding an E. coli degradation reaction, which can be experimentally realized by engineering bacteria to induce lysis, can drive the population of the two species to predator-prey like dynamics, i.e. sustained oscillations. These oscillations can, however, be destroyed in noisy environments. Overall, we demonstrated that such kind of heterogeneous synthetic ecosystems could exhibit interesting dynamics.

As demonstrated here and in different studies[18], the development of synthetic microbial consortia using species with different characteristics (e.g. different growth rate or volume) yields systems with intriguing dynamics, such as oscillations. These systems could have various potential applications such as the delivery of two different drugs in dissimilar time intervals[23].

Our mathematical model may potentially drive the experimental design of microbial consortia with a heterogeneous population. This and similar mathematical models can further be used to predict interspecies bioreactor dynamics under numerous conditions, with differing chemical signals, and employing various population control mechanisms. Engineered interspecies system have substantial implications for complex chemical synthesis as well as future biorefinery design and optimization. Thus, the dynamics analysis presented herein may be used as the basis for the in vivo design of such promising synthetic ecosystems.

Authors contributions

KB developed the model, carried out the simulations and wrote the manuscript. DB participated in the design of the study. CSD conceived of the study and participated in its design. YK conceived of the study, participated in its design and coordination and helped to write the manuscript. All authors read and approved the final manuscript.

References

  1. Andrianantoandro E, Basu S, Karig D, Weiss R: Synthetic biology: new engineering rules for an emerging discipline. Mol Syst Biol 2006, 2: 2006-0028.

    Article  Google Scholar 

  2. Khalil A, Collins J: Synthetic biology: applications come of age. Nat Rev Genet 2010,11(5):367.

    Article  CAS  Google Scholar 

  3. Elowitz M, Leibler S: A synthetic oscillatory network of transcriptional regulators. J Biol Chem 1999, 274: 6074-6079. 10.1074/jbc.274.10.6074

    Article  Google Scholar 

  4. Tuttle L, Salis H, Tomshine J, Kaznessis Y: Model-driven designs of an oscillating gene network. Biophys J 2005,89(6):3873-3883. 10.1529/biophysj.105.064204

    Article  CAS  Google Scholar 

  5. Stricker J, Cookson S, Bennett M, Mather W, Tsimring L, Hasty J: A fast, robust and tunable synthetic gene oscillator. Nature 2008,456(7221):516-519. 10.1038/nature07389

    Article  CAS  Google Scholar 

  6. Gardner T, Cantor C, Collins J: Construction of a genetic toggle switch inEscherichia coli. Nature 2000, 403: 339-342. 10.1038/35002131

    Article  CAS  Google Scholar 

  7. Win M, Smolke C: Higher-order cellular information processing with synthetic RNA devices. Science 2008,322(5900):456. 10.1126/science.1160311

    Article  CAS  Google Scholar 

  8. Ramalingam K, Tomshine J, Maynard J, Kaznessis Y: Forward engineering of synthetic bio-logical AND gates. Biochem Eng J 2009,47(1-3):38-47. 10.1016/j.bej.2009.06.014

    Article  Google Scholar 

  9. Anderson J, Voigt C, Arkin A: Environmental signal integration by a modular AND gate. Mol Syst Biol 2007, 3: 133.

    Article  Google Scholar 

  10. Bayer T, Smolke C: Programmable ligand-controlled riboregulators of eukaryotic gene expression. Nat Biotechnol 2005,23(3):337-343. 10.1038/nbt1069

    Article  CAS  Google Scholar 

  11. Win M, Smolke C: A modular and extensible RNA-based gene-regulatory platform for engineering cellular function. Proc Nat Acad Sci 2007,104(36):14283. 10.1073/pnas.0703961104

    Article  CAS  Google Scholar 

  12. Volzing K, Biliouris K, Kaznessis Y: ProTeOn and ProTeOff, new protein devices that inducibly activate bacterial gene expression. ACS Chemical Biology 2011,6(10):1107-1116. 10.1021/cb200168y

    Article  CAS  Google Scholar 

  13. Sotiropoulos V, Kaznessis Y: Synthetic tetracycline-inducible regulatory networks: computer-aided design of dynamic phenotypes. BMC Syst Biol 2007, 1: 7. 10.1186/1752-0509-1-7

    Article  Google Scholar 

  14. Miller M, Bassler B: Quorum sensing in bacteria. Annu Rev Microbiol 2001, 55: 165-199. 10.1146/annurev.micro.55.1.165

    Article  CAS  Google Scholar 

  15. Choudhary S, Schmidt-Dannert C: Applications of quorum sensing in biotechnology. Appl Microbiol Biotechnol 2010,86(5):1267-1279. 10.1007/s00253-010-2521-7

    Article  CAS  Google Scholar 

  16. Pai A, Tanouchi Y, Collins C, You L: Engineering multicellular systems by cell-cell communication. Curr Opin Biotechnol 2009,20(4):461-470. 10.1016/j.copbio.2009.08.006

    Article  CAS  Google Scholar 

  17. Tanouchi Y Smith RP You L: Engineering microbial systems to explore ecological and evolutionary dynamics. Curr Opin Biotechnol 2012,23(5):791-7. 10.1016/j.copbio.2012.01.006

    Article  Google Scholar 

  18. Balagaddé F, Song H, Ozaki J, Collins C, Barnet M, Arnold F, Quake S, You L: A synthetic Escherichia coli predator–prey ecosystem. Mol Syst Biol 2008, 4: 187.

    Article  Google Scholar 

  19. You L, Cox R, Weiss R, Arnold F: Programmed population control by cell-cell communication and regulated killing. Nature 2004,428(6985):868-871. 10.1038/nature02491

    Article  CAS  Google Scholar 

  20. Hu B, Du J, Zou R, Yuan Y: An environment-sensitive synthetic microbial ecosystem. PloS one 2010,5(5):e10619. 10.1371/journal.pone.0010619

    Article  Google Scholar 

  21. Balagadde F, You L, Hansen C, Arnold F, Quake S: Long-term monitoring of bacteria undergoing programmed population control in a microchemostat. Science 2005,309(5731):137. 10.1126/science.1109173

    Article  CAS  Google Scholar 

  22. Prindle A, Samayoa P, Razinkov I, Danino T, Tsimring L, Hasty J: A sensing array of radically coupled genetic/biopixels/’. Nature 2011, 481: 39-44. 10.1038/nature10722

    Article  Google Scholar 

  23. Brenner K, You L, Arnold F: Engineering microbial consortia: a new frontier in synthetic biology. Trends Biotechnol 2008,26(9):483-489. 10.1016/j.tibtech.2008.05.004

    Article  CAS  Google Scholar 

  24. Shou W, Ram S, Vilar J: Synthetic cooperation in engineered yeast populations. Proc Nat Acad Sci 2007,104(6):1877. 10.1073/pnas.0610575104

    Article  CAS  Google Scholar 

  25. Basu S, Mehreja R, Thiberge S, Chen M, Weiss R: Spatiotemporal control of gene expression with pulse-generating networks. Proc Nat Acad Sci USA 2004,101(17):6355. 10.1073/pnas.0307571101

    Article  CAS  Google Scholar 

  26. Basu S, Gerchman Y, Collins C, Arnold F, Weiss R: A synthetic multicellular system for programmed pattern formation. Nature 2005,434(7037):1130-1134. 10.1038/nature03461

    Article  CAS  Google Scholar 

  27. Brenner K, Karig D, Weiss R, Arnold F: Engineered bidirectional communication mediates a consensus in a microbial biofilm consortium. Proc Nat Acad Sci 2007,104(44):17300. 10.1073/pnas.0704256104

    Article  CAS  Google Scholar 

  28. Anesiadis N, Cluett W, Mahadevan R: Dynamic metabolic engineering for increasing bioprocess productivity. Metab Eng 2008,10(5):255-266. 10.1016/j.ymben.2008.06.004

    Article  CAS  Google Scholar 

  29. Weber W, Daoud-El Baba M, Fussenegger M: Synthetic ecosystems based on airborne inter-and intrakingdom communication. Proc Nat Acad Sci 2007,104(25):10435. 10.1073/pnas.0701382104

    Article  CAS  Google Scholar 

  30. Smith C, Song H, You L: Signal discrimination by differential regulation of protein stability in quorum sensing. J Mol Biol 2008,382(5):1290-1297. 10.1016/j.jmb.2008.08.009

    Article  CAS  Google Scholar 

  31. Webster T, Dickson R: Direct selection of Saccharomyces cerevisiae resistant to the antibiotic G418 following transformation with a DNA vector carrying the kanamycin-resistance gene of Tn903. Gene 1983,26(2-3):243-252. 10.1016/0378-1119(83)90194-4

    Article  CAS  Google Scholar 

  32. You Y, Marella H, Zentella R, Zhou Y, Ulmasov T, Ho T, Quatrano R: Use of bacterial quorum-sensing components to regulate gene expression in plants. Plant Physiol 2006,140(4):1205. 10.1104/pp.105.074666

    Article  CAS  Google Scholar 

  33. McAdams H, Arkin A: Stochastic mechanisms in gene expression. Proc Nat Acad Sci 1997,94(3):814. 10.1073/pnas.94.3.814

    Article  CAS  Google Scholar 

  34. Elowitz M, Levine A, Siggia E, Swain P: Stochastic gene expression in a single cell. Science 2002,297(5584):1183. 10.1126/science.1070919

    Article  CAS  Google Scholar 

  35. Kærnm M, Elston T, Blake W, Collins J: Stochasticity in gene expression. Nat Rev Genet 2005, 6: 451-464. 10.1038/nrg1615

    Article  Google Scholar 

  36. Gillespie D: A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. J Comput Phys 1976,22(4):403-434. 10.1016/0021-9991(76)90041-3

    Article  CAS  Google Scholar 

  37. Salis H, Kaznessis Y: An equation-free probabilistic steady-state approximation: dynamic application to the stochastic simulation of biochemical reaction networks. J Chem Phys 2005, 123: 214106. 10.1063/1.2131050

    Article  Google Scholar 

  38. Salis H, Kaznessis Y: Accurate hybrid stochastic simulation of a system of coupled chemical or biochemical reactions. J Chem Phys 2005, 122: 054103. 10.1063/1.1835951

    Article  Google Scholar 

  39. Kaznessis Y: Models for synthetic biology. BMC Syst Biol 2007, 1: 47. 10.1186/1752-0509-1-47

    Article  Google Scholar 

  40. Gibson M, Bruck J: Efficient exact stochastic simulation of chemical systems with many species and many channels. J Phys Chem A 2000,104(9):1876-1889. 10.1021/jp993732q

    Article  CAS  Google Scholar 

  41. Gillespie D: The chemical Langevin equation. J Chem Phys 2000, 113: 297. 10.1063/1.481811

    Article  CAS  Google Scholar 

  42. Charlebois DA, Intosalmi J, Fraser D, Kærn M: An Algorithm for the Stochastic Simulation of Gene Expression and Heterogeneous Population Dynamics. Commun Comput Phys 2011, 9: 89-112.

    Google Scholar 

  43. Biliouris K, Daoutidis P, Kaznessis Y: Stochastic simulations of the tetracycline operon. BMC Syst Biol 2011, 5: 9. 10.1186/1752-0509-5-9

    Article  Google Scholar 

  44. Salis H, Kaznessis Y: Numerical simulation of stochastic gene circuits. Comput Chem Eng 2005,29(3):577-588. 10.1016/j.compchemeng.2004.08.017

    Article  CAS  Google Scholar 

  45. De Jong H: Modeling and simulation of genetic regulatory systems: a literature review. J Comput Biol 2002, 9: 67-103. 10.1089/10665270252833208

    Article  CAS  Google Scholar 

  46. Hasty J, McMillen D, Isaacs F, Collins J, et al: Computational studies of gene regulatory networks: in numero molecular biology. Nat Rev Genet 2001,2(4):268-279.

    Article  CAS  Google Scholar 

  47. Charlebois D, Perkins T, Kaern M: Stochastic gene expression and the processing and propagation of noisy signals in genetic networks. In Information Processing and Biological Systems. Volume 11. Edited by: Niiranen S, Ribeiro A. Springer-Verlag, Berlin Heidelberg; 2011:89-112.

    Chapter  Google Scholar 

  48. Kloeden P, Platen E, Schurz H: Stochastic differential equations. Numer Solution SDE Through Comput Exp 1994, 1: 63-90.

    Article  Google Scholar 

  49. Milo R, Jorgensen P, Moran U, Weber G, Springer M: BioNumbers the database of key numbers in molecular and cell biology. Nucleic Acids Res 2010,38(suppl 1):D750.

    Article  CAS  Google Scholar 

  50. Sundararaj S, Guo A, Habibi-Nazhad B, Rouani M, Stothard P, Ellison M, Wishart D: The CyberCell Database (CCDB): a comprehensive, self-updating, relational database to coordinate and facilitate in silico modeling of Escherichia coli. Nucleic Acids Res 2004,32(suppl 1):D293.

    Article  CAS  Google Scholar 

  51. Hui C: Carrying capacity, population equilibrium, and environment’s maximal load. Ecol Modell 2006,192(1-2):317-320. 10.1016/j.ecolmodel.2005.07.001

    Article  Google Scholar 

  52. Kambam PK SL Henson MA: Design and mathematical modelling of a synthetic symbiotic ecosystem. IET Syst Biol 2008, 2: 33-38. 10.1049/iet-syb:20070011

    Article  Google Scholar 

  53. Collins C, Arnold F, Leadbetter J: Directed evolution of Vibrio fischeri LuxR for increased sensitivity to a broad spectrum of acyl-homoserine lactones. Mol Microbiol 2005,55(3):712-723.

    Article  CAS  Google Scholar 

  54. Tomshine J, Kaznessis Y: Optimization of a stochastically simulated gene network model via simulated annealing. Biophys J 2006,91(9):3196-3205. 10.1529/biophysj.106.083485

    Article  CAS  Google Scholar 

  55. Brooker RJ: Genetics: Analysis and Principles. Janice Roerig-Blong, McGraw Hill; 2009.

    Google Scholar 

  56. Blake W, Kærn M, Cantor C, Collins J: Noise in eukaryotic gene expression. Nature 2003,422(6932):633-637. 10.1038/nature01546

    Article  CAS  Google Scholar 

  57. Lorenzo P, Susanna Z, Manuel L, Maria C, Paolo M: Characterization of a synthetic bacterial self-destruction device for programmed cell death and for recombinant proteins release. J Biol Eng 2011, 5: 8. 10.1186/1754-1611-5-8

    Article  Google Scholar 

Download references

Acknowledgements

This work was supported by a grant from the University of Minnesota Initiative for Renewable Energy and the Environment (Rl-001-11), a grant from the National Institutes of Health (American Recovery and Reinvestment Act grant GM086865) and a grant from the National Science Foundation (CBET-0644792). Partial funding for this work was also provided through a grant from the Synthetic Ecology Program of the Biotechnology Institute and from the Presidential Initiative on Biocatalysis at the University of Minnesota.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Yiannis N Kaznessis.

Additional information

Competing interests

The authors declare that they have no competing interests.

Electronic supplementary material

12918_2011_953_MOESM1_ESM.pdf

Additional file 1:Stochastic differential equations used to simulate the behavior of the synthetic ecosystem. This file contains the equations and the species used to stochastically simulate the behavior of the synthetic ecosystem. It also includes 100 trajectories of the evolution of all the species when S. cerevisiae and E. coli coexist (In support of Figure2). (PDF 276 KB)

Authors’ original submitted files for images

Rights and permissions

Open Access This article is published under license to BioMed Central Ltd. This is an Open Access article is distributed under the terms of the Creative Commons Attribution License ( https://creativecommons.org/licenses/by/2.0 ), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and permissions

About this article

Cite this article

Biliouris, K., Babson, D., Schmidt-Dannert, C. et al. Stochastic simulations of a synthetic bacteria-yeast ecosystem. BMC Syst Biol 6, 58 (2012). https://doi.org/10.1186/1752-0509-6-58

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/1752-0509-6-58

Keywords