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 Figures
A (yeast cells) and
B (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 file
1.
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 Figure
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 Table
). 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. Figures
A and
B 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· 10
4 nM
−1(blue). As observed in Figure
, 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].
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 Figures
A and
B. 100 trajectories of the two species population for the different Gc concentrations are provided in Figures
C and
D.
As Figure
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 Figure
. The concentrations used in Figure
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.
Figures
C and
D 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 Table
.
| Table 3Standard 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
cmax 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. Figures
A (
S. cerevisiae) and
B (
E. coli) show average population sizes for different
cmax values.
As expected, an increase in cmax 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 cmax to find this minimum threshold under which the ecosystem goes extinct. According to our simulation results, the minimum cmax 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 Table
. Notably, the larger the reactor capacity, the higher is the deviation among the different population trajectories.
| Table 4Standard 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. Figure
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(Figures
A,
B), 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 Figures
C,
D 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 Figures
E and
F.
The bifurcation diagram describing our ecosystem’s oscillatory behavior is presented in Figure
. 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 Figure
B 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 Figure
A. 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.
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. Figure
shows a comparison between deterministic and stochastic simulations. For d = 0.50 h
−1(Figures
A and
B), 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. Figures
B and
C 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.
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.