As shown in Box 1A
, two E. coli
populations (‘predator' and ‘prey') communicate and regulate each other's density through quorum sensing (QS), a mechanism by which many bacteria broadcast and respond to changes in their local population density (Camilli and Bassler, 2006
). Two QS modules, LuxI/LuxR from Vibrio fischeri
and LasI/LasR from Pseudomonas aeruginosa
, are used to enable two-way communication between the predator and prey populations. When the prey density is low, the predator cells die due to constitutive expression of a suicide gene (ccdB
). In prey cells an acyl-homoserine lactone (AHL), 3OC6HSL, is synthesized by LuxI. As the prey density increases, 3OC6HSL accumulates in the culture and, at sufficiently high concentrations, it is bound by the transcriptional regulator LuxR in the predator cells leading to increased expression of an antidote gene (ccdA
) to rescue the predators. The predator cells in turn produce another AHL, 3OC12HSL by LasI that permeates into the prey cells where it is bound by LasR, which activates expression of the killer ccdB
gene, inducing ‘predation.' This system satisfies the broader definition of predation for a two-species ecosystem, where one species (prey) suffers from the growth of the second (predator) and the latter benefits from the growth of the former (Taylor, 1984
). However, it differs from the canonical predator–prey system in two aspects. First, instead of acting as a food source, the prey provides an ‘antidote' to programmed death of the predator. Second, there is competition between predator and prey for nutrients in a co-culture, which is generally absent in natural predator–prey systems.
A synthetic predator-prey ecosystem.
A synthetic predator–prey ecosystem. (A
) The system consists of two types of engineered bacteria controlling each other's survival and death via two different QS signals. Outer boxes represent cell walls. Arrows represent activation or production; blunt arrows represent inhibition or killing. CcdB (‘B') is controlled by Plac/ara-1
in the predator and PluxI
in the prey. CcdA (‘A') is controlled by PluxI
. The lux
system is shown in blue and the las
system in green. The QS genes in the predator are controlled by PLtetO-1
, and those in the prey by Plac/ara-1
. Filled circles represent 3OC6HSL and filled diamonds represent 3OC12HSL. The predator–prey interaction is activated by IPTG. See main text and Supplementary information
for more details. (B
) Bifurcation analysis of the model in terms of the killing constant dc
)) and the growth rate kc
)) (left panel). The model construction is described in equations S14–S17. A line connecting Hopf bifurcation points demarcates three regions: (i) prey domination, (ii) predator–prey oscillation and (iii) predator domination. Typical population dynamics for specific parameter sets are illustrated for each region (right panel).
We employed ordinary differential equations to model the major kinetic events during the functioning of this circuit (see Supplementary information
and Figure S5
for details). The model predicts extinction, coexistence or oscillatory dynamics (Box 1B
) depending on rates of cell growth, cell death controlled by AHL, as well as production and degradation of AHL. According to bifurcation analysis (Box 1B
, left panel), the larger the CcdB killing rate constant (dc1
), the higher the likelihood of robust oscillatory behavior: increasing dc1
significantly broadens the parameter region that gives rise to oscillation (region bound by the loci of Hopf bifurcation points). Our experiments indicated that the original LacZα-CcdB fusion protein, which was used in our population control circuit (You et al, 2004
), might lack the necessary toxicity (data not shown). To overcome this limitation, we constructed a new LacZα′-CcdB protein by deleting 32 amino acids in LacZα, which resulted in an eight-fold increase in potency compared to the full-length protein (Supplementary Figure S2
). The LacZα′-CcdB fusion protein is referred to as ‘CcdB' in this paper.
It was previously shown that both active LasR and active LuxR can initiate gene expression at the luxI
) (Gray et al, 1994
) (Supplementary Figure S3
). Therefore, we used PluxI
to drive both the ccdA
gene in the predator cells (regulated by LuxR) and the ccdB
gene in the prey cells (regulated by LasR). To expedite rescue of the predator by its survival signal (3OC6HSL), the corresponding QS system was placed under the control of PLtetO-1
promoter (Lutz and Bujard, 1997
), which is constitutively active in the MG1655 E. coli
strain. This strategy simultaneously accelerates 3OC12HSL synthesis and thus prey killing. To ensure sufficient killing of the predator, we implemented the CcdB induction circuits in a high-copy-number plasmid (p15A origin of replication with a copy number of ~20 per cell). This synthetic ecosystem is activated by isopropyl-β-D
-thiogalactopyranoside (IPTG), which induces ccdB
expression in the predator and the QS system in the prey. To optimize protein expression levels and synchronize the timing of the circuit components, we tested many combinations of promoters, origins of replication and bacterial host strains, and narrowed our focus to one configuration (Supplementary Figure S4A
). Unless noted otherwise, the predator circuit was implemented in the E. coli
strain MG1655, the prey circuit in the strain Top10F′. Finally, a destabilized green fluorescent protein gene (gfpuv-lva
) was introduced into the predator (Supplementary Figure S4A
) to capture dynamic changes and differentiate the predator from the prey cell populations in mixed cultures.
We verified the basic function of each population independently using IPTG and the corresponding AHL supplied exogenously (). In the OFF culture containing neither IPTG nor 3OC6HSL, predator cells grew to a relatively high density in LBK medium (). Growth was significantly inhibited in the ON culture containing 1 mM of IPTG, which fully induced ccdB expression. However, predator cells induced with IPTG were rescued by 3OC6HSL (100 nM), which when bound by LuxR can initiate ccdA expression. Rescued predators (with IPTG and 3OC6HSL) grew slightly slower than those in the OFF condition (), possibly due to the metabolic burden involved in expression of circuit components, incomplete neutralization of CcdB by CcdA, or both. As shown in , prey cells were able to grow in the OFF culture (without inducers) as well as medium containing only IPTG but perished when both IPTG and 3OC12HSL were present. The IPTG-induced prey population grew slower compared to that of the OFF culture. This minor growth retardation could be a result of the metabolic burden of the circuit components, basal level CcdB expression due to crosstalk between 3OC6HSL and LasR, or both.
Figure 1 Individual growth behaviors (without interactions) of (A) predator and (B) prey cells in liquid media. For each condition, 6 ml LBK medium containing chloramphenicol and kanamycin was inoculated with a single bacterial colony and was divided into three (more ...)
This synthetic predator–prey system represents one of the most complicated synthetic circuits reported to date. Long-term characterization of such circuits, which is essential to the verification of the designed dynamics, presents a major technological challenge. Measurement of circuit performance in conventional macroscopic (test tube) cultures hinted at oscillatory dynamics (data not shown). But insufficient temporal resolution due to the arduous and time-consuming process involved in acquiring such data made these observations inconclusive. In addition, the short-lived stability of the circuit in the macro-sized populations (Desai et al, 2007
) could only provide a fleeting glimpse of the programmed behavior: function was lost before the circuit behavior could be appreciably observed. To this end, we resorted to a recently developed microchemostat platform (Balagadde et al, 2005
); a high-throughput screening technology that can inexpensively perform rapid characterization of synthetic circuits under a matrix of conditions with unprecedented capabilities including long-term, non-invasive measurements of microbial population properties under steady-state conditions with single cell resolution and advanced automation. Equally important, the microchemostat's micro-sized population (~102
cells versus at least ~109
in macroscale cultures) reduces the number of cell division events per unit time and hence slows down microbial evolution (Desai et al, 2007
). This aspect facilitates monitoring of programmed behavior of bacterial populations for hundreds of hours despite strong selection pressure to evade population control (Balagadde et al, 2005
Furthermore, the challenge of real-time monitoring of community composition was met by an engineering revision to the microchemostat reactor channel to enable single-cell-resolved fluorescence quantification (Supplementary Figure S1
). In addition, we reduced reactor volume from ~16 to ~9 nl to further reduce the bacterial population size in each reactor, which has been shown to help stabilize a population controller during long-term culturing (Balagadde et al, 2005
). This strategy also allowed us to increase the number of parallel reactors per chip from 6–14 for higher throughput.
illustrates the typical oscillatory dynamics of predator and prey populations in the microchemostat with a period of ~180 h (IPTG=5 μM, dilution rate (D)=0.1125 h−1). Initially, when the prey density was low, predator growth was inhibited due to constitutive expression of CcdB. Conversely, the prey grew to a high density (~35 000 cells per reactor) in ~10 h. The high prey density led to rescue of the predator population, which recovered from ~5000 cells per reactor to >20 000 cells per reactor in ~10 h. The increased predator density in turn resulted in killing of the prey cells, leading to a drastic reduction in the prey density (from ~35 000 cells per reactor to just a few cells per reactor). The decrease in the prey density led to insufficient rescue of predator cells and the subsequent decline in predator density at ~35 h. The prey population recovered ~10 h later. In summary, the system yielded oscillatory dynamics between the predator and prey populations that clearly resemble those observed in natural predator–prey ecosystems.
Figure 2 Long-term characterization of system dynamics using the microchemostat. (A) Typical oscillatory dynamics of the system with a period of ~180 h (IPTG=5 μM, dilution rate=0.1125 h−1). (B) System dynamics for varying IPTG induction (more ...)
The parallel design of microchemostat reactors greatly facilitates the investigation of circuit response to changing circuit parameters. illustrates a typical set of dynamics of the synthetic predator–prey ecosystem for various circuit induction levels. Without induction (IPTG=0), predator and prey had negligible circuit-driven interactions. Since the predator (MG1655 strain) grows faster than the prey (Top10F′ strain), the prey population was eventually (after ~50 h) driven to extinction (the fluctuations of the prey density around 0 cell per reactor were numerical artifacts of imaging processing; see Supplementary information
for details). Higher IPTG levels (
5 μM) did elicit oscillatory dynamics between the predator and the prey populations.
These results suggest the crossing of a Hopf bifurcation point using the IPTG level as the control parameter (also see Supplementary Figures S6 and S7A
). This aspect is illustrated with bifurcation analysis () that accounts for the effects of IPTG on 3OC6HSL production by the prey and CcdB expression by the predator. At sufficiently low IPTG levels, there is minimal programmed predation between two populations and the primary interaction is competition. As a result, the predator population dominates due to its growth advantage, leading to the washout of the prey. Increasing IPTG beyond a critical level enables crossing of a Hopf bifurcation point and elicits oscillatory dynamics.
The model further predicts two salient features of the predator behavior at high IPTG induction levels. First, the response dynamics are much slower compared to those at low IPTG induction levels (, inset). Second, these dynamics involve minimal predator densities that are close to zero. Both aspects are qualitatively consistent with experimental data (top panels of ; Supplementary Figure S7A
), although the data do not exclude the possibility that the system may eventually approach steady state. Because of the miniaturized volume of the microchemostat reactor (~9 nl), the predator (or prey) population is prone to falling below the detection limit (~75 cells per reactor) during an oscillation trough (, row 1, columns 1 and 2; from ~10 to ~80 h) and in severe cases, getting washed out (, row 1, column 3).
Similarly, we explored the impact of dilution rate (D
) on the system dynamics. illustrates coexistence dynamics at D
, suggestive of an initial phase of long-period oscillations. However, an increase in D
(to 0.2 h−1
) led to a phase of damped oscillatory dynamics with much shorter periods (~50 h). We observed similar response with a different circuit configuration: predator (Top10F′) and prey (Top10F′) (; also see Supplementary Figure S4B
). These populations initially coexist at a low dilution rate (D
). They switched to damped oscillations with short periods (~30 h) upon an increase in D
(to 0.23 h−1
). Further increase in D
(to 0.31 h−1
) at ~250 h led to predator washout. Bifurcation analysis qualitatively accounts for the experimental observations that an increase in the dilution rate can lead to a significant decrease in the period of oscillations (, inset). Simulations further show that an increase in D
can shift a slow, sustained oscillation (, bottom left panel) to a fast, damped oscillation (, bottom middle panel). Predator cells will be washed out if D
is too large (, bottom right panel).
Figure 3 Dependence of systems dynamics on dilution rate (D). (A) Experimental dynamics of predator and prey populations at different D in the microchemostat: for the pair of predator (MG1655) and prey (Top10F′). (B) For the pair of predator (Top10F′) (more ...)
These modeling results qualitatively account for the experimental observations of the dependence of systems dynamics on dilution rate. We note that the system's response to a shift in D
is highly sensitive to the prevailing conditions of operation. Upon switching D
from 0.11 to 0.2 h−1
, the system either exhibits a damped oscillation or approaches a steady state (the third row, Supplementary Figure S7A
, where IPTG=50 μM). The sensitivity of the damped oscillation to perturbations of a parameter or initial conditions can qualitatively explain this phenomenon (Supplementary Figure S7C
, insets). Since damped oscillations are transient dynamics, small variations in the control of the microchemostat dilution rate or different initial phases of the predator and prey populations can significantly impact the oscillation behavior.
Drawing inspiration from the nature, synthetic biologists have created circuits that operate in single cells (Elowitz and Leibler, 2000
; Gardner et al, 2000
; Atkinson et al, 2003
; Fung et al, 2005
; Levskaya et al, 2005
), in individual populations (Kobayashi et al, 2004
; You et al, 2004
; Basu et al, 2005
; Anderson et al, 2006
) and between populations (Brenner et al, 2007
). Our study represents a major advance beyond recent efforts in constructing synthetic ecosystems (Shou et al, 2007
; Weber et al, 2007
), in terms of system complexity and resolution of characterization. It provides a framework for exploring how naturally existing cellular components can be assembled to program and coordinate highly complex cellular behavior. In addition, synthetic ecosystems such as ours may serve as well-defined systems for exploring evolutionary and ecological questions regarding, for example, the generation and maintenance of biodiversity (Chesson, 2000
; Nowak and Sigmund, 2004
). Conventionally, typical experiments involving microorganism systems are limited to using nutrient concentration and dilution rate as the sole control parameters (Fussmann et al, 2000
; Becks et al, 2005
). In comparison, our system allows direct manipulation of intrinsic parameters, such as growth rate, death rate and strength of cell–cell communication. In this system, there is a clear and explicit mapping between the genetic constructs and the corresponding population dynamics, which may facilitate future studies of how molecular interactions are manifested in the temporal patterns of population dynamics, a central theme in ecology (Bohannan and Lenski, 2000
; Foster et al, 2007
). The generic nature of our system design makes it portable to other ecological interactions, including mutualism, competition, commensalism and amensalism (May, 1974
). The virtually unlimited configurations that are possible with these basic elements will allow us to further examine the interplay between the environment, gene regulation and population dynamics. With additional control over population mixing or segregation, it will be possible to program bacterial populations to mimic development and differentiation in multicellular organisms.