Our goal in this paper was to map the environmental fitness landscape, and to identify aspects of stochastic gene expression that are necessary for establishing the environment-fitness connection. This required developing a simple modeling framework to account for the mutual dependency of cell growth and gene expression, so that the behaviors of cells with multiple phenotypes could be better understood and controlled. To allow for predictive quantitative modeling and properly test the modeling predictions, it was necessary to minimize the biological complexity inherent in natural gene networks. Thus, we chose to focus on cell populations carrying an inducible synthetic gene circuit, in the spirit of the booming field of synthetic biology 
. The synthetic gene circuit controlled a bifunctional antibiotic resistance protein, which provided protection from Zeocin, and could also be visualized in single cells due to its fluorescence 
After determining fitness and subpopulation fractions in two sets of well-defined environments, we used this information to predict fitness and subpopulation fractions in a combination of these environments. This is non-trivial when phenotypic variation is present, because the cell population's average fitness is an aggregate of individual cell fitness values that are transient due to phenotype switching. Thus, to map the environmental fitness landscape, we defined two different types of fitness (transient cellular fitness and overall cell population fitness), which relate to each other as the gene expression in individual cells relates to the population average of gene expression. A similar distinction has existed for genetically diverse populations 
, - but here we defined these types of fitness for genetically homogenous populations of cells that carry the same noisy bistable synthetic gene circuit.
Instantaneous fitness was affected by two different factors, namely the growth-inhibitory effect of Zeocin and rtTA-associated toxicity. We also used the inducer ATc to separately control cellular memory, defined as the average lifetime of a gene expression state (or the inverse of spontaneous switching rates) in a constant
. We found that the instantaneous fitness and cellular memory jointly shape the overall gene expression distribution and determine overall cell population fitness in various environments in a non-trivial manner. Specifically, we argued that the synthetic gene network forces cells to explore an intracellular potential analogous to Waddington's landscape 
under the influence of noise. While exploring the gene expression landscape, cells are also forced to move on a non-genetic version of Sewall Wright's fitness landscape 
, and face selection through the instantaneous fitness associated with their gene expression. Mapping the environmental fitness landscape requires integrating the movement on these two different landscapes, entirely defined by the cellular memories and instantaneous fitnesses associated with specific cellular states. The interplay between cellular memory and noisy cellular fitness results in distributions and fitness values () that disagree with the naïve expectations illustrated in . In summary, a true understanding of the environment-fitness connection can benefit from a Darwinian perspective applied to heritable nongenetic variation 
. Mapping and computationally predicting the environment-dependent fitness landscape for this synthetic gene circuit might provide a potential roadmap to follow for other synthetic gene circuits and even for natural gene networks, to identify fitness-optimizing conditions relevant for microbial drug resistance, persistence, and virulence.
We strove to develop mathematical models that were as simple as possible, but no simpler, when describing population level gene expression. In isolation, the original goals of these models were to describe mean gene expression levels, stochastic transitions between cell states, rtTA-associated toxicity, and cell growth. We also studied how pairwise interactions between these biological components shape population level fluorescence distributions. Finally, we synthesized these modeling approaches into comprehensive population-level stochastic simulations of gene expression. By modifying previous stochastic models 
to describe gene activation rather than repression and to account for the cost of inducer-bound activator, cell division, and non-genetic memory we were able to reproduce the overall experimental distributions of PF cells by computer simulations (see Fig. S12 in Text S1
). In summary, non-genetic fitness differences between various gene expression states need to be carefully considered when studying instantaneous population-level gene expression measurements based on intracellular biochemical kinetics. Modeling approaches that do not incorporate fitness differences will lead to biased parameter estimates and misinterpretation of intracellular dynamics.
For example, an important implication of our results is that temporal averages over individual cell lineages can be dramatically different from the population means due to slight changes in their individual fitnesses (). Consider the bimodal population of PF cells with an approximately 1
1 ratio of low and high expressors. Using live cell imaging under the microscope in a microfluidic chamber with controlled environment 
, we have tracked individual cell lineages from this population (see Fig. S12C in Text S1
, and Video S1
). Although cell populations were initially composed of 45% high expressors, over 30 hours the cell lineages spent 63% of the time in the high expression state. As the period of observation approaches the time of cellular memory, we expect that the cell lineages will converge to a percentage of~97.5% in the high expression state. Thus, the population-average of gene expression measured at any time differs from the mean estimated by tracking individual cell lineages over time (). The same is probably true for other statistical measures of gene expression. This is important as some research groups have been obtaining single cell-level data from large cell populations at a given time 
, while other groups have been collecting gene expression data by tracking single live cells over time 
. Neither of these approaches is ideal, as both ignore information about division rates in expanding cell lineages. Our results indicate that gene expression statistics obtained from these different approaches will in general disagree, and may need to be reinterpreted in light of the fitness differences caused by gene expression.
Our work differs from earlier studies aimed to determine stochastic switching rates that optimize fitness in fluctuating environments, some of which were exclusively theoretical 
, while others included experiments, but did not map the fitness landscape at increasing levels of sustained stress 
. Seeking to predict optimal switching rates, these earlier studies assumed that switching rates could be altered without a fitness cost. By contrast, adjusting the memory in our system to favor of the high expression state involved an increasing fitness cost. Consequently, the switching rates that optimize fitness in our system differ from the prediction of earlier bet hedging studies, namely that maximal memory of the high expression state should optimize fitness for very long periods of drug exposure 
. In fact, we found the exact opposite. Minimal memory of the high expression state optimized fitness in our system, as long as a minuscule high expressor subpopulation existed to salvage the population after the onset of stress. This different optimum is caused by the cost associated with increasing the memory of the high expression state. While higher memory could salvage a larger fraction of the population immediately after stress onset, it would be unfavorable on the long run due to the increased rtTA toxicity (see Section 9 and Fig. S11 in Text S1
Is there a similar fitness cost in natural systems for increasing the memory of a stress-tolerant state? The answer appears to be “yes”. First, bistability of stress-tolerance and stress resistance states are likely wide-spread in microbes 
, due to implicit, growth rate-mediated positive feedbacks 
resulting from the expression cost or toxicity associated with various defense mechanisms. Examples include the toxin component of toxin-antitoxin systems 
expression for type II persisters 
, or the TetA protein from the tetRA
tetracycline resistance operon 
. While defense and fitness cost are caused by the same protein in these natural systems, this is not unlike our synthetic gene network, in which the sources of cost (rtTA
) and defense (yEGFP::zeoR
) are expressed from identical promoters, and are therefore identically regulated, - almost as if these genes were fused or were part of an operon. Increasing memory in such implicit, growth-coupled positive feedback systems implies an elevated level of toxic protein, which causes slower growth and thereby a higher fitness cost. Therefore, toxicity, toxin expression level and memory of high toxin expression are practically synonymous in these systems. It will be highly important to study at the single cell level how toxin expression coupled with toxicity controls the emergence of bistability in many microbial toxin-antitoxin systems. This was exemplified by a recent pioneering study on the hipBA
system in Escherichia coli
, where tuning the level of the HipA toxin resulted in the coexistence of normally growing and growth-arrested cells. The coexistence of such distinct phenotypes implies bistability somewhere in the network. While a growth rate-mediated positive feedback 
was not considered by 
, it would be important to investigate if its inclusion can cause bistability 
The potential implications of the “sweet spot” observed in our experiments are intriguing, and may relate to the extremely low fraction of persister cells observed in microbial populations 
. If increasing the memory of the drug-tolerant state is costly, it is advantageous to dedicate as small a fraction of the population as possible to the persister state. This minute persister population will then be able to reproduce during drug treatment at a minimal fitness cost. This is true as long as persisters have extensive memory and nonzero division rates like type II persisters 
, which may decrease or vanish as the fraction of persister cells and the memory of the persister state increases. Our experimental model and computational approach can therefore be useful in studying the properties of microbial populations with small fraction of resistant cells and bet hedging strategies.
In addition, our findings have important implications for interpreting the outcome of biological measurement techniques focusing on the average gene expression in the population of cells (such as gel blot assays or microarrays). Specifically, apparent “upregulation” following stress exposure may be due to nongenetic selection of a particular pre-existing subpopulation, and may have nothing to do with actual gene regulation. For example, at the “sweet spot” Zeocin treatment caused a drastic apparent increase in average gene expression in the population that did not involve any increase in protein levels in individual cells, but was simply due to the increase of the proportion of high-expressing cells in the population. This effect was due to the selection of high expressors in the presence of Zeocin. A similar nongenetic selection mechanism may also be implicated in the emergence of microbial or cancer drug resistance, where rare pre-existing drug-tolerant cells could be selected during chemotherapy, causing additional resistance to subsequent rounds of treatment. Consistent with these studies, our results suggest that even without underlying genetic changes, gene expression variability can generate a stable, selectable subpopulation of rare survivor cells 
that maintains resistance over hundreds of cell generations, and may serve as a reservoir of increasingly drug resistant mutants.
The PF circuit design can be easily extended up to mammalian systems and down to bacterial systems, due to the simplicity of the circuit design and the availability of equivalent circuit components in other organisms. Growing evidence supports the feasibility of transferring synthetic gene circuits across organisms 
, which might someday enable translational applications in the clinic 
. For example, a potential application for synthetic gene circuits that consist of a self-activating regulator controlling a drug resistance gene may be engineering the microbiome 
. All parts of the human body host immense microbial communities, with cell counts that exceed by an order of magnitude the number of our own cells 
. The microbiome can impact human physiology, and is required for vital functions of the human body. Consequently, rebalancing the microbiome may facilitate the reestablishment of homeostasis within the human body. Such rebalancing may be necessary after antibiotic treatment, which – in addition to the targeted pathogens – may drastically perturb the gut microbiome. For this reason, crucial members of the human microbiome could someday be engineered to survive or even proliferate during antibiotic treatment. Gene circuits similar to the one described above operating near the “sweet spot” may be capable of this task by ensuring near-optimal growth in both the absence and presence of antibiotic.
Another potential use of gene circuits similar to the one described above might be for in vitro
studies of metastable states of stem cells or adult progenitor cells. It is known that these cell types can stochastically transition between various metastable states, some more prone to differentiation than others 
. These stochastic transitions and their implications for differentiation may be difficult to study in vivo
. On the other hand, in vitro
conditions may drastically perturb or abolish these metastable subpopulations. Gene circuits similar to PF controlling a metastable-state specific gene (in addition to a drug resistance gene) may reestablish the desired subpopulation fractions in vitro
. Moreover, applying selective pressure by introducing a drug into the medium could adjust the stability of these transient states independently of their subpopulation fractions. Such a system could someday provide a constant, desired yield of transient-amplifying or differentiated cells for further studies or perhaps regenerative medicine 
Finally, this work highlights the importance of selective environments for characterizing the behavior of synthetic gene circuits, but also has implications for studying natural gene networks. Diverse environments can broaden the functionality of regulatory modules that respond to specific external signals. Yet, understanding the environmental response of gene networks and their host cell populations is not trivial. As the size of natural and synthetic gene networks demanding quantitative description increases, a new challenge will be to describe their behavior in various well-defined environments. It will be interesting to see how quantitative modeling can tackle the combinatorial complexity associated with tuning multiple environmental factors imposed on regulatory networks.