## Abstract

The lactose operon regulation in *Escherichia coli* is a primary model of phenotypic switching, reminiscent of cell fate determination in higher organisms. Under conditions of bistability, an isogenic cell population partitions into two subpopulations, with the operon's genes turned on or remaining off. It is generally hypothesized that the final state of a cell depends solely on stochastic fluctuations of the network's protein concentrations, particularly on bursts of lactose permease expression. Nevertheless, the mechanisms underlying the cell switching decision are not fully understood. We designed a microfluidic system to follow the formation of a transiently bimodal population within growing microcolonies. The analysis of genealogy and cell history revealed the existence of pre‐disposing factors for switching that are epigenetically inherited. Both the pre‐induction expression stochasticity of the lactose operon repressor LacI and the cellular growth rate are predictive factors of the cell's response upon induction, with low LacI concentration and slow growth correlating with higher switching probability. Thus, stochasticity at the local level of the network and global physiology are synergistically involved in cell response determination.

## Visual Overview

### Synopsis

The bacterium *Escherichia coli*, like many other microorganisms can use different sugars as a carbon source and uses some of these sugars in preference to others. For example, when grown in the presence of both lactose and glucose, the bacteria first consume glucose and use lactose only when glucose is exhausted. To this end, the enzymes necessary for lactose uptake and metabolism, grouped in one transcriptional unit called the *lac* operon (*lacZ*, *Y*, *A* encoding for the lactose degrading enzyme (β‐galactosidase), permease and transacetylase, respectively) are produced only in the absence of glucose and in the presence of lactose or its analogs, such as the non‐metabolizable analog thiomethyl‐β‐galactoside (TMG). In the absence of such inducers, the transcription of the operon is inhibited by the repressor molecule LacI. This inhibition is relieved by the inducers, which bind and inactivate LacI, initiating an amplifying feedback loop through the expression of the permease that ensures a high influx of inducer to maintain the operon's expression in the ‘on’ state. This phenomenon of adaptive enzyme production has been widely studied since its discovery by Jacques Monod and François Jacob and is one of the most famous and best characterized examples of transcriptional gene regulation. *E. coli* lactose operon is also a paradigm of cellular differentiation. Indeed, in the presence of an intermediate concentration of TMG, an isogenic bacterial population is divided in two subpopulations of cells with the operon's genes either turned on or remaining off. The differentiation step is generally hypothesized to depend on fluctuations in expression of the operon's proteins. Nevertheless, it is still poorly characterized. On the basis of experimental and theoretical approaches, we explored the determinants of cell fate in this system.

We designed a microfluidic device allowing the observation of single cells growing within a microcolony under conditions that can be changed at will. We used this setup to study phenotypic variability in the lactose operon induction under conditions leading to a transient bimodality of *lac* expression in the population. We used an *E. coli* strain modified to express the yellow fluorescent protein (YFP) and the cyan fluorescent protein (CFP), both under the control of a promoter regulated by LacI (P_{L}lacO1). Therefore, yellow and cyan fluorescence intensities both represent the concentration of active repressor molecules and indirectly, the expression state of the lactose operon. Microcolonies originating from a single cell were grown in the microfluidic device and followed by time‐lapse microscopy. During the first generations of growth, cells were grown in the absence (or with a very low concentration) of inducer and after several generations, TMG was introduced at intermediate concenteration into the medium and maintained thereafter. In the absence of TMG, cells exhibit an overall weak fluorescence yet with significant variations between cells that were shown to correspond in part to the variability in the intracellular concentration of active LacI molecules. Upon induction, transient bimodality is observed, as the cells are divided between two subpopulations of bright and dim fluorescence.

We found a strong clustering of induced cells within their genealogical trees, indicating a substantial epigenetic inheritance. This epigenetic inheritance can be traced back up to two generations prior induction, suggesting that some molecular determinants of cell fate are epigenetically inherited with a short‐range memory lasting around two divisions.

The promoter used to control fluorescence proteins expression is sensitive to small variations in active LacI concentration. Thus, in the absence of inducer, these variations result into correlated variations of YFP and CFP levels. We used the arithmetical mean of yellow and cyan fluorescence intensities to estimate the concentration of active LacI in the cells. We found that the cells exhibiting a low LacI concentration before induction are more likely to be induced upon TMG introduction. Likewise, the slowly growing cells were found to have a higher switching probability than the fast‐growing ones. We used a multivariate analysis based on a generalized linear model to estimate the correlations of pre‐induction LacI concentration and growth rate with the switching probability (Figure 5C). This analysis confirms that both LacI concentration and growth rate are correlated with the switching probability and demonstrates that even though LacI concentration and growth rate can be linked, their correlations with the switching probability represent independent effects. Together, these effects can account for 90% of the observed switching events.

To gain a better understanding of the possible influence of LacI expression fluctuations and growth rate on the switching probability of a cell, we used a model consisting in a system of differential equations and describing the dynamics of the lactose utilization network. In this model, LacI concentration controls the basal level of expression of the operon and the growth rate acts through the dilution of intracellular components. According to this model, depending on both LacI concentration and growth rate, a cell can be in a monostable or bistable regime. Therefore, monostable and bistable cells can coexist in the population due to parameters’ variability. In addition, for cells in the bistable regime, the size of the minimal LacY burst necessary to trigger induction increases with LacI concentration and growth rate. Thus, in agreement with our experimental results, these two variables control the sensitivity of the cell to permease bursts and therefore influence its switching probability.

We thus found pre‐disposing factors governing the lactose operon switching in a regime of transient bimodality. Some factors, such as LacI and LacY concentrations result from stochasticity at the local level of the network. On the contrary, growth rate variability represents variations in the cell global physiology. Therefore, the effects of local stochasticity are coupled with the influence of the global physiology, demonstrating the importance of considering the embedding of a particular genetic network in the whole cellular physiology to understand fully its dynamics.

Under conditions of bistable induction of

*Escherichia coli*lac operon, epigenetic patterns of sublineages of ‘on’ and ‘off’ cells originate from distinguishable ancestors up to two generations before induction.We found two switching pre‐disposing factors, namely low repressor levels and slow growth, demonstrating that stochasticity in gene expression and global physiology synergistically determine the single‐cell responses.

A quantitative model where growth rate acts through simple dilution of intracellular content and repressor level controls the basal activity of the operon demonstrates that both growth rate and repressor concentration influence the cell switching ability.

## Introduction

The ability of living cells to differentiate underlies the development of multicellular organisms through cellular speciation. It was shown that feedback loops often govern such differentiation processes by maintaining cells in different states (Becskei *et al*, 2001; Losick and Desplan, 2008). Such phenotype switching mechanisms were hypothesized to be crucial as well for unicellular organisms such as bacteria, to generate phenotypic diversity as a way to adapt to a randomly fluctuating environment (Thattai and van Oudenaarden, 2004; Kussell and Leibler, 2005; Dubnau and Losick, 2006; Raj and van Oudenaarden, 2008; Veening *et al*, 2008).

The initiation of differentiation processes that lead to specific gene expression patterns can be either deterministic or stochastic. As example of the former, the asymmetric division of *Caulobacter crescentus* invariantly produces a swarmer and a stalked cell with different fates and genetic expression patterns (Holtzendorff *et al*, 2004). Similarly, the distinct chemoreceptor expression levels of the two *Caenorhabditis elegans* gustatory neurons are determined by their lateral position (Johnston *et al*, 2005). In contrast, numerous systems exhibit stochastic cell fate determination, such as generation of photoreceptors in *Drosophila melanogaster* (Losick and Desplan, 2008) or competence and sporulation in *Bacillus subtilis* (Maamar *et al*, 2007; Veening *et al*, 2008). From a mechanistic perspective, the differentiation process of cells involves two levels. The local level corresponds to the functional modules responsible for the cell fate determination (e.g. intrinsic noise in gene networks) and the global level relates to the context where the differentiation is achieved (e.g. environment, physiology, extrinsic noise, etc.).

The lactose utilization network of *Escherichia coli* is a paradigm for the study of gene expression regulation and differentiation (Monod and Jacob, 1961). The genes responsible for lactose uptake and metabolism in *E. coli* (*lacZ*, *Y*, *A* coding for the lactose degrading enzyme (β‐galactosidase), permease and transacetylase, respectively) are grouped in a single transcriptional unit, the *lac* operon. The expression of this operon is regulated by the LacI repressor protein, which inhibits transcription from the P_{lac} promoter. This inhibition is relieved by inducers such as allolactose or non‐metabolizable lactose analogues, such as thiomethyl‐β‐galactoside (TMG). Early work from Novick and Weiner (1957) showed that when a clonal population of *E. coli* is treated with an intermediate concentration of TMG, two stable subpopulations emerge, with cells expressing β‐galactosidase either at a maximal or a negligible rate. This bistable behavior has been extensively studied both experimentally and theoretically ever since (Carrier and Keasling, 1999; Vilar *et al*, 2003; Ozbudak *et al*, 2004; van Hoek and Hogeweg, 2007). It was shown to be governed by the positive feedback loop in the network: TMG relieves the repression of the LacY permease, which in turn allows further TMG uptake. It is generally hypothesized that, under these conditions, the lactose operon expression state is determined by stochastic fluctuations of protein expression (Novick and Weiner, 1957; Ozbudak *et al*, 2004). In particular, induction is assumed to rely on a burst of transcription of the lac operon, leading to an increase in cellular LacY concentration (Novick and Weiner, 1957; Choi *et al*, 2008).

Here, we study phenotypic variability in growing, monoclonal populations of bacteria under conditions of transient bimodality of the lactose operon expression. By following individual cell behavior through successive generations, we emphasize the epigenetic inheritance of cell fate determinants and look for pre‐induction characteristics that would constitute predictive factors for the cell switching decision. Given the small number of LacI molecules (around 10) per cell (Müller‐Hill, 1996), its cellular concentration is likely to undergo strong cell‐to‐cell variations. These variations may affect the probability of occurrence and the size of LacY bursts as well as the size of the minimal LacY burst required for cell induction. Moreover, as the lactose network is embedded within the global cellular physiology, in itself subject to substantial cell‐to‐cell variations, we further address whether physiological parameters, such as cellular growth rate variability (Stewart *et al*, 2005; Reshes *et al*, 2008) influence the cell's switching ability. Our results, gained from coupling experimental and modeling approaches, suggest that indeed both the physiological state of the cell and its content in repressor constitute predictive factors of its response upon induction.

## Results

### Real‐time characterization of the lactose operon expression in a growing microcolony using a microfluidic device

To follow the expression profile of the *lac* operon, we used the sequenced MG1655 *E. coli* strain, modified to express the yellow fluorescent protein (YFP) and the cyan fluorescent protein (CFP) both under the control of the inducible promoter P_{L}lacO1 (strain LCY1, Figure 1A) (Lutz and Bujard, 1997), which consists of the lactose operon operator O1 and the P_{L} promoter of lambda phage (see Supplementary Information S1). In the absence of any *lac* inducer, the transcription of the fluorescence protein genes as well as of the *lac* operon is repressed by LacI. Conversely, in the presence of sufficient inducer concentration, LacI repression is inactivated, allowing expression of the fluorescence proteins (Supplementary Figure S1) and the *lac* operon genes. Thus, levels of both fluorescence markers represent active LacI concentration, and indirectly, lactose operon expression.

To observe single cells within growing microcolonies under controlled conditions that can be changed at will, we developed a microfluidic device (see Supplementary Information S2) where cells are confined to 2D growth between a glass cover slip and a thin agarose gel layer (Figure 1B; Supplementary Figure S2). Above the agarose gel, a microfluidic channel delivers the growth medium by diffusion of soluble compounds through the gel with a time scale of homogenization <1 min (see Supplementary Information S2.1). This device is simple to fabricate and to use, potentially compatible with the use of any species of bacteria as well as other cell types, and shown here to be effective in studying single‐cell responses to environmental changes.

### Real‐time follow‐up of transient bimodality in lactose operon expression

Microcolonies originating from a single cell were grown in the microfluidic device in minimal medium with succinate as a carbon source and followed by time‐lapse microscopy to obtain phase contrast and fluorescent images (Supplementary Video 1). During the initial growth phase, cells were grown either in the absence of or with a very low TMG concentration (5 μM; see Supplementary Information S3, Supplementary Table S1 and Supplementary Figure S3). After 5–6 generations, TMG (29 μM) was introduced into the flow (time *t*_{I}) and maintained until the end of the experiment. We followed the yellow fluorescence of individual cells through time (Figure 2B) and both cyan and yellow fluorescence images were taken at time *t*_{I} (Supplementary Figure S4). At this time, cells exhibited an overall weak fluorescence, yet with obvious variations between cells (Figure 2A, inset). The variations in CFP and YFP expression were highly correlated (Supplementary Figure S4). Upon TMG introduction, two types of cell responses were observed, where fluorescence intensity increased either strongly or slightly (Figure 2B and C). After two divisions, the fluorescence of most cells seems close to reaching an equilibrium (Figure 2B), of either strong or weak intensity (Figure 2A and C; Supplementary Information S1). The resulting fluorescence distribution is thus bimodal and a fluorescence threshold can be defined (see Materials and methods) to assign to each cell its most probable induction state (Figure 2C).

It is worth noting that we study a regime of transient bimodality of the population, as opposed to the more classically studied regime where the population exhibits a stationary bimodality (Ozbudak *et al*, 2004) (Supplementary Information S1). In the latter, at the population level, the switching rate from ‘off’ to ‘on’ state is compensated by the growth rate differential between ‘off’ and ‘on’ populations, cells expressing fully the operon bearing the metabolic cost of this expression. Therefore, some cells constantly switch from ‘off’ to ‘on’ state but as the ‘on’ population grows more slowly, the proportion of ‘on’ and ‘off’ cells can be stable. In our experiments, to observe a statistically relevant number of switching events (within the 2 h of induction window available in our experimental setup, limited by second layer microcolony growth and out sizing of the microscope field), we use induction conditions leading to a high switching rate. Under these conditions, at the population level, the bimodality is transient as the switching rate is too high to be compensated by growth rate differential (Supplementary Figure S1). Nevertheless, at the single‐cell level, this regime is a true bistable regime where a cell can either be in a low or high equilibrium level of expression, as demonstrated by the transient bimodality (Figure 2C; Supplementary Figure S5) and the saturations approached by individual cell fluorescence traces after 2 h of induction (Figure 2B, red and green lines).

### Epigenetic inheritance of cell fate determinants

We used the phase contrast images of the time‐lapse movies to identify each and every cell of the growing colonies (Lindner *et al*, 2008), and measured their mean yellow fluorescence level. This enabled us to reconstruct each microcolony lineage ‘history’ tree, containing both lineage information and yellow fluorescence intensity as a function of time for all individual cells (Figure 3A).

The TMG‐induced fluorescence propagates through the lineage. Clusters of highly fluorescent cells in which the *lac* operon is induced are clearly observed (red lines in Figure 3A). Quantitative analysis of this clustering is depicted in Figure 3B. For each cell present at time *t*_{I} (TMG introduction), the proportion of induced cells in its final progeny was calculated. The distribution of these proportions exhibits two peaks at 0 and 1 (Figure 3B, left), suggesting that most of the initial cells have in their progeny either almost only induced cells or almost only uninduced cells. This bimodality is not created by a specific distribution of the number of descendants per cell, as concluded from a comparison with a distribution created by random re‐sampling (Figure 3B, right; Supplementary Information S4), which is devoid of this bimodal characteristic.

The strong clustering of induced cells in the tree is indicative of a substantial epigenetic inheritance in the system (Supplementary Information S4). Once a mother cell switches and begins to express the LacY permease, its progeny will inherit the proteins and is likely to continue the switching process. But in addition to this trivial epigenetic inheritance, the molecular determinants of the response might also be inherited, leading to correlated switching decisions of non‐induced related cells. In support of this is the fact that the first common ancestor of clusters often divides before introduction of TMG (Figure 3A; Supplementary Figure S3, black circles). To test the epigenetic inheritance of cell fate determinants, we compared the responses of related cells before induction (Figure 4). We found that sister cells whose mother divided before TMG introduction exhibit a correlated behavior as the mean fluorescence of the progenies are correlated (Spearman's rank correlation ρ=0.61, *P*‐value=10^{−7}; Figure 4B). Therefore, some determinants of cell fate are epigenetically inherited and lead to correlated switching decisions of sister cells. To see how many generations this epigenetic memory can last, we also compared the responses of cells sharing a more distant common ancestor. To avoid statistical biases created by using a cell several times (as it has for instance several cousins), we compared the mean fluorescence of the progeny of two sister cells drawn from the populations of mothers or grandmothers of the cells present at TMG introduction (Figure 4A). We found that the first cousins still exhibit a correlated response (Spearman's rank correlation ρ=0.57, *P*‐value=10^{−7}; Figure 4B), whereas this correlation is lost among second cousins (Spearman's rank correlation ρ=0.02, *P*‐value=0.9; Figure 4B, inset). Thus, it seems that there is a short‐range memory in the determinants of cell fate that lasts around two divisions.

### LacI concentration influences switching probability

Stochasticity in gene expression can be categorized into two broad classes: extrinsic and intrinsic (Elowitz *et al*, 2002). Extrinsic noise corresponds to fluctuations in the concentrations of polymerases, ribosomes or regulatory factors that would cause variations of expression of a given gene between one cell and another but not between two identical genes in the same cell. Intrinsic noise corresponds to stochasticity in the molecular events involved in the transcription and translation of genes and thus governs variation between the expression of two identical genes in the same cell.

We estimated LacI concentration indirectly, from the expression levels of the YFP and CFP genes. Thus, intrinsic noise in expression of these genes could substantially reduce the accuracy of our estimation. To minimize the effect of intrinsic noise, we introduced a variable corresponding to the arithmetical mean of the normalized yellow and blue fluorescence at time *t*_{I} of TMG introduction (*F*_{I}=1/2(YFP/〈YFP〉+CFP/〈CFP〉; Supplementary Figure S4). Therefore, *F*_{I} estimates the correlated variations of YFP and CFP levels. To test whether *F*_{I} is a predictive factor of the subsequent switching behavior, we divided the ancestors at *t*_{I} into two groups (Figure 5A), according to the proportion of induced cells in their progeny. Ancestors having more than half of their progeny induced have an expression level 50% higher (*t*‐test *P*‐value<0.001) than the other ancestors. Thus, cells exhibiting a higher value of *F*_{I} have a higher probability to switch.

The correlated variations in YFP and CFP levels have two potential origins. It can derive either from a variation in the global protein synthesis ability of the cell, based for example on a variation in RNA polymerase or ribosome concentrations, or from a variation in LacI concentration. Indeed, even a small variation in active LacI concentration is likely to have an impact on the fluorescence level controlled by P_{L}lacO1, as this level is continuously increasing when TMG concentration increases between 0 and 12 μM (Supplementary Figure S1), concentrations corresponding to the monostable non‐induced regime. Therefore, this promoter is very sensitive to variations in active LacI concentration. To evaluate the relative importance of the two sources of fluctuations of *F*_{I}, global and LacI dependent, we repeated the same induction experiment with the strain LCY2, carrying the CFP gene under the control of a constitutive promoter (P_{R} lambda promoter) and the YFP gene under the control of P_{L}lacO1 (as before). As CFP expression is constitutive, its extrinsic variations stem only from global sources, whereas YFP expression is in addition directly sensitive to the level of LacI in the cell. In agreement, in the absence of induction, YFP expression is much more heterogeneous than CFP expression in LCY2 (Figure 5B; Supplementary Figure S7A): the noise in YFP level, defined as the ratio of the standard deviation over the mean is 0.31, whereas the noise in CFP level is only 0.14. Conversely, when fully induced, YFP expression exhibits the same noise as CFP expression (0.14), and YFP and CFP expressions are highly correlated (Supplementary Figure S7B; Spearman's rank correlation ρ=0.8, *P*‐value<10^{−10}). This strong correlation demonstrates that intrinsic noise in CFP expression is limited and that CFP fluorescence is therefore a relevant measure of global variations in protein synthesis ability. The correlation is partially lost under repressed conditions (Supplementary Figure S7A; Spearman's rank correlation ρ=0.4, *P*‐value<10^{−10}), due to the additional noise caused in YFP expression by LacI repression. Therefore, under repressed conditions, a substantial amount of the variations in YFP expression is caused by LacI fluctuations. We used a generalized linear model to analyze the respective influence of CFP and YFP initial levels on the switching probability of a cell (see Materials and methods for a description of the analysis). We found that the initial YFP level was a better predictor of the response of the progeny than initial CFP level. In the generalized linear model, the coefficient for the corresponding reduced centered variable was higher for YFP (2.9) and statistically relevant (*P*‐value<10^{−5}) as compared with the CFP fluorescence (coefficient: 0.9; *P*‐value>0.05). From this analysis, we conclude that LacI cellular concentration is influencing the probability of lac operon switching.

### LacI concentration can act on switching probability through the modulation of LacY transcription burst sizes

LacI concentration at induction can directly affect the switching probability, by determining the level of inhibition of the *lac* operon and thus the size of the minimal LacY burst that would trigger induction. Alternatively, it can have an indirect effect through the determination of LacY initial concentration: lower cellular levels of LacI might result in more frequent or larger bursts of LacY expression. In our system, the fluorescent reporters monitor only the dynamics of the LacI‐dependent transcription initiation with no direct information on the expression level of LacY. Moreover, the synthetic promoter used to control the fluorescence expression is leakier than the one of the *lac* operon. Therefore, to test whether high levels of fluorescence before induction correspond also to high levels of LacY, we assessed independently its concentration. Although direct LacY level detection is not easily performed (Choi *et al*, 2008), quantification of the β‐galactosidase enzyme encoded by the *lacZ* gene within the same operon is straightforward and should reflect the adjacent *lacY* expression level. Cells grown without inducer were sorted according to their yellow fluorescence into two subpopulations, differing in their mean expression levels by a factor of 1.45 (see Materials and methods). The β‐galactosidase activity of the population exhibiting the highest mean fluorescence was 50% higher than that of the other population. A lower cellular level of LacI is therefore correlated with a higher level of lactose operon genes expression, as it may lead to either more frequent or larger transcription bursts. Thus, the effect of LacI stochasticity on the *lac* operon switching probability may be mediated at least in part through its control of LacY burst size or frequency. From the molecular perspective, it seems rather unlikely that LacI concentration would affect the frequency of bursts, as a burst result from a unimolecular event of dissociation of the repressor from its operator site. On the contrary, the possible effect on burst size seems clear, as the LacI association to its operator is concentration dependent. Indeed, it was shown in Choi *et al* (2008) that the size of the bursts directly depends on the TMG concentration in the medium, through its effect on intracellular LacI activity. It is thus likely that when the active LacI concentration is lower in the cell, once the repressor is dissociated from the operators, it takes a longer time for another LacI molecule to bind the operator, which results in a larger transcription burst.

### The lactose utilization network within the physiological context: growth rate influences cell response

To calculate the growth rate of a cell in the microscopy experiments, we measured the length of the cell at every time points using the image analysis software developed by Stewart *et al* (2005). As the diameter of the cell is constant during growth, the growth rate is given by an exponential fit to the length over time (see Materials and methods). To test whether growth rate has an impact on the switching decision of a cell, we retrospectively divided the ancestors at time *t*_{I} into two groups, with respect to the proportion of induced cells in their progeny. The initial growth rate is 15% smaller in the subpopulation of ancestors that have more than half of their final progeny induced when compared with their complementary subpopulation (*t*‐test *P*‐value<0.01). This difference is significantly larger (*P*‐value=0.05) than the difference in growth rate caused by the cost of fully expressing the operon's genes (6% difference between uninduced and fully induced cells under our conditions, in agreement with previously published data (Dekel and Alon, 2005)). Thus, there is a substantial correlation between the growth rate of a cell and its decision to switch: slower growing cells have a higher propensity to switch.

### Low LacI concentration and slow growth independently increase switching probability

We experimentally found two different variables correlating with the switching probability: cells having a higher initial fluorescence *F*_{I} or exhibiting slower growth had a higher probability to switch and thus a higher proportion of induced cells in their progeny. Nevertheless, growth rate and *F*_{I} could be linked, for example if the growth rate influences the intracellular concentration of LacI. In this case, the effect of growth rate could be indirect and redundant with the effect of LacI concentration (see Supplementary Information S5 for further discussion). To quantify the proper effect of each variable, we used a generalized linear model (see Materials and methods). This statistical analysis has several major advantages over a more classical correlation study. First, it allows to disentangle the proper effect of each explanatory variable when these variables are correlated. Indeed, the model tests if when one variable is known, the knowledge of the second one brings additional information. The effect of each variable, independently of the other, can thus be quantified. Second, the genealogical structure of the data is taken into account in the model thus avoiding the biases that would be created in a simple correlation study by the fact that all ancestor cells do not have the same number of descendants.

The linear model confirms that the proportion of induced cells in the progeny is positively correlated to *F*_{I} (*P*‐value<10^{−5}; Supplementary Table S1A) and independently negatively correlated to the growth rate (*P*‐value<10^{−4;} Supplementary Table S1A) and allows an estimation of the switching probability as a function of these two variables (Figure 5C). Our data shows that cells having more than half of their progeny induced are clustered in a region of low LacI levels (high *F*_{I}) and small growth rates. Conversely, cells having less than half of their progeny induced are clustered in a region of high LacI levels and high growth rates (Figure 5C). Likewise, the switching probability estimated with the generalized linear model is a decreasing function of initial LacI concentration and growth rate. According to this estimation, for a given growth rate, initial LacI content can vary the switching probability from 0 to 1. Similarly, for a cell with an initial fluorescence *F*_{I} of 1, switching probability would vary from 60% for a low growth rate of 0.005 min^{−1} to 1% for a high growth rate of 0.02 min^{−1}. Using the model, the initial fluorescence and growth rate of a cell allow us to predict if its progeny will be as a majority induced or non‐induced with a 93% accuracy. In all, 99% of the cells that have more than half of their progeny induced are predicted to be in this category and 60% of the others are correctly predicted to have less than half of their progeny induced. Thus, the initial LacI concentration and growth rate are independent predictive factors controlling the cell response upon induction.

### A model of single‐cell dynamics demonstrates the effect of growth rate and LacI concentration on a cell's sensitivity to LacY permease transcription bursts

To investigate how the growth rate could directly influence the decision of a cell to switch, we adapted the model of Ozbudak *et al* (2004) to describe the lactose utilization network at the single‐cell level. As discussed above, LacI concentration may act on switching probability through the modulation of LacY bursts size. However, LacI concentration could also have a direct effect in changing the cell's sensitivity to a given LacY burst. Our modeling approach also allowed us to investigate whether the data are quantitatively compatible with such a direct influence of LacI concentration.

The model consists in a system of three differential equations governing TMG intracellular concentration and LacY and YFP cellular concentrations (see Supplementary Information S6). The third equation determining YFP concentration is used only for parameter estimation and we therefore consider here only the two equations determining TMG intracellular concentration (*x*; in dimensionless unit) and LacY concentration (*y*):

TMG enters the cell at a rate proportional to the permease concentration (with a rate of TMG uptake per LacY molecule β) and is depleted by dilution due to cell growth (van Hoek and Hogeweg, 2007) (with a cellular growth rate μ). Likewise, LacY is diluted, and degraded in a first‐order reaction with a rate γ (van Hoek and Hogeweg, 2007), and its production is an increasing function of TMG intracellular concentration (see Supplementary Information S6 for further information), with maximal value α and minimal value α/ρ. The repression factor ρ is a linear function of the cellular LacI concentration.

To model the cell‐to‐cell variations before induction, we consider distributions in the three factors that we experimentally determined to be involved in the switching decision. Individual growth rate and LacI concentration are considered as parameters (through μ and ρ; see Supplementary Information S6), whereas LacY initial variation is modeled as a variation in initial conditions. The distributions were derived from the microscopy experiments, and the other parameters of the system were measured or calculated from additional experiments (see Supplementary Information S6 for parameter calculation; Supplementary Table S2 for parameter estimation). The model does not include any adjustable parameter, all parameters being measured or calculated from experiments.

To investigate the influence of LacI concentration and growth rate on the switching probability, we built bifurcation diagrams, which show the possible equilibria of the system as a function of LacI concentration (Figure 6A) or growth rate (Figure 6B). In Figure 6A, the bifurcation parameter is the repression factor ρ, which is linearly related to LacI concentration, and the growth rate is fixed. When ρ is smaller than a minimal threshold (left dashed black line), there is only one stable fixed point (blue line). Thus, there is only one possible equilibrium for the cells, corresponding to an induced state that is reached regardless of LacY initial concentration. Likewise, if ρ is higher than a maximal threshold (right dashed black line), the system exhibits only one stable fixed point, corresponding to the non‐induced state. When ρ is between these two threshold values, the system exhibits two stable fixed points and is driven by initial LacY concentration. Thus, for a given ρ in this range, if the initial LacY concentration is high (above the black line), the cell will reach the induced state, whereas if the initial LacY concentration is low (below the black line), the cell will stay in the non‐induced state. Therefore, a cell can be in a monostable regime induced or non‐induced or in a bistable regime depending on its intracellular LacI concentration. Furthermore, in the bistable regime, the minimal LacY initial concentration necessary to induce the cell (black line) increases with LacI concentration. The range of ρ calculated from our experimentsis between 100 and 750. For this range of ρ values, the cells can be either in the monostable‐induced regime or in the bistable regime.

Likewise, the effect of growth rate is shown in Figure 6B. It is qualitatively similar to the effect of LacI concentration: depending on its growth rate, a cell can be in a monostable regime induced or non‐induced or in the bistable regime. Moreover, in the bistable regime, the minimal LacY initial concentration necessary to induce the cell (black line) increases with growth rate. The cellular growth rates observed in our experiments range between 0.010 and 0.017 min^{−1}. For this range of values, the cells can be either in the monostable‐induced regime or in the bistable regime.

In our model, the protein production rate α is constant and does not change with growth rate. What determines the growth rate at the single‐cell level is unknown, but its fluctuations may be correlated to fluctuations in protein production rate. In particular, it is known that at the population level, ribosome and RNA polymerase concentrations change with growth rate (Nomura *et al*, 1984; Bremer and Dennis, 1987; Klumpp and Hwa, 2008). This property might be valid at the single‐cell level and consequently, we also did a bifurcation analysis in which the parameter α was proportional to growth rate (α=α_{0} × μ). The effect of growth rate on switching probability is also valid in this case (see Supplementary Figure S8A).

In our model, the growth rate acts through the dilution of intracellular TMG and LacY permease, but we do not include such a dilution effect for LacI. LacI is negatively autoregulated, as the protein can bind to an operator site located at the end of its own gene, leading to truncated transcripts. Because of this autoregulation, the mathematical relationship between growth rate and LacI concentration is unknown. Nevertheless, when simple LacI dilution is modeled without taking the autoregulation into account, the negative influence of growth rate on switching probability still holds (see Supplementary Figure S8B). In addition, negative autoregulation has been shown to provide stability and increased robustness to fluctuations in biochemical parameters (Becskei and Serrano, 2000). It is thus expected to decrease the sensitivity of LacI concentration with respect to fluctuations in growth rates, and therefore the dilution effect.

Thus, according to our model, both LacI concentration and growth rate have a strong influence on the size of the smallest LacY burst required to trigger induction. Therefore, they both control the sensitivity of a cell to a given LacY burst.

## Discussion

We designed a microfluidic setup (Figure 1B) allowing to follow single cells growing in microcolonies under a dynamically controlled environment. We used this device to study the formation of phenotypic variability in bacterial microcolonies upon induction of the lactose operon in a bistable regime. The genealogical context of the cells indicates that closely related cells show a correlated behavior upon induction (Figure 3). In particular, the responses of sister cells and first cousin cells are substantially correlated (Figure 4). As this correlation is lost between second cousin cells, the epigenetic memory of cell fate determinants seem to be on the order of two generations. A similar phenomenon of epigenetic inheritance was found in a stochastic switch between two epigenetic states in the yeast *Saccharomyces cerevisiae* (Kaufmann *et al*, 2007), where it was suggested to be the consequence of bursts in expression of a regulatory protein. In the model developed by Kaufmann *et al* (2007), for the same average level of the regulatory protein, a higher burst size leads to a stronger correlation between related cells. After division of a cell, as long as no burst has occurred in any of the daughter cells, the behavior of these cells should stay highly correlated. As the burst size increases, the period between bursts increases, leading to longer periods of correlated behavior between related cells. Similarly, bursts in LacY or LacI expression could both explain the epigenetically inherited switching pre‐disposition in our system. Nevertheless, correlated behavior between related cells could also emerge from inheritance of physiological parameters. As growth rate has an influence on the cell's response and growth rate between related cells are correlated (for sister cells, Spearman's rank correlation ρ=0.7, *P*‐value<10^{−10}), the observed inheritance could as well be due to an inheritance of the global physiological state revealed by cellular growth rate.

We found pre‐disposing factors of the lactose operon switching, under conditions of transient bimodality of the population (Figure 5). Pre‐dispositions in cells response have already been studied in eukaryotic systems (Cohen *et al*, 2008) as well as in prokaryotic systems, in particular in two other model systems of bistability, the lambda phage lysis/lysogeny decision and the sporulation process in *B. subtilis*. In *B. subtilis* sporulation, epigenetic inheritance lasting up to two generations was found but no pre‐disposing factors could be determined (Veening *et al*, 2008). Interestingly, in the lambda lysis/lysogeny decision, it has been shown that the size of the cell before infection was a strong predictor of cell fate (St‐Pierre and Endy, 2008).

Together, our experimental and theoretical approaches suggest that LacI concentration acts both directly on LacY transcription bursts size and indirectly on the cell's sensitivity to these bursts in presence of TMG. LacY transcription bursts are triggered by the dissociation of the repressor molecule from the operator site of the promoter. Once this dissociation has happened, transcription can take place until a new LacI molecule binds to the operator. If the active LacI concentration is high, this binding of a new molecule is likely to happen sooner than for a low repressor concentration. Thus, the transcription period and therefore the size of the burst are expected to be negatively correlated to the cell's LacI content. A cell with a low LacI content is therefore expected to have a higher LacY content independently of the presence of inducer in the medium. But the repressor concentration could also act on the sensitivity to a given LacY concentration in the presence of TMG. Indeed, once a LacY burst occurs, some TMG molecules enter the cell. If the number of active LacI molecules inside the cell is low, their depletion through their binding to TMG molecules is expected to be more dramatic. In this case, active LacI depletion could in turn lead to increased LacY burst size or even frequency if a free TMG molecule manage to bind to the repressor bound to DNA and cause its dissociation from DNA.

Variability in LacI and LacY concentrations can emerge from stochasticity in genes expression at the local level of the lactose utilization network. We also found that phenotypic variability at the more global level of the cell's physiology has a function in a cell's decision to switch on the lactose operon. Indeed, the growth rate of a cell is a predictor of its subsequent response upon induction, slowly growing cells being more readily induced than faster growing cells (Figure 5). According to our model, the fact that a fast‐growing cell dilutes its molecular content, in particular its intracellular inducer, faster than a slow‐growing cell could be sufficient to explain the negative correlation between growth rate and switching probability. Thus, based on experiments and modeling, we found that local stochasticity in gene expression at the network level and global physiology synergistically determine the single‐cell responses. Therefore, both our results and the results of St‐Pierre and Endy (2008) show that to fully understand the dynamics and properties of a particular network, it is of great importance to consider its embedding in the whole cellular physiology.

The growth rate of a bacterium directly reflects its fitness and adaptation to the environment. Thus, the lactose operon may offer an example of a fitness‐dependent switching mechanism, where cells that have a slower growth rate and consequently a lower fitness are more susceptible to switching. A mechanism of fitness‐dependent epigenetic state determination was recently proposed and studied in a synthetic bistable switch (Kashiwagi *et al*, 2006). In that study, it was suggested that without any signal transduction and sensing machineries, a cell could adapt to its environment by a mechanism named fitness‐induced attractor selection. According to this model, when the cell is in a non‐adaptive state, the metabolism is lowered and the effect of stochasticity in gene expression is increased, which renders the non‐adaptive state unstable. However, in our data, there is no difference in expression variability between cells of smaller versus bigger growth rates (see Supplementary Information S7). Contrary to this mechanism, we propose that in the lactose operon bistable switch, the effect of growth rate on cell response to the inducer is based on the network architecture and on the dilution rates of the molecules involved. The lactose utilization network includes a sensing and signaling machinery. However, in conditions of bistability, the epigenetic state of the cell is no longer determined by the external conditions, and the switching probability seems to be influenced by growth rate, demonstrating a fitness‐dependent epigenetic determination. A mechanism where the slowly growing, less fit cells could more readily change their genetic expression, thus exploring the space of epigenetic states to optimize their fitness, could promote the adaptation of genetic expression to the environment. Growth rate would thus be used by the cell as a measure of environmental adaptation quality. In particular, in the case of the lactose operon, a slow‐growth rate may indicate a poor nutrient utilization capacity, hence the benefit of expressing alternative nutrient utilization systems.

The major difference between the transient bimodality under study and the stationary bimodality classically studied is at the population level and not at the single‐cell level (see above). Indeed at the single‐cell level, both regimes are characterized by the existence of two equilibrium states, ‘on’ and ‘off’, and a transition probability from ‘off’ to ‘on’ state. Nevertheless, as we use a high inducer concentration leading to a high transition probability, the mechanisms of switching could be different from those acting in the low switching rate conditions (lower inducer concentration). According to our model, the effects of LacI concentration and growth rate on the switching probability should not depend on inducer concentration and therefore can reasonably be expected to concern also the common regime of stationary bimodality. In this regime, some rare but important events like errors in LacI transcription have been shown to trigger induction (Gordon *et al*, 2009). This phenomenon can be seen as an extreme case of the effect of LacI concentration on the switching probability.

The different regimes under which the bimodality of the population is stationary or only transient are therefore both regimes of bistability (meaning existence of two equilibrium states), with low or high switching rates. The transiently bimodal regime could be of greater biological relevance, as it is leading to a highly variable response on a short‐time scale. With a lower switching rate, the variability can be maintained stably in the population but its formation requires that the cells grow exponentially in a constant environment on a long time scale (Ozbudak *et al*, 2004), conditions that *E. coli* is unlikely to meet in its natural environment.

## Materials and methods

### Bacterial strains and culture conditions

The two strains used in this study are constructed from the M22 (MC4100 galK∷P_{L}lacO1‐CFP intC∷P_{L}lacO1‐YFP) and MRR (MC4100 galK∷λPR‐CFP intC∷λPR ‐YFP) strains (kind gift from Elowitz *et al* (2002)) by P1 transduction into MG1655. The strain LCY1 is MG1655 galK∷P_{L}lacO1‐CFP intC∷P_{L}lacO1‐YFP and the strain LCY2 is MG1655 galK∷λPR‐CFP intC∷P_{L}lacO1‐YFP. intC and galK are located on the chromosome symmetrically from the origin of replication (Elowitz *et al*, 2002).

Cells were grown in M9 minimal medium supplemented with succinate (0.4%), casamino acids (0.01%), MgSO_{4} (3 mM), thiamine (30 mM), and uracil (10 mM), at 30°C with agitation. For microscopy experiments, overnight cultures were first diluted 1000‐fold into 5 ml of fresh medium. Four hours later, the exponentially growing cells were diluted 20‐fold and inoculated into the microfluidic device.

### Microfluidic device fabrication and use

Soft lithography was used to fabricate the mould for the channels of 300 μm width and 110 μm height cast in PDMS (see Supplementary Information S2 for fabrication).

The agarose gel (1.5 wt%) was prepared as follows: a desired weight of agarose powder (Qbiogene, QA‐Agarose TM) was added to the succinate Minimal Medium whose exact composition is described above. The mixture was boiled to dissolve the agarose, and allowed to cool to 45°C. A drop of this solution was placed in a mould (height 130 μm) preheated at 45°C and covered with a glass slide.

The microfluidic device was assembled as follows: after cooling down to room temperature, the agarose slab was carefully peeled off from the mould, and 3 μl of bacterial solution was deposited on the top of the agarose slab. After total liquid absorption, the agarose slab was placed on a cover slide (freshly treated with air plasma) to confine bacteria between the activated glass slide and the agarose slab. The freshly plasma‐activated PDMS device was then bound to the glass slide to get the microfluidic channel in the central part of the agarose slab. After 15 min at 30°C, the bond between PDMS and glass was strong enough to support the pressure of the flow through the microchannel.

### Microscopy and image analysis

The microfluidic setup was placed in a temperature controlled (Life Imaging Services GmbH Cube and Box) automated inverted microscope (Zeiss Axiovert 200 M) at 30°C. Up to four fields located under the microfluidic channel and containing a single cell were identified and followed by automated time‐lapse acquisition (MetaMorph microscope control software, Universal Imaging), as described earlier (Stewart *et al*, 2005; Lindner *et al*, 2008), with minor modifications (see Supplementary Information S3.2).

The custom analysis software (BHV) (Stewart *et al*, 2005) was used to identify and follow the cells and reconstruct their lineage using the phase contrast images (Lindner *et al*, 2008) and measure their fluorescence (average of pixel intensity) through the induction experiment. Growth rates were determined by an exponential fit to the cells’ length over time. Calculation of growth rate at induction was limited to not more than 10 min after TMG introduction (time window in which the increase in fluorescence is not detectable) to avoid bias from lower growth rate of induced cells.

### Statistical analysis of the influence of chemical noise and growth rate on switching probability

The final state of each cell (induced or non‐induced) was assessed from its yellow fluorescence intensity 2 h after TMG introduction. Using the ‘mclust’ function of the R open source software, the final fluorescence distribution was fitted with a mixture of two Gaussian distributions, thus allowing the definition of a fluorescence threshold above which a cell should be considered as induced.

We used a generalized linear model with a binomial distribution and a logistic link to analyze the effect of different variables on the cell switching probability. The number of cells induced in the progeny of an ancestor cell is modeled as a binomial variable with parameters *n*, the number of cells in the progeny, and *p* given by: *p*=(e^{ax+by+c})/(1+e^{ax+by+c}) with *x*=(*X*−〈*X*〉)/*sd*(*X*) and *y*=(*Y*−〈*Y*〉)/*sd*(*Y*), *X* and *Y* being the explicative variables. The use of the reduced centered variables *x* and *y* allows comparing the respective influence of the different explicative variables by comparing *a* and *b*.

### Cell sorting and β‐galactosidase assay

Bacteria were grown overnight at 30°C then diluted 100‐fold and further incubated for 3 h at 30°C. Cells were then placed at 4°C and sorted in two subpopulations with a cell sorter (FACS‐Aria with Diva software (Becton Dickinson)) according to their yellow fluorescence intensity (the gate was set at the median fluorescence of the whole population). After sorting, the two populations were composed of 20 000 000 cells each and showed a 45% difference in fluorescence. A β‐galactosidase assay was performed on each sorted subpopulation, as described earlier (Miller, 1992). The activity was calculated from the kinetics of *o*‐nitrophenol production, based on four time points.

## Acknowledgements

We thank Didier Chatenay and Stéphane Douady for advice and helpful discussions; Simon Baeriswyl and Médéric Diard for expertise in molecular biology; Marie‐Florence Bredèche and Marjorie Selva for technical assistance; Eric Stewart, Marina Elez, Xavier Manière, and Branimir Lukic for critical reading of the paper. This work was supported by an HFSP, The French National Research Agency (ANR), the Centre National de la Recherche Scientifique (CNRS) grants, and EURYI award (FT).

*Author contributions:* LR, YC, FT, DB, and ABL designed the research. LR conceived and performed the experiments, data analysis, and modeling. LR and GP analyzed the model. LR wrote the article with contributions from GP, DB, and ABL.

## Conflict of Interest

The authors declare that they have no conflict of interest.

## Supplementary Information

Video 1. LCY1 microcolony growth [msb201012-sup-0001.mov]

Supplementary Information

Supplementary text, Supplementary Tables S1 and 2, Supplementary Figures S1–S11, Video 1 Legends [msb201012-sup-0002.pdf]

## References

This is an open‐access article distributed under the terms of the Creative Commons Attribution License, which permits distribution, and reproduction in any medium, provided the original author and source are credited. This license does not permit commercial exploitation without specific permission.

- Copyright © 2010 EMBO and Macmillan Publishers Limited