Bistability and oscillations in co-repressive synthetic microbial consortia

Abstract

Background: Synthetic microbial consortia are conglomerations of genetically engineered microbes programmed to cooperatively bring about population-level phenotypes. By coordinating their activity, the constituent strains can display emergent behaviors that are difficult to engineer into isogenic populations. To do so, strains are engineered to communicate with one another through intercellular signaling pathways that depend on cell density.

Methods: Here, we used computational modeling to examine how the behavior of synthetic microbial consortia results from the interplay between population dynamics governed by cell growth and internal transcriptional dynamics governed by cell-cell signaling. Specifically, we examined a synthetic microbial consortium in which two strains each produce signals that down-regulate transcription in the other. Within a single strain this regulatory topology is called a “co-repressive toggle switch” and can lead to bistability.

Results: We found that in co-repressive synthetic microbial consortia the existence and stability of different states depend on population-level dynamics. As the two strains passively compete for space within the colony, their relative fractions fluctuate and thus alter the strengths of intercellular signals. These fluctuations drive the consortium to alternative equilibria. Additionally, if the growth rates of the strains depend on their transcriptional states, an additional feedback loop is created that can generate oscillations.

Conclusions: Our findings demonstrate that the dynamics of microbial consortia cannot be predicted from their regulatory topologies alone, but are also determined by interactions between the strains. Therefore, when designing synthetic microbial consortia that use intercellular signaling, one must account for growth variations caused by the production of protein.

Graphical abstract

Keywords

synthetic biology / microbial consortia / quorum sensing / relaxation oscillations

INTRODUCTION

A major goal of synthetic biology is the construction of genetic circuits that endow cells and organisms with novel functions. Synthetic gene circuits provide the basis for technologies such as gene therapy [1], biofuel and biopharmaceutical production [2], and have a range of environmental applications [3]. In addition, they allow for unprecedented control of biological systems thus opening new avenues in biological research [4]. The majority of currently available synthetic gene circuits have been built within a single strain and operate at the single-cell level. However, to realize the full potential of synthetic biology we need to be able to design consortia of interacting cells and organisms. Cooperating cells can specialize and assume different responsibilities within a consortium [5]. This allows bacterial consortia to be more efficient, and have a wider range of functions than monocultures. In such consortia, the signals within and between bacterial populations shape the response of genetic networks within cells. The activity of the population, in turn, arises from the coordinated activity of individuals [6, 7]. Consortia can thus perform computations and make decisions that are far more sophisticated than those of a single bacterium [8, 9].

To understand the behaviors of naturally occurring microbial consortia and to engineer synthetic consortia for practical applications, it is necessary to develop mathematical and computational models that describe their behavior. Such multi-scale models must simultaneously describe transcriptional dynamics within cells, interactions between cells due to cell-to-cell communication, and population-level dynamics that arise as different cell types compete for limited resources [10, 11]. To further complicate matters, each of these levels of organization is linked to the others. Transcriptional dynamics within single cells are affected by intercellular signaling molecules. The concentrations of signaling molecules, in turn, are determined by gene network activity and total strain density. Strain density is also affected by transcriptional dynamics, as protein production can affect the growth rate of cells [12, 13, 14]. Therefore, to understand and predict the dynamics of microbial consortia, one must consider the dynamic interplay of multiple levels of organization.

Here we introduce a class of models describing the dynamics of synthetic microbial consortia in which two strains transcriptionally repress each other. We first introduce a deterministic model to describe the average behavior of cells within each population, their global interactions, and the resulting emergent dynamics. Furthermore, we show how to extend this model to include stochastic effects due to small population size and small molecular counts within each cell. This stochastic model has the deterministic model as its mean field approximation while it is able to capture fluctuations within cells and across strains as well as strain extinctions.

The co-repressive microbial consortium is a generalization of the well-known synthetic toggle switch that operates in a single bacterium [15] (see Figure 1A). In single cells, the co-repressive toggle switch can exhibit transcriptional bistability with two mutually exclusive gene expression states: i) expression of tetR and repression of lacI; or ii) expression of lacI and repression of tetR. The analogous synthetic co-repressive microbial consortium, shown in Figure 1B, might be constructed with two strains using two orthogonal quorum sensing systems (here cinI/R and rhlI/R) [16, 17], and two transcriptional inverters [18]. When one strain is active, it produces an intercellular signal that transcriptionally represses the opposing strain. Unlike its bistable single strain counterpart, the co-repressive consortium may exhibit more complicated behaviors. In particular, we show that the co-repressive consortium can oscillate if the growth rates of the strains depend on their transcriptional state, which can occur when heterologous protein is produced [12, 13].

While we only consider the dynamics of co-repressive microbial consortia, our modeling approach can be extended to any consortia of fixed size when spatial effects are negligible. Both our deterministic and stochastic models can be easily modified to describe different gene circuits and cell-to-cell interactions. We thus provide a general framework for modeling and analyzing the interplay between population and gene circuit dynamics that drives the behavior of microbial consortia.

RESULTS

The dynamics of the two-strain co-repressive consortium

We first consider the deterministic dynamics of a two-strain consortium growing in a small, well-mixed turbidostat, such as a microfluidic trap [19]. Since such traps have fixed volume, we assume that the total population size and cellular volume of the consortium remain constant. To maintain a fixed population size whenever a cell divides, we assume that a randomly chosen cell from the consortium exits the chamber. For simplicity, we assume that all cells in the consortium have equal size. The fraction of strain 1 within the chamber is then defined by the ratio r=V_{\mathrm{1}}/V=n_{\mathrm{1}}/N where V_{\mathrm{1}} is the volume occupied by strain 1, V is the total volume, n_{\mathrm{1}} is the number of cells in strain 1, and N is the total number of cells. Note that the ratio r is treated as a real number that can vary between \mathrm{0} and \mathrm{1}. Moreover, denoting the volume occupied by strain 2 by V_{\mathrm{2}} and the number of cells in strain 2 by n_{\mathrm{2}}, we have V_{\mathrm{1}}+V_{\mathrm{2}}=V and n_{\mathrm{1}}+n_{\mathrm{2}}=N. Therefore, the fraction of strain 2 in the chamber is given by \mathrm{1}-r=V_{\mathrm{2}}/V=n_{\mathrm{2}}/N.

We model the dynamics of a single gene within each strain. We assume that the production of an enzyme that creates a signaling molecule is repressed by the presence of the signaling molecule produced by the other strain. Further, the concentration of a signaling molecule is assumed to be spatially homogeneous and directly proportional to the total concentration of its enzyme within the culture (i.e., the product of the average intracellular concentration of the enzyme and the number of cells containing it). Therefore, the time evolution of the average concentrations x_{\mathrm{1}} and x_{\mathrm{2}} of the corresponding enzymes can be described by

(1)\begin{array}{l}{\dot{x}}_{\mathrm{1}}=\alpha f_{\mathrm{1}}(x_{\mathrm{2}},r)-{\beta }_{\mathrm{1}}{x}_{\mathrm{1}},\\ {\dot{x}}_{\mathrm{2}}=\alpha f_{\mathrm{2}}(x_{\mathrm{1}},r)-{\beta }_{\mathrm{2}}{x}_{\mathrm{2}},\end{array}

where

(2)f_{\mathrm{1}}(x_{\mathrm{2}},r)={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+{((\mathrm{1}-r)N{x}_{\mathrm{2}}/\theta )}^n}},f_{\mathrm{2}}(x_{\mathrm{1}},r)={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+{(rN{x}_{\mathrm{1}}/\theta )}^n}}.

Here \alpha is the maximal production rate of the proteins, \theta is a scaling parameter that determines half-maximal repression, and {\beta }_{\mathrm{1}} and {\beta }_{\mathrm{2}} are the growth rates of strains 1 and 2, respectively. The proteins are assumed to be stable and decrease in concentration only through cellular growth and division at a rate proportional to the growth rates. It is assumed that the volume is measured in units of single cells so that the total cellular volume is V=N. Moreover, we assume that there is a linear relation between the expression of the signaling molecule and protein within a single cell. Therefore, protein production in strain \mathrm{2} is repressed in proportion to the total protein signal V_{\mathrm{1}}{x}_{\mathrm{1}}=rV{x}_{\mathrm{1}}=rN{x}_{\mathrm{1}} from strain \mathrm{1}. Similarly, the protein production in strain \mathrm{1} is repressed by the total protein signal (\mathrm{1}-r)N{x}_{\mathrm{2}} from strain \mathrm{2}. Note that this approximation for the amount of signaling molecule assumes that: i) the transient dynamics of the signaling molecule are fast with respect to changes in the corresponding protein concentration x_i, and ii) the quasi-equilibrium concentration of the signaling molecule is linearly proportional to the amount of enzyme making it. The first approximation is generally valid provided that the growth chamber is small enough and that diffusion across cell walls is fast. The second approximation is valid provided that the presence of the signaling molecule does affect the enzyme’s ability to make it.

Next, we describe a deterministic model of the dynamics of the population ratio r assuming that the number of cells within each strain can be described by a birth-death process. A new cell is born at cell division, while a “death” occurs when a cell is removed from the chamber. Recall that n_{\mathrm{1}} and n_{\mathrm{2}} are the number of cells in strains \mathrm{1} and \mathrm{2}, respectively, so that n_{\mathrm{1}}+n_{\mathrm{2}}=N is constant. Left on their own, the two strains would grow exponentially with rates {\beta }_{\mathrm{1}}{n}_{\mathrm{1}} and {\beta }_{\mathrm{2}}{n}_{\mathrm{2}}, respectively. To keep the total population size constant, we set the total rate at which cells are removed from the population to {\beta }_{\mathrm{1}}{n}_{\mathrm{1}}+{\beta }_{\mathrm{2}}{n}_{\mathrm{2}}. If all cells are equally likely to be removed, the probability that a cell is removed from strain i is n_i/N, giving a death rate of ({\beta }_{\mathrm{1}}{n}_{\mathrm{1}}+{\beta }_{\mathrm{2}}{n}_{\mathrm{2}}){\displaystyle \frac{n_i}{N}}. The deterministic birth-death process of the strain \mathrm{1} can then be described by the differential equation

(3){\dot{n}}_{\mathrm{1}}={\beta }_{\mathrm{1}}{n}_{\mathrm{1}}-({\beta }_{\mathrm{1}}{n}_{\mathrm{1}}+{\beta }_{\mathrm{2}}{n}_{\mathrm{2}}){\displaystyle \frac{n_{\mathrm{1}}}{N}}.

Substituting n_{\mathrm{1}}=rN and n_{\mathrm{2}}=(\mathrm{1}-r)N into Equation (3), we obtain the logistic equation

(4)\dot{r}=({\beta }_{\mathrm{1}}-{\beta }_{\mathrm{2}})r(\mathrm{1}-r),

for the ratio r of the strain \mathrm{1} in the chamber.

Bistability in the absence of metabolic loading

There are many reasons why protein production might influence the growth rate of a cell. For instance, heterologous protein expression may slow growth due to metabolic loading [12, 13], or transcription of a suicide gene may be linked to one of the two states [20, 21]. For simplicity, we will refer to these phenomena collectively as “metabolic loading”—a burden imposed by the production of heterologous protein.

In the absence of metabolic loading, the growth rates of the two strains are not affected by the production of the enzymes and hence will remain constant. If the two strains grow at the same rate {\beta }_{\mathrm{1}}={\beta }_{\mathrm{2}}, Equation (4) implies that the ratio r will remain fixed at its initial value. Then we can treat r as a parameter in Equation (1) and the equilibria x_{\mathrm{1}}^{*} and x_{\mathrm{2}}^{*} are given by

(5)\begin{array}{l}{\beta }_{\mathrm{1}}{x}_{\mathrm{1}}^{*}=\alpha f_{\mathrm{1}}(x_{\mathrm{2}}^{*},r),\\ {\beta }_{\mathrm{2}}{x}_{\mathrm{2}}^{*}=\alpha f_{\mathrm{2}}(x_{\mathrm{1}}^{*},r).\end{array}

Indeed, the equilibria are the solutions of x_{\mathrm{1}}^{*}=h_{\mathrm{1}}(h_{\mathrm{2}}(x_{\mathrm{1}}^{*})), where h_i=\alpha f_i/{\beta }_i. Since h_{\mathrm{1}} and h_{\mathrm{2}} are monotonically decreasing sigmoidal functions, h_{\mathrm{1}}\circ h_{\mathrm{2}} is a monotonically increasing sigmoidal function and it intersects the diagonal in 1, 2, or 3 points which correspond to the equilibria.

Figure 2A shows a typical example of how x_{\mathrm{1}}^{*} changes with the ratio r. If the ratio r is high or low, there exists a unique equilibrium (solid line) and one of the strains dominates the trap; the dominant strain will be expressed and the opposite strain repressed. For mid values of r, the system is bistable with two stable equilibria (solid lines) and an unstable equilibrium in the middle (dashed line). This region of bistability is bounded by bifurcations at which two of the equilibria disappear in a saddle-node collision. That is, the consortium behaves like a toggle switch for a range of strain ratios. Unlike switches that operate on the level of single cells, bistability in the present case depends on the ratio of the strains in the trap.

Figure 2B shows simulations of Equation (1) with equal growth rates {\beta }_{\mathrm{1}}={\beta }_{\mathrm{2}} for two different sets of initial conditions where the ratio r is chosen from the bistable domain. If the initial conditions fall below the gray dashed separatrix (given by the stable manifold of the unstable equilibrium denoted by red ○), the system approaches the equilibrium denoted by blue ♦. The simulations starting above the separatrix approach the other stable equilibrium denoted by blue ■. Note that the slope of the separatrix is equal to \mathrm{1} for ratio r=\mathrm{0.5}, and decreases with decreasing r.

When the difference {\beta }_{\mathrm{1}}-{\beta }_{\mathrm{2}} between the growth rates of the two strains is small, r changes slowly according to Equation (4). Let us assume that

(6)\begin{array}{l}{\beta }_{\mathrm{1}}={\beta }_{\mathrm{0}}(\mathrm{1}+\epsilon ),\\ {\beta }_{\mathrm{2}}={\beta }_{\mathrm{0}},\end{array}

where \epsilon reflects the relative difference between the growth rates of the two strains. We only consider\epsilon \ge \mathrm{0}, since if the roles of strains \mathrm{1} and \mathrm{2} can be reversed. When \epsilon is small, Equations (1), (2), and (4) form a slow-fast system. Thus the concentrations x_{\mathrm{1}} and x_{\mathrm{2}} track a stable equilibrium until it disappears in a sadddle-node bifurcation [22] (see Figure 2A). Once an equilibrium disappears, the system jumps to the other stable equilibrium, and a switch in expression levels occurs. For example, as shown in Figure 3B, 3C, when \epsilon >\mathrm{0} and r is initially close to \mathrm{0}, the concentration x_{\mathrm{1}} initially approaches the lower equilibrium value (repressed state) and stays close to it as r increases slowly due to the larger growth rate of the first population. After the ratio r passes through the critical value at which the lower equilibrium disappears, x_{\mathrm{1}} switches to the higher equilibrium (expressed state). Figure 3B, 3C show such transitions for small and intermediate differences in growth rates. Notice that the transition occurs earlier in time as \epsilon increases. On the other hand when the growth rates are equal, no transition occurs as the ratio r remains constant (Figure 3A). While the dynamics for different values of \epsilon >\mathrm{0} are similar, the slow-fast approach is valid only when \epsilon is small.

When the two strains have identical growth rates, the two-strain microbial consortium behaves like a single-strain toggle switch: gene expression in the two strains is bistable for a range of population ratios r. However, if one strain has a larger growth rate, the opposite strain eventually disappears from the trap. Even before it is driven to extinction, the less numerous strain becomes fully repressed.

Impact of metabolic loading on population toggle dynamics

We next investigate the dynamics of the co-repressive consortium in the presence of metabolic loading, i.e., when the growth rates of the two strains depend on their transcriptional states. If a balance of population sizes is necessary to maintain a particular behavior, then any change in the growth rates can affect the dynamics of the consortium. To gain insight as to how metabolic loading in the two-strain toggle consortium can lead to relaxation oscillations, assume that the cells with higher expression rate experience an increased metabolic burden, and thus grow slower. This, in turn, allows cells in the repressed state to increase their relative population size. Once these cells dominate the trap, they are no longer repressed. As they reach high expression levels, they experience higher metabolic load and the process repeats.

To demonstrate such relaxation oscillations in the co-repressive consortium, we assume that the growth rates of both strains depend on the rate of expression of each gene, f_{\mathrm{1}} and f_{\mathrm{2}} defined in Equation (2). In particular, Equation (6) is replaced by

(7)\begin{array}{l}{\beta }_{\mathrm{1}}={\beta }_{\mathrm{0}}(\mathrm{1}+\epsilon )(\mathrm{1}-\rho f_{\mathrm{1}}(x_{\mathrm{2}},r)),\\ {\beta }_{\mathrm{2}}={\beta }_{\mathrm{0}}(\mathrm{1}-\rho f_{\mathrm{2}}(x_{\mathrm{1}},r)),\end{array}

where the parameter determines the impact of the metabolic load on the growth rates, such that \rho \approx \mathrm{0} corresponds to a low and \rho \approx \mathrm{1} corresponds to a high impact. The growth rate of a strain is therefore largest when gene expression is at its minimum. Here the parameter \epsilon determines the difference between the maximal growth rates of the two strains.

Figure 3 shows solutions of Equations (1), (4), and (7) for different values of the parameters \epsilon and \rho. When \epsilon =\mathrm{0}, the system exhibits oscillations for all values of \rho >\mathrm{0}, see Figure 4A‒4C (recall that for \rho =\mathrm{0} case, we have the bistable system shown in Figure 2A). The frequency and amplitude of these oscillations increase with the impact of the metabolic loading \rho. Figure 5A‒5C show the bifurcation diagrams while using \rho as the bifurcation parameter for the case \epsilon =\mathrm{0}. In Figure 5A, three unstable equilibria are shown as a function of \rho. Figure 5B shows the value of the population ratio r for these unstable equilibria. Figure 5C shows the amplitude of the oscillations as a function of \rho.

Figure 4D‒4F shows numerical simulations of Equations (1), (4), and (7) for different \rho values where \epsilon =\mathrm{0.25} is kept fixed. We observe that when \epsilon >\mathrm{0}, oscillations occur only when metabolic loading is sufficiently large. The bifurcation diagrams for the case \epsilon =\mathrm{0.25} are shown in Figure 5D‒5F. When \rho is close to zero the equilibrium r^{*}=\mathrm{1} is stable. An increased metabolic load, \rho, leads to the appearance of a stable equilibrium for which (red line in Figure  5D, 5E) that emerges via a transcritical bifurcation indicated by \times in Figure  5D, 5E. This equilibrium is stable over a small range of the parameter \rho, then undergoes a supercritical Hopf bifurcation (marked by *) leading to stable oscillations. Note that if we continue the equilibrium shown by the red branch to lower \rho values, the corresponding population ratio r^{*} is larger than \mathrm{1} and physically not meaningful. Therefore, this branch section is not shown here. The amplitude of the periodic orbit is shown in Figure 5F. The middle equilibrium (red branch) is also included in Figure 5F as a steady state with zero amplitude. The amplitude of the oscillations arising from the Hopf bifurcation grows explosively over a small parameter range after the bifurcation. We conjecture that this is due to a canard explosion [23].

Small population effects

In smaller traps, or confined geometries, the number of cells may be sufficiently small for stochastic fluctuations in population size and gene expression to become appreciable. We therefore describe a full stochastic model of the dynamics of intracellular proteins, as well as the birth and removal of cells in the two strains. We work under the same assumptions as in the previous section: the number of cells, n_{\mathrm{1}}+n_{\mathrm{2}}=N, in the trap is fixed, and is maintained by following each birth by a removal of a randomly chosen cell. We again neglect spatial effects.

We model protein dynamics separately in each cell. Let x_{i,j} be the number of proteins in cell j of strain i, so that i\in \{\mathrm{1},\mathrm{2}\}, and j\in \{\mathrm{1},\mathrm{2},\mathrm{...},n_i\}. We therefore have a set of N birth processes modeling protein production

(8)\begin{array}{llllll}\varnothing & \stackrel{\alpha f_{\mathrm{1}}}{\to }& x_{\mathrm{1},\mathrm{1}}& \varnothing & \stackrel{\alpha f_{\mathrm{2}}}{\to }& x_{\mathrm{2},\mathrm{1}}\\ \varnothing & \stackrel{\alpha f_{\mathrm{1}}}{\to }& x_{\mathrm{1},\mathrm{2}}& \varnothing & \stackrel{\alpha f_{\mathrm{2}}}{\to }& x_{\mathrm{2},\mathrm{2}}\\ & ⋮& & & ⋮\\ \varnothing & \stackrel{\alpha f_{\mathrm{1}}}{\to }& x_{\mathrm{1},n_{\mathrm{1}}}& \varnothing & \stackrel{\alpha f_{\mathrm{2}}}{\to }& x_{\mathrm{2},n_{\mathrm{2}}}.\end{array}

The rate of protein production in each cell in strain i is \alpha f_i, where

f_{\mathrm{1}}={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+{\left({\displaystyle \frac{x_{\mathrm{2},\mathrm{1}}+x_{\mathrm{2},\mathrm{2}}+\cdots +x_{\mathrm{2},n_{\mathrm{2}}}}{\theta }}\right)}^n}},\hspace{1em}

and

(9)\hspace{1em}{f}_{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+{\left({\displaystyle \frac{x_{\mathrm{1},\mathrm{1}}+x_{\mathrm{1},\mathrm{2}}+\cdots +x_{\mathrm{1},n_{\mathrm{1}}}}{\theta }}\right)}^n}},

cf. Equation (2). Note that we do not model protein degradation explicitly, as we assume that proteins are relatively stable. Approximately fixed concentrations are maintained by dilution through division, as explained further below.

We separately model the division (birth) and removal (death) of cells from the trap. The rates of these processes are as described in the derivation of Equation (3). The number of cells n_{\mathrm{1}} and n_{\mathrm{2}} in each strain is governed by the coupled birth and death processes

(10)\begin{array}{l}\varnothing \stackrel{{\beta }_{\mathrm{1}}{n}_{\mathrm{1}}}{\to }\mathrm{One\; cell\; in\; strain\; 1\; divides\; and}\{\begin{array}{l}\mathrm{a\; cell\; from\; strain\; 1\; is\; removed\; with\; probability}{n}_{\mathrm{1}}/N,\hfill \\ \mathrm{a\; cell\; from\; strain\; 2\; is\; removed\; with\; probability}{n}_{\mathrm{2}}/N,\hfill \end{array}\hfill \\ \varnothing \stackrel{{\beta }_{\mathrm{2}}{n}_{\mathrm{2}}}{\to }\mathrm{One\; cell\; in\; strain\; 2\; divides\; and}\{\begin{array}{l}\mathrm{a\; cell\; from\; strain\; 1\; is\; removed\; with\; probability}{n}_{\mathrm{1}}/N,\hfill \\ \mathrm{a\; cell\; from\; strain\; 2\; is\; removed\; with\; probability}{n}_{\mathrm{2}}/N.\hfill \end{array}\hfill \end{array}

The two population sizes therefore follow Moran dynamics [24, 25]: whenever a cell divides, another randomly chosen cell is removed.

We use Gillespie algorithm to sample trajectories from the processes described by Equations (8) and (10). At each step of the Gillespie algorithm, the possible events and their probabilities are given by these two equations. There are N+\mathrm{2} possible events that can occur in each step of the algorithm: a birth (division) in either strain, accompanied by the removal of a random cell, and a birth (formation) of a protein in one of the N cells. The rate at which a division occurs in strain i is given by {\beta }_i{n}_i where {\beta }_i are given by Equation (7), i.e.,

(11)\begin{array}{l}{\beta }_{\mathrm{1}}={\beta }_{\mathrm{0}}(\mathrm{1}+\epsilon )(\mathrm{1}-\rho f_{\mathrm{1}}),\\ {\beta }_{\mathrm{2}}={\beta }_{\mathrm{0}}(\mathrm{1}-\rho f_{\mathrm{2}}),\end{array}

where the protein production rates, f_{\mathrm{1}} and f_{\mathrm{2}}, are defined in Equation (9).

When a cell divides, its cellular material, including all proteins, is divided between the two daughter cells. Let the index pair (i,j) correspond to cell j in strain i, so that i\in \{\mathrm{1},\mathrm{2}\} and j\in \{\mathrm{1},\mathrm{2},\dots ,n_i\}. If an event corresponds to division in strain \mathrm{1}, we pick a random cell in that strain with index (\mathrm{1},b), to be divided into two. If the cell to be removed comes from the same strain, say it has index (\mathrm{1},d), we replace cell (\mathrm{1},d) with a new cell, and divide the proteins from cell (\mathrm{1},b) into two groups by sampling from a binomial distribution. One of the two groups is assigned to cell (\mathrm{1},b), while the other is transferred to the new cell with index (\mathrm{1},d). In this case, the cell counts n_{\mathrm{1}} and n_{\mathrm{2}} remain constant. If one of the cells in the strain 2 is chosen to be removed, say cell (\mathrm{2},d), we add an additional cell with index (\mathrm{1},n_{\mathrm{1}}+\mathrm{1}) to strain 1. We again partition the proteins from the dividing cell (\mathrm{1},b) between the mother and daughter cell. We then remove the cell (\mathrm{2},d) from strain 2 and renumber the remaining cells in strain 2 to close the resulting gap, that is x_{\mathrm{2},j}\leftarrow x_{\mathrm{2},j+\mathrm{1}}, for j=d,d+\mathrm{1},\dots ,n_{\mathrm{2}}. In this second case, n_{\mathrm{1}} increases by 1 and n_{\mathrm{2}} decreases by 1. The algorithm is equivalent if a cell divides in strain 2. Also a modeling assumption in this algorithm is that the birth and removal of cells happen concurrently, not allowing for the cell that is just born to be removed.

Next we explore the effects of stochasticity on the dynamics of the two-strain toggle consortia.

Stochastic dynamics in the absence of metabolic loading

With no metabolic loading, i.e., when \rho =\mathrm{0}, the growth rates of the two strains are constant. The slower growing strain is more likely to disappear from the trap. However, even if the two growth rates are equal, random fluctuations eventually lead to the extinction of one strain.

When the two strains of bacteria have equal growth rates, the probability of a strain dominating the whole population is equal to the initial proportion of that strain in the trap [25]. For example, the probability that the strain \mathrm{1} will take over is n_{\mathrm{1}}(\mathrm{0})/N. Figure 5 shows simulations of the stochastic model described by Equations (8)‒(11) with \rho =\mathrm{0} and when the two strains have equal growth rates, i.e., \epsilon =\mathrm{0}, for different population sizes N. In Figure 6A, 6B, the blue and green curves show the mean number of proteins in strains \mathrm{1} and \mathrm{2}, respectively. The gray curves represent the solutions of corresponding deterministic system described by Equations (1) and (4) with the same parameters. When the population size is smaller the effects of random fluctuations are more pronounced. This can be observed when comparing Figure 6A and 6B. In Figure 6A, where the total cell number N is larger, the protein concentrations x_{\mathrm{1}} and x_{\mathrm{2}} follow the deterministic model (shown by the gray curves) fairly close so that x_{\mathrm{1}} approaches the repressed state and x_{\mathrm{2}} approaches the expressed state. The ratio r oscillates around the initial value of r=\mathrm{0.4}. In Figure 6B, initially the genes in strain \mathrm{2} are expressed while the genes in strain \mathrm{1} are repressed. However, as the total cell population is smaller, the fluctuations in the population ratio r are larger and the likelihood of a switch taking place in the gene expression states is bigger. Therefore, we observe more switches between the gene expression states in Figure 6B. This behavior cannot be predicted by the deterministic model. In larger cell populations, the fluctuations in the population ratio are not strong enough and a switch in the gene expression states is less likely to occur. Furthermore, in a stochastic model random fluctuations always drive one population to extinction. However, in the deterministic model when growth rates are equal, no extinctions occur as the ratio remains constant.

Figure 6C, 6D show the time evolution of the population ratio r obtained by running 50 simulations of the stochastic model from Equations (8)‒(11). For N=\mathrm{200}, the variability in the population ratio is observed to be smaller than that observed in for N=\mathrm{40}. In the bottom plots in Figure 6C, 6D, the standard deviation of the population ratio ({\sigma }_r) and the mean ({\mu }_r) are shown versus time. In the case that the total population size is smaller (panel D), the variance of population ratio r grows faster with time.

When the two strains have different growth rates, the likelihood that the strain with a lower growth rate goes extinct increases. Figure 6 shows the simulations of the stochastic model from Equations (8)‒(11) for cases \epsilon =\mathrm{0.1} and \epsilon =\mathrm{0.25} each with population sizes N=\mathrm{200} and N=\mathrm{40}. Simulations of the deterministic model using the same parameters are also shown as gray curves. As predicted by the deterministic model, starting from \mathrm{0.1}, the population ratio r increases and crossing a critical value a switch happens in the gene expression states. However, the time at which this switch occurs depends on the random fluctuations in the population ratio and could be either before or after the time predicted by the deterministic model. It can also be the case that strain \mathrm{1} goes extinct (r becomes \mathrm{0}) and no switches occur. When \epsilon is larger (Figure 7C, 7D), the switch occurs faster and the stochastic dynamics are closer to those of the deterministic model.

The impact of metabolic loading on stochastic growth dynamics

We next simulate the stochastic model given by Equations (8) ‒ (11) taking into account the effect of metabolic loading on the growth rates of the bacterial strains by setting in Equation (11). As Figure 7 shows, the oscillations occur when the metabolic load is taken into account. However, when the total population size is small, the likelihood that one of the strains goes extinct increases. Therefore, as shown in Figure 7B, the oscillations can stop early due to extinction.

Next we ask whether metabolic loading affects the extinction time. Our deterministic analysis showed that metabolic loading increases the frequency of the oscillations in the population ratio r. Therefore, for smaller populations, one expects that the chance of extinction increases as the population ratio gets close to \mathrm{0} or \mathrm{1} more frequently. To test this hypothesis numerically, for each of the cases \rho =\mathrm{0} (no metabolic loading), \rho =\mathrm{0.5}, and \rho =\mathrm{0.9}, we run 500 simulations of the stochastic model from Equations (8)‒(11). The histograms of the extinction times are shown in Figure 8. Note that the extinction time is determined when either one of the populations goes extinct, i.e., the time at which r=\mathrm{0} or r=\mathrm{1}. We see that the distribution of the extinction times gets narrower as the metabolic load increases. For higher metabolic loads, the change in the mean and variance of the distribution is more pronounced.

DISCUSSION

We have shown that population growth can significantly alter the dynamics of synthetic microbial consortia. Differential growth between the strains which constitute the consortium, whether due to random fluctuations or changes in growth rate due to protein production, can lead to an imbalance in population sizes and alter the strength of signals between cells. Further, when the growth rates of cells are directly affected by protein production, “hidden” feedback loops arise that can change the dynamical landscape of the consortium. For instance, the regulatory structure of the two-strain co-repressive consortium forms a positive feedback loop. Positive feedback loops generally do not permit oscillatory solutions. However, we showed that growth rate changes due to metabolic load can create a hidden negative feedback loop that acts on a slow time scale. Therefore, the entire system has a fast positive feedback loop (due to signaling) and a slow negative feedback loop (due to metabolic loading)—the hallmark of relaxation oscillators [26]. A two-strain toggle consortium where the growth rates of the strains are affected by the metabolic load on the cells can therefore exhibit relation oscillations, a behavior that cannot be observed in the single-strain toggle switch. It is also shown that under the circumstances where the strains have an identical growth rate not affected by metabolic loading, the consortium shows bistability similar to the single-strain counterpart.

The perturbations due to the random partitioning of proteins at the time of cell division can have a strong effect on internal cell dynamics [27]. With metabolic loading internal and external fluctuations are even more strongly coupled. Variations in the ratio between the strains can change the expression within each cell. On the other hand, internal fluctuations within cells can affect growth, and thus the ratio between strains. Our model captures this interaction of fluctuations across scales, and can be extended to describe more details of local and global processes, or different dynamical behaviors.

In particular, our analysis could be extended to include spatial effects. Such effects will be most important for consortia that are not well-mixed or are large enough to limit intercellular signaling. As strains grow and compete for space within the colony, their spatial arrangement within the colony will change in time. Therefore, the regulatory “topology” of such a system will depend on both time and space, significantly complicating resulting dynamics. Any model that accurately recapitulates such a situation must include the internal dynamics of proteins within cells, the spatiotemporal dynamics of intercellular signals, the growth rate dynamics of the strains, and the time-dependent rearrangement of boundaries between cell types.

References

[1] Xie, Z., Wroblewska, L., Prochazka, L., Weiss, R. and Benenson, Y. (2011) Multi-input RNAi-based logic circuit for identification of specific cancer cells. Science, 333, 1307–1311.

[2] Zhang, F., Carothers, J. and Keasling, J. D. (2012) Design of a dynamic sensor-regulator system for production of chemicals and fuels derived from fatty acids. Nat. Biotechnol., 30, 354–359.

[3] Masiello, C. A., Chen, Y., Gao, X., Liu, S., Cheng, H.-Y., Bennett, M. R., Rudgers, J. A., Wagner, D. S., Zygourakis, K. Z. and Silberg, J. J. (2013) Biochar and microbial signaling: production conditions determine effects on microbial communication. Environ. Sci. Technol., 47, 11496–11503.

[4] Sprinzak, D. and Elowitz, M. B. (2005) Reconstruction of genetic circuits. Nature, 438, 443–448.

[5] Wintermute, E. H. and Silver, P. A. (2010) Dynamics in the mixed microbial concourse. Genes Dev., 24, 2603–2614.

[6] Chen, Y., Kim, J. K., Hirning, A. J., Josić K. and Bennett, M. R. (2015) Emergent genetic oscillations in a synthetic microbial consortium. Science, 349, 986–989.

[7] González, C., Ray, J. C., Manhart, M., Adams, R. M., Nevozhay, D., Morozov, A. V. and Balázsi, G. (2015) Stress-response balance drives the evolution of a network module and its host genome. Mol. Syst. Biol., 11, 827

[8] Regot, S., Macia, J., Conde, N., Furukawa, K., Kjellen, J., Peeters, T., Hohmann, S., de Nadal, E., Posas, F. and Sole, R. (2011) Distributed biological computation with multicellular engineered networks. Nature, 469, 207–211.

[9] Kong, W., Celik, V., Liao, C., Hua, Q. and Lu, T. (2014) Programming the group behaviors of bacterial communities with synthetic cellular communication. Bioresour. and Bioprocess., 1, 24

[10] Kanakov, O., Laptyeva, T., Tsimring, L. and Ivanchenko, M. (2016) Spatiotemporal dynamics of distributed synthetic genetic circuits. Physica D, 318– 319, 116–123.

[11] Blanchard, A. E., Liao, C. and Lu, T. (2016) An ecological understanding of quorum sensing-controlled bacteriocin synthesis. Cell. Mol. Bioeng., 9, 443–454.

[12] Tan, C., Marguet, P. and You, L. (2009) Emergent bistability by a growth-modulating positive feedback circuit. Nat. Chem. Biol., 5, 842–848.

[13] Scott, M., Gunderson, C. W., Mateescu, E. M., Zhang, Z. and Hwa, T. (2010) Interdependence of cell growth and gene expression: origins and consequences. Science, 330, 1099–1102.

[14] Nevozhay, D., Adams, R. M., Van Itallie, E., Bennett, M. R. and Balázsi, G. (2012) Mapping the environmental fitness landscape of a synthetic gene circuit. PLoS Comput. Biol., 8, e1002480

[15] Gardner, T. S., Cantor, C. R. and Collins, J. J. (2000) Construction of a genetic toggle switch in Escherichia Coli. Nature, 403, 339–342.

[16] Miller, M. B. and Bassler, B. L. (2001) Quorum sensing in bacteria. Annu. Rev. Microbiol., 333, 1315–1319.

[17] Wu, F., Menn, D. J. and Wang, X. (2014) Quorum-sensing crosstalk-driven synthetic circuits: from unimodality to trimodality. Chem. Biol., 21, 1629–1638.

[18] Tabor, J. J., Salis, H. M., Simpson, Z. B., Chevalier, A. A., Levskaya, A., Marcotte, E. M., Voigt, C. A. and Ellington, A. D. (2009) A synthetic genetic edge detection program. Cell, 137, 1272–1281.

[19] Bennett, M. R. and Hasty, J. (2009) Overpowering the component problem. Nat. Biotechnol., 27, 450–451.

[20] You, L., Cox, R. S. III, Weiss, R. and Arnold, F. H. (2004) Programmed population control by cell-cell communication and regulated killing. Nature, 428, 868–871.

[21] Balagaddé F. K., Song, H., Ozaki, J., Collins, C. H., Barnet, M., Arnold, F. H., Quake, S. R. and You, L. (2008) A synthetic Escherichia coli predator-prey ecosystem. Mol. Syst. Biol., 4, 187

[22] Hek, G. (2010) Geometric singular perturbation theory in biological practice. J. Math. Biol., 60, 347–386.

[23] Krupa, M. and Szmolyan, P. (2001) Relaxation oscillation and canard explosion. J. Differ. Equ., 174, 312–368.

[24] Moran, P. A. P. (1958) Random processes in genetics. Math. Proc. Camb. Philos. Soc., 54, 60–71.

[25] Nowak, M. A. (2006) Evolutionary Dynamics: Exploring the Equations of Life.Brighton: Harvard University Press

[26] van der Pol, B. (1926) LXXXVIII. On “relaxation-oscillations”. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 2, 978–992

[27] Veliz-Cuba, A., Gupta, C., Bennett, M. R., Josić K. and Ott, W. (2016) Effects of cell cycle noise on excitable gene circuits. Phys. Biol., 13, 066007

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg