Structure and Properties of the Network with Mutually Inhibitory Operons
To study adaptive responses to environmental changes without signal transduction machinery, we constructed a synthetic gene network that can exhibit bistability and monitored its ability to switch between two stable attractor states in response to evnironmental change. Gene regulatory circuits that are based on mutually inhibiting operons have previously been reported 
. However, unlike existing work on synthetic networks (for example, 
), for our purposes, the network need to express not only state-reporting fluorescence proteins, but also proteins with selectable phenotypes. To obtain a bistable gene network that is linked to metabolic phenotype switching we used the following two mutually inhibitory operons (construct pALL7, ). Operon1 is composed of the tetA
, Lac Repressor gene (lacI
, green fluorescence protein (GFP) gene (egfp
and a mutant glutamine synthetase (GLS-H) gene (gls-h
; the original post-translational regulation was eliminated) 
. Operon2 is composed of the trc
, Tet Repressor gene (tetR
, red fluorescence protein (RFP) gene (dsred.T4
, and the mouse dihydrofolate reductase (mDHFR) gene (mdhfr
. The two fluorescent proteins are introduced to monitor the transcriptional level of the two operons. GLS-H and mDHFR confer a metabolic phenotype to each of the two attractor, since these enzymes compensate for two distinct conditions of nutrient depletion.
The Plasmid Structures of the Mutual Inhibition Network
As a basic property of the circuit architecture, the cells harboring pALL7 can be in a monostable and a bistable behavioral regimes, depending on culture conditions. ( and ). After several serial overnight cultures in Medium N (see Experimental procedures), which does not impose restrictions on essential nutrients, the cells proliferated sufficiently fast so that the gene products of the two operons, including the repressors, were kept low due to dilution. Consequently mutual inhibition was too weak and did not produce bistability. In this monostable regime cells were reproducibly distributed around low levels of expression for both operons (blue dots). Two-color flow cytometry analysis () indicate low levels of expression of Operon 1 (green fluorescence) and Operon 2 (red fluorescence) for each cell. We refer to this state of weak expression from both operons as Attractor W, which is expected to appear at the low total concentration of the two repressors relative to the dissociation constants for their interacton with their promoter binding sites.
The Three Possible Attractors Generated by the Mutual Inhibition Network under Different Conditions
In contrast, when the total concentration of gene products is high, the network was shifted into the bistable regime, exhibiting the two Attractors 1 and 2, as expected (); either strong expression of Operon 1 with strong repression of Operon 2 (Attractor 1), or vice versa (Attractor 2). This was achieved by slowing down cell proliferation using 170 µg/ml nalidixic acid 
which reduces the specific growth rate by 30~40%. In this bistable regime, the balanced state of Attractor W is no longer stable. When either of the operons happens to express its encoded repressor at a slightly higher level than the other, the other operon is slightly suppressed, which in turn decreases the concentration of the repressor for the former operon, leading to a further increase in expression of the former. Thus, the cells are tipped into either of the two self-stabilizing Attractors 1 and 2, as originally intended by the circuit design.
Shift in Gene Expression Specifically toward the Adaptive Attractors
We next studied the “adaptive” response of the network to external changes by exposing the cells to culture conditions that require the presence of either enzyme (GLS-H or mDHFR) whose mutually exclusive expression is associated with the two attractors. Thus, we ask whether cells can find the “adaptive attractor” that copes with the nutrient condition (Scheme in ). For this purpose, we used two environmental conditions to implement the respective nutrient depletion: Medium M lacks glutamine, so that cells are required to syntheszie it to keep up cellular activity. Cells carrying pALL7 can overcome glutamine depletion if glutamine synthetase (GLS-H) in Operon 1 is stably expressed, that is, when they are in Attractor 1. Conversely, Medium T consists of Medium N plus trimethoprim lactate, which causes tetrahydrofolate depletion 
in the host cells. In this case the host cells carrying pALL7 can overcome tetrahydrofolate depletion if mDHFR in Operon 2 is expressed, which is active when cells are in Attractor 2. In summary, the Attractor 1 (2) is adaptive in mdeium M (T), respectively.
We put the cells carrying pALL7 to around Attractor W in the monostable regime by overnight culture in Medium N and then transferred them to either Medium M or T supplemented with 170 µg/ml nalidixic acid. In the flow-cytometry measurements, we observed that each cell underwent a unidirectional shift in gene expression pattern toward the attractor expressing the adaptive enzyme in response to respective nutrient depletion (). To monitor the time evolution of the shift we plotted the number of cells per unit volume against the green fluorescence/red fluorescence (G/R) ratio which reflects relative transcriptional activities from the two operons. At 0.5 h after transfer to either Medium M or T, the cells initially retained gene expression patterns, staying around Attractor W. Thereafter, in Medium M the distribution of cells shifted toward Attractor 1. This attractor is adaptive in medium M since glutamine synthetase in Operon 1 is expressed and compensates for the glutamine depletion in this medium. Conversely, the population transferred to Medium T was pulled toward Attractor 2. This shift is also adaptive since the high mDHFR expression compensates for tetrahydrofolate depletion in Medium T. Note that the control transfer to Medium N supplemented with the same concentration of nalidixic acid resulted in a stochastic shift of gene expression toward that of either Attractor 1 or 2 (). Therefore, the nutrient depletions in Medium M and T caused a selective, unidirectional shift toward the corresponding adaptive attractors.
The Initial Time Course of Adaptive Response in Gene Expression
Adaptive Response to Changing Environments
We next examined how cells adjust their gene expression adaptively to fluctuating environments. Cells carrying pALL7 were subjected to serial overnight culture with sequential medium changes in two different orders (). In these experiments, nalidixic acid was omitted from all three media to minimize the risk of unexpected genomic rearrangements during long-term cultivation. When cultured in Medium N, most of the cells exhibited a neutral position of G/R, i.e., were in Attractor W. Following transfer to Medium T (), the cells displayed a low G/R, indicating strong expression of Operon 2, encoding mDHFR, which compensates for the depletion of tetrahydrofolate in this medium. After return to Medium N, the cells again displayed the neutral position. The cells showed successive changes in G/R: neutral→high→neutral, when the medium was changed in the order N→M→N; consistently, the high G/R ratio here indicated expression from Operon 1, which encodes glutamine synthetase that compensates for the lack of glutamine in Medium M. Further transfer to Medium T caused the G/R ratio to shift back to the low value, indicating that the cells retained the ability to adjust their gene expression adaptively to Medium T. When the medium was changed in a different order (), the same three distributions appeared in response to the three media.
Attractor Selection in Changing Environments
Microscopic analysis confirmed that at the single cell level these reproducible distributions with the peaks of high, neutral, and low G/R correspond to Attractor 1, Attractor W, and Attractor 2, respectively (). The cells in Medium N expressed GFP and RFP weakly, indicating Attractor W. On the other hand, cells in Medium T showed strong expression of RFP and repression of GFP compared to those in Medium N, indicating Attractor 2. Similarly, cells in Medium M exhibited strong expression of GFP and repression of RFP, consistent with Attractor 1. The alternative appearance of the three attractors in each medium was also supported by mRNA quantification; the average numbers of mRNA copies transcribed from Operon 1 and Operon 2 per cell were 10.0 and “undetectable”, respectively, in Medium M, whereas they were “undetectable” and 6.0, respectively, in Medium T, and “undetectable” and 0.3, respectively, in Medium N. These results clearly show that the cells chose the adaptive attractors reproducibly in fluctuating environment.
Microscopic Examination of Cells Cultured Overnight in Three Different Media
The Response Is Not Due to Proliferation of Fit Cells
A fundamental question of cell state switching is whether (i) an external signal that causes a commitment to the expression of a particular phenotype does so by somehow instructing the gene transcription apparatus to express the appropriate (set of) gene(s) in all cells, or (ii) the signal merely promotes survival and expansion of the few cells that “happen” to express that desired phenotype 
. In other words, if the switching is only due to the scenario (ii), the two attractors are selected with equal probability first but only the cells in the adaptative attractor will proliferate, causing a shift in the population distribution. We found that the scenario (ii) alone cannot explain the observed macroscopic shift toward adaptive attractor, but the scenario (i) indeed is necessary, as demonstrated by careful monitoring of the time course in which cells redistribute between each of attractors, as explianed in the following.
At 0.5–2 h, when the total cell concentration had hardly increased as indicated by the minimal increase in the total areas under the curve (AUC) of distribution (), the cell distribution began to shift toward the adaptive attractor. To better represent the shift in population (change in histogram of cell number for a given G/R ratio) shows the increment in cell number within each bin in the histogram up to each sampling time, obtained by subtratcing the cell number from that at the previous sampling time point in . In Medium M, the subpopulation of cells with high G/R ratio clearly increased within the time period from 0.5 to 2 h (1.5×107 cells/l in the G/R ratio between 100.3 and 10), while that with low G/R ratio decreased (1.2×107 cells/l in the G/R ratio between 10−0.5and 100.3). The small difference in the total AUC between the downward peak for cells with low G/R ratio and the upward peak for cells with high G/R ratio indicates little cell growth from 0.5 to 2 h. Since nutrient starvation for less than 5 h caused little cell rupture (data not shown), the two allmost equal AUCs of downward and upward peaks indicate that each cell underwent a unidirectional alteration in its gene expression pattern from Attractor W toward Attractor 1 during this time period. Taken together, this kinetic analysis suggests that the observed unidirectional shift until 2 h can be attributed not to the prolification of cells randomely tipped toward the adaptive attractor (scenario (ii)) but to a deterministic change in gene expresssion of each cell (scenario (i)).
In contrast, for the period from 2 to 5 h, the change in gene expression toward Attractor 1 and the growth of cells at high G/R occurred simultaneously, as indicated by the smaller total area of the downward peak for cells with low G/R ratio compared with that of the upward peak for those with high G/R ratio. From 5 to 7.5 h, as indicated by the small downward peak for cells with the low G/R ratio, the change in gene expression toward Attractor 1 was almost complete, while cells with the high G/R ratio continued to grow. Thus, the time course analysis indicates that after exposure to Medium M the cells that happened to have been tipped toward the adaptive attractor at the beginning were pulled to these attractors and started growing in the favorable condition (scenario (ii)) while the other cells that happened to have been tipped to the non-adaptive attractor escaped from them and switched to the adaptive attractor in a scenario (i) process. Upon exposure to Medium T, the cells showed the same time course of attractor selection but in the opposite direction. These unidirectional changes in gene expression suggest that cells are capable of response adaptively to environmental changes through attractor selection. Thus, within the limits of measuring the time course of gene expression for each cell, both mechanisms appear to play a role. It remains to be determined to what extent each of two scenarios contribute to the observed adaptive attractor selection.
Fitness-Induced Gene Expression without Signal Transduction Machinery
Given that the mechanism of regulation of selecting the attractor cannot be solely explained by selective growth advantage of randomly swicthed cells, one may suspect that adaptive attractor selection in response to nutrient depletion may be due to some direct but ‘hidden’ instructive signal transduction that may somehow be “hard-wired” in the host cell's genome, connecting nutrient situation to the respective cis-region of the two metabolic genes. To exclude this possibility, we carried out a ‘promoter swap’ experiment. The mdhfr
genes were swapped in the pALL7 plasmids, so that the same regulatory circuit was maintained (plasmid pALL8, Figure S1A
). This new plasmid was introduced into the same host cells, E. coli
OSU1. The transformed cells were subjected to a series of overnight cultures in essentially the same media with a slight modification (legend of Figure S1
). As shown in Figure S1B
, attractor selection took place in the same way with the mutually inhibitory network of the pALL8 plasmid as with the pALL7 plasmid but only with the difference that with pALL8, the correlation between the media and the strongly expressed operon was the reverse of that observed with pALL7.
We did not observe a non-adaptive response with pALL8, in that for instance, depletion of tetrahydrofolate would have activated Operon 2, which contains the promoter that drives mDHFR in pALL7 but glutamin synthetase in pALL8. Thus, the promoter swap experiments corroborate the principle that expression of operons was dependent on the fitness conveyed by the encoded enzyme to the environmental condition and was not due to some unexpected fixed functional relation (such as an unknown signaling pathway) between the DNA sequence of these nominally generic promoter regions and the environmental signal.
In summary, we have shown experimentally that by introducing a gene network with bistable attractors expressing phenotypes that are sensitve to the environment, most of the cells changed the attractor state and associated gene expression in response to the environment, so as to mount an adaptive response by producing the enzymes necessary to compensate for the nutrient depletion.
Theoretical Model for Adaptive Attractor Selection
How can we understand the above adaptive attractor selection without specific signal transduction machinery? Here, we propose that such selection is a rather general process allowing cells to grow by selecting and maintaining stochastic gene expression patterns. The basis for the alternative expression patterns in the present gene network is a deterministic bistable system. Thus, we propose here first a model for bistability for our network that links attractor state with a metabolic phenotype.
We analyzed a standard model of mutually inhibitory operons which can generate bistability 
, but introduced terms to capture the phenotypic consequence of a stable state, which determines the adaptation to an environment, and a noise term:
1 and m
2 are the concentrations of the mRNAs or their protein products (which are here lumped together into one term since the phenotype switch of the population occurs at a slower time scale than transcription), transcribed from Operon1 and Operon2, respectively. S
) and D
) are the rate coefficients of synthesis and degradation and/or dilution due to the cell volume growth, respectively. Importantly, they depend on A
, which represents cellular activity. η1
represent independent white noise in gene expression.
By setting the equations (1) without noise to dm
0, one can obtain the fixed-point solutions for the system. This yields a single attractor of m
2* for S/D<0.5
(corresponding to Attractor W in the monostable regime in the experiment in the absence of nalidixic acid). For S
>0.5, there are two attractors satisfying (m
*) and (m
*) (see Supporting Information). Without the noise term, the initial condition with m
2 is attracted to the attractor located in the region of m
2 and vice versa.
To capture the phenotypic consequence that allows adaptation we introduce the variable A
“cellular activity”, which is increased when cells approach the adaptive attractor, expressing the genes that allows survival and optimal growth in a given environment. It is not easy to define “cellular activity” on a material basis, as it is a complex function of the concentrations of ATP and other chemicals. However, cellular activity can be lumped together and represented as the variable A
that increases monotonously with cell growth rate, particularly for unicellular organisms.
The first important postulate in the model is that both synthesis S
) and degradation and/or dilution D
) are increasing functions of activity A
(albeit in distinct ways) which in turn is correlated with the nutrition condition and growth rate. Since the rate of metabolic processes increase with activity A
, the amplitudes of the synthesis rate of mRNA or its protein product, S
), is expected to increase with A
. As the degradation rates of LacI and TetR are much smaller than the specific growth rate of cells in this work 
, the amplitude of D
) is mostly determind by the dilution due to cell growth and thus increases with A
The different attractors show different activities. The adaptive attractor has higher activity, i.e., larger values of S(A) and D(A). In other words, the deterministic part of Eq. (1) (i.e., all terms without the noise terms) take larger (smaller) values for an adaptive (non-adaptive) attractor, respectively. (Recall S(A) and D(A) are increasing functions of the activity A).
We now discuss selection of adaptive attractor. First we note that such selection of attractors is not possible in the above deterministic dynamics alone. However, gene expression is always accompanied with considerable random fluctuations, or noise 
, represented by the term η
which can account for noise-driven transitions between attractor states 
The second basic postulate in the model is that the noise amplitude is independent of the activity A, or at least it does not vanish with the decrease of the activity. Specifically, in a first approximation we assume a constant amplitude of noise. However, the specific form of noise in the Langevin Eq. (1) is not important, as long as a considerable amount of noise remains when the network is in the non-adaptive attractors, even if noise strength may depend on (m1, m2), or activity to some extent.
In general, for noise arising from chemical reactions associated with the activity, its strength should increase with the activity. However, even if the activity vanishes, fluctuation in gene expression may remain due to basal (housekeeping) biochemical processes in the cell. Indeed, recent measurement of noise on a variety of environmental conditions suggests that considerable amount of noise remains independent of the growth speed of a cell that is correlated with cellular activity 
. Hence, it is natural to assume that some part of noise is activity independent, even though the actual form of noise strength is difficult to predict at the present stage. In our experiment, we found that considerable noise remains even for a cells in the “non-adaptive condition” where they are in a state with low activity. When cells located around Attractor 1 in Medium M were transferred to Medium T (), most of the cells moved to Attractor 2, while a small fraction of cells remained at Attractor 1. The observation that the latter cells, though they had very low activity around the non-adaptive attractor, maintained a similar scattering in fluorescence intensity (Figure S2B
,) suggests that noise is present even in states with low activity. On the other hand, we found that the relative contribution of noise to the deterministic metabolic rate decreased as cells approached the adaptive attractors. In , the distribution of the G/R ratio became narrower as cells shifted their gene expression toward the adaptive attractors. The relative variances in G/R ratio (estimated as full width at half-maximum of the histogram) decreased by 72% (in Medium T) or 57% (in Medium M) by the time of 7.5 h, suggesting that the relative magnitude of the noise term decreased while that of the deterministic metabolic term became large with increasing activity.
We now describe how, given the above two basic postulates, the adaptive attractor is selected. Let us first consider the case in which an environmental change causes a marked decrease in A because the gene expression pattern is inappropriate (i.e., the cell is in the non-adaptive attractor). Then, the deterministic metabolic rate in Eq. (1) will be so small that it will approach the same magnitude as that of the noise term, η. The dynamics of gene expression will therefore be dominated mostly by random fluctuations. This is true as long as the network is in a region with the state of low activity. On the other hand, when the network moves into a region with high acttivity A, i.e., to the adaptive attractor, then the metabolic rate inceases, and the deterministic part becomes much larger than the noise term. Consequently, the dynamics of the system will be governed by the deterministic part of Eq. (1). Thus, whichever gene expression state the network starts with, gene expression will continue to fluctuate until it arrives at the adaptive attractor which is more stable against the noise because of the relatively lager metabolic rate of the first and second terms in Eq. (1).
With a simple form of D(A) and S(A), the above dynamics for attractor selection was confirmed by numerical simulation of Eq. (1) (). (for details on the mathematical model, see Supporting Information).
Simulation of Dynamics of Networks with Mutually Inhibitory Operons
The present scheme of adaptive attractor selection generally works for a system with dx/dt
, with f(x, A
),…} where fi
) is the metabolic rate of mRNA or protein for gene i. In general, if |f(x, A
)| increases with A
, the attractor with small A
is disturbed by noise so that an attratcor with high A
Switching between attractors has been studied in dynamical systems as noise-induced selection of attractors 
. In these previous studies, however, the direction of switching was not adaptive. In contrast, in the case presented here, due to the coupling of the attractor states with the cellular activity which influences the rate of the turnover of the deterministic state variables, an adaptive attractor with its higher cellular activity is inevitably selected in the presence of noise.