Location via proxy:   [ UP ]  
[Report a bug]   [Manage cookies]                
Advertisement

  • Loading metrics

Open Access

Peer-reviewed

Research Article

Noise Expands the Response Range of the Bacillus subtilis Competence Circuit

  • Andrew Mugler ,

    Contributed equally to this work with: Andrew Mugler, Mark Kittisopikul

    Affiliation Department of Physics and Astronomy, Purdue University, West Lafayette, Indiana, United States of America

    ⨯
  • Mark Kittisopikul ,

    Contributed equally to this work with: Andrew Mugler, Mark Kittisopikul

    Affiliations Department of Biophysics, University of Texas Southwestern Medical Center, Dallas, Texas, United States of America, Division of Biological Sciences, University of California San Diego, California, United States of America

    ⨯
  • Luke Hayden,

    Affiliations Department of Physics and Astronomy, Purdue University, West Lafayette, Indiana, United States of America, Division of Natural Sciences, Indiana Wesleyan University, Marion, Indiana, United States of America

    ⨯
  • Jintao Liu,

    Affiliation Division of Biological Sciences, University of California San Diego, California, United States of America

    ⨯
  • Chris H. Wiggins,

    Affiliation Department of Applied Physics and Applied Mathematics, Columbia University, New York, New York, United States of America

    ⨯
  • Gürol M. Süel ,

    gsuel@ucsd.edu (GMS); awalczak@lpt.ens.fr (AMW)

    Affiliation Division of Biological Sciences, University of California San Diego, California, United States of America

    ⨯
  • Aleksandra M. Walczak

    gsuel@ucsd.edu (GMS); awalczak@lpt.ens.fr (AMW)

    Affiliation Laboratoire de Physique Théorique, CNRS, Université Pierre et Marie Curie and École Normale Supérieure, Paris, France

    ⨯
  • Article
  • Authors
  • Metrics
  • Comments
  • Media Coverage

Abstract

Gene regulatory circuits must contend with intrinsic noise that arises due to finite numbers of proteins. While some circuits act to reduce this noise, others appear to exploit it. A striking example is the competence circuit in Bacillus subtilis, which exhibits much larger noise in the duration of its competence events than a synthetically constructed analog that performs the same function. Here, using stochastic modeling and fluorescence microscopy, we show that this larger noise allows cells to exit terminal phenotypic states, which expands the range of stress levels to which cells are responsive and leads to phenotypic heterogeneity at the population level. This is an important example of how noise confers a functional benefit in a genetic decision-making circuit.

Author Summary

Fluctuations, or “noise”, in the response of a system is usually thought to be harmful. However, it is becoming increasingly clear that in single-celled organisms, noise can sometimes help cells survive. This is because noise can enhance the diversity of responses within a cell population. In this study, we identify a novel benefit of noise in the competence response of a population of Bacillus subtilis bacteria, where competence is the ability of bacteria to take in DNA from their environment when under stress. We use computational modeling and experiments to show that noise increases the range of stress levels for which these bacteria exhibit a highly dynamic response, meaning that they are neither unresponsive, nor permanently in the competent state. Since a dynamic response is thought to be optimal for survival, this study suggests that noise is exploited to increase the fitness of the bacterial population.

Citation: Mugler A, Kittisopikul M, Hayden L, Liu J, Wiggins CH, Süel GM, et al. (2016) Noise Expands the Response Range of the Bacillus subtilis Competence Circuit. PLoS Comput Biol 12(3): e1004793. https://doi.org/10.1371/journal.pcbi.1004793

Editor: Sergei Maslov, University of Illinois at Urbana-Champaign, UNITED STATES

Received: August 30, 2015; Accepted: February 5, 2016; Published: March 22, 2016

Copyright: © 2016 Mugler et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Data Availability: Data and code available at http://github.com/amugler/noise-bacillus.

Funding: MK was supported by NIH biophysics training grant T32GM008297. LH was supported by NSF REU grant PHY-1460899. GMS was supported by the National Institute of General Medical Sciences Grant R01 GM088428 and the National Science Foundation Grant MCB-1450867. Work in Paris was supported by MCCIG grant no. 303561. This work was also supported by the San Diego Center for Systems Biology (NIH Grant P50 GM085764). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing interests: The authors have declared that no competing interests exist.

Introduction

Snapshots of bacterial populations often reveal large phenotypic heterogeneity in the gene expression states of its composite individuals. Such phenotypic heterogeneity in a clonal population of bacterial cells in a single environment has significant consequences for how well the organisms can adapt and survive. On the one hand, a population with little or no heterogeneity may allow for all cells to take advantage of certain optimal conditions to which the population is exposed. In this case, heterogeneity is suboptimal and therefore detrimental to fitness. On the other hand, numerous recent studies have shown that heterogeneous populations allow for cells to account for uncertainty in future environmental conditions [1–6]. In this case, heterogeneity is beneficial to fitness. A straightforward way to maintain high phenotypic heterogeneity is for each cell to exhibit a dynamic response. This allows each cell in its turn to transition among the various states of the population, e.g. via switching, pulsing, or oscillatory dynamics. The heterogeneity is intrinsically encoded in each cell, and is often enhanced by, or even entirely due to, stochasticity, or “noise”, at the molecular level [2–4, 6].

The ability of molecular noise to cause stochastic phenotype changes has been demonstrated in a number of biological systems. In the context of enzymes, several studies have explored how intrinsic noise due to low numbers of molecules, or even a single molecule, can have dramatic effects through the amplified actions of a few enzymes [7, 8]. Moreover, studies of bacterial operons, including in the context of bacterial persistence, have suggested that stochasticity could be encoded in the interactions between genes in a genetic regulatory network by ensuring that certain operon states are exposed to low numbers of molecules [9–11]. Recently, a theoretical study has demonstrated the conditions for when deterministic approaches to modeling genetic circuit dynamics break down, due to amplified effects of rare events caused by a small number of regulators [12]. Together, these works suggest that phenotypic heterogeneity could be rooted in low-molecule-number noise, and that this noise could in turn be encoded in the architecture of genetic regulatory networks.

The competence response of the gram-positive bacterium Bacillus subtilis provides a striking example of dynamically maintained phenotypic heterogeneity. Under stress, B. subtilis undergoes a natural and transient differentiation event, termed competence, that allows the organism to incorporate exogenous genes into its genome. Previous studies have shown that entry into the competent state is controlled by a genetic circuit that that can be tuned to one of three dynamical regimes [13]: an excitable regime at low stress levels, where cells rarely and transiently enter the competent state; an oscillatory regime at intermediate stress, where cells oscillate in and out of the competent state; and a mono-stable regime at high stress, where cells remain in the competent state. Importantly, oscillatory (and repeatably excitable) dynamics lead to phenotypic heterogeneity, since cells are dynamically transitioning in and out of the competent state (see Fig 1A). This heterogeneity is especially important to the survival of B. subtilis: if no cells respond, competence is not exploited, and the population may succumb to the stress. On the other hand, if all cells are permanently in the competent state, this can also be fatal to the population, since competence has been shown to reduce the cell growth rate and prevent cell division due to the inhibition of FtsZ [14, 15]. Therefore, maintaining a dynamic competence response, and therefore a heterogenous population, is thought to be crucial to survival under stress.

thumbnail
Download:
Fig 1. Schematic illustrating phenotypic heterogeneity and the effects of noise.

(A) When all cells in a population exhibit either no response (left) or a high response (right), then the population is homogenous. In contrast, if individual cells exhibit a dynamic response (middle), this leads to a heterogenous population, with a fraction f of cells in the responsive state at any given time. (B) Intrinsic noise affects the dynamics of the response. Without intrinsic noise, the viable stress level response range is narrow as indicated by the black dashed lines since it is limited to deterministic dynamics. However, for the B. subtilis competence response, we find in this study that noise expands the viable response range: the range of stress levels over which f remains neither 0 nor 1. This expansion is illustrated by the blue dashed lines indicating the extent where f, solid blue line, remains between 0 and 1. f permits a heterogeneous population when between 0 and 1.

https://doi.org/10.1371/journal.pcbi.1004793.g001

The competence circuit includes a negative feedback loop, which is known to play a critical role in controlling the exit from competence [15], the duration of competence [13], the variability of the duration of competence [16], and the integration of signals into competence that permit the existence of a transient competent state [17]. However, the effects of noise on the dynamics of the competence response are only partially understood. Previous work has shown that noise can trigger excitations into the competent state when the circuit is tuned to the excitable regime [15]. This transition into the competent state was shown to be governed by intrinsic noise rather than extrinsic noise [18]. Later work showed further that these excitations have a high variability in their duration, and that this variability is directly linked to the architecture of the competence circuit [16]. In particular, this work employed an analogous synthetic excitable circuit, termed SynEx, to provide evidence that the duration variability is due to intrinsic noise from low molecule numbers in the native circuit. Additionally, a similar strain, SynExSlow, was created that had a mean competence duration comparable to the native strain at the cost of additional complexity in order to reduce the efficiency of a proteolytic negative feedback loop [16]. However, the ability of this intrinsic noise to trigger sustained or repeatable excitations has not yet been quantified. Moreover, the generic effects of intrinsic noise on the three dynamic regimes, and how these effects translate to the physiological function of B. subtilis at the population level, are unknown.

Here, using stochastic modeling and quantitative fluorescence microscopy, we study the effects of intrinsic noise on the competence dynamics and the ensuing population heterogeneity of B. subtilis. We uncover a novel effect of noise that goes beyond architecture-dependent stochastic effects in a single cell. Specifically, we find that at both low and high stress levels, noise prevents cells from becoming unresponsive or indefinitely responsive to the stress, and instead allows cells to respond dynamically. These effects expand the range of stress levels over which the population of cells maintains a heterogeneous response distribution, which is critical to the population viability (see Fig 1B). The use of efficient numerical methods and stochastic simulation at several levels of model complexity allows us to elucidate the mechanisms behind these effects. A central prediction from our modeling is that these effects are rooted in noise arising from low numbers of molecules. We verify this prediction using quantitative fluorescence microscopy by comparing the population response of native B. subtilis with that of synthetic mutants harboring the less-noisy SynExSlow circuit. Taken together, these results constitute a fundamental example of how noise can increase the functionality of a phenotypic response.

Results

Entry of B. subtilis cells into the competent state occurs at high expression levels of the ComK protein. This protein activates a set of downstream genes allowing for the uptake of DNA [15]. ComK is typically expressed at a basal level, and stress in the environment alters the level of expression. In our genetic circuit design, as described below, increasing the stress level is mimicked by inducing comK expression using an increasing amount of a lactose analogue, Isopropyl β-D-1-thiogalactopyranoside (IPTG), in the environment.

In the native competence circuit, ComK activates its own expression, and represses the expression of another protein, ComS. ComS and ComK compete to be rapidly degraded by the MecA protein complex [15] (see Fig 2A, bottom). Therefore, high concentrations of ComS hinder the degradation of ComK, effectively providing positive feedback to ComK by allowing ComK levels to build up. These interactions are summarized in Fig 2A (top). Intrinsic fluctuations in comK mRNA have been shown to be another possible source of variations in ComK proteins [18].

thumbnail
Download:
Fig 2. Architectures and model parameters of the native and SynEx circuits.

The top row summarizes the regulatory interactions, while the bottom row depicts the model details. (A) In the native circuit, ComK is produced with the induction rate αk and activates its own expression with Hill function parameters βk, kk, and h. ComS is expressed at the basal rate αs and is repressed by ComK with Hill function parameters βs, ks, and p. ComK and ComS are degraded at rates λk and λs, respectively, and, additionally, both compete for binding to the degradation enzyme MecA. MecA degrades ComK and ComS with maximal rates δk and δs, respectively, and with Michaelis-Menten constants Γk and Γs, respectively. (B) In the SynEx circuit, ComK is produced with the induction rate αk and activates its own expression with Hill function parameters βk, kk, and h. MecA is expressed at the basal rate αm and is activated by ComK with Hill function parameters βm, km, and p. ComK and MecA are degraded at rates λk and λm, respectively, and MecA enzymatically degrades ComK with rate δ.

https://doi.org/10.1371/journal.pcbi.1004793.g002

In the SynEx circuit, as described in [16], the repression of ComS by ComK is removed by gene knockout. Then, the expression of MecA is placed under the control of ComK. This causes ComK to activate MecA, which in turn represses ComK via active protein degradation (see Fig 2B, bottom). These interactions are summarized in Fig 2B (top). Note that in the native circuit, ComK represses its own activator (ComS), while in the SynEx circuit, ComK activates its own repressor (MecA).

Both the native and SynEx circuits have architectures characteristic of molecular oscillators. Therefore we expect both circuits to allow for a dynamic response of each individual in a population. However, the main difference is that in the native circuit, when ComK levels are high, ComS levels are low, which leads to large amounts of intrinsic noise. In contrast, in the SynEx circuit, when ComK levels are high, MecA levels are also high, corresponding to less intrinsic noise. Previous work showed that this difference in architecture causes the native circuit to display a broad range of competence durations, whereas the SynEx circuit displays a relatively narrow range of competence durations [16]. However, the effects of noise and architecture on the ranges of dynamic response and the ensuing population heterogeneity in these systems remained unknown.

Noise expands the response range

To elucidate the effects of noise in each of the native and SynEx circuits (Fig 2), we develop a stochastic model of each circuit, which includes noise, and then compare each to its deterministic analog, which does not include noise. As described in Materials and Methods, we develop the stochastic models at several levels of complexity to investigate the robustness of our findings to our modeling assumptions, and we solve each model using a combination of efficient numerical solution and stochastic simulation. We first describe the behavior of the deterministic models. As shown in Fig 3A, a standard linear stability analysis of the deterministic model for each circuit reveals three dynamical regimes, depending on the value of the control parameter, the ComK induction rate αk. At low induction, each circuit is excitable, resulting in a transient differentiation event into and out of the competent (high-ComK) state. At intermediate induction, each circuit is oscillatory, periodically entering and exiting the competent state. At high induction, each circuit is mono-stable, staying in the competent state indefinitely; relaxation to the non-competent state does not occur. These three dynamical regimes have been confirmed in experimental studies of the native competence circuit [13].

thumbnail
Download:
Fig 3. Stochastic modeling of B. subtilis competence.

(A) The deterministic model of each circuit (see Eqs 6 and 7) exhibits three dynamic regimes (excitable, oscillatory, and mono-stable), depending on the ComK induction rate αk, which models stress level. (B) The stochastic model (see Eqs 1–5) reveals the ensuing distribution of ComK levels in each of the three dynamic regimes (excitable, oscillatory, and mono-stable). The fraction of the distribution in the responsive state f (determined by the inflection points, see Materials and Methods) is shaded. (C) Whereas the deterministic model exhibits sharp transitions between the dynamic regimes (dashed lines), the stochastic model exhibits a continuous dependence of f on induction rate. We see that for both circuits, stochasticity extends the viable response range (0 < f < 1) beyond the transitions predicted by the deterministic model, in both directions, by the factors given above the arrows (see Materials and Methods). Parameters are as in [16] and are given in S1 Text. In A and B, from left to right, the values of the control parameter are αk = {0.072, 1.15, 36}/hour (native) and αk = {0.036, 1.8, 36}/hour (SynEx). In A, from left to right, the initial conditions are ComK molecules and ComS molecules (native), and ComK molecules and MecA molecules (SynEx); in the excitable regime (left), the initial conditions are chosen to demonstrate the single, transient excitation.

https://doi.org/10.1371/journal.pcbi.1004793.g003

We find that these deterministic dynamics are reflected in the steady-state solutions to the minimal stochastic models. As shown in Fig 3B, the three types of dynamics correspond to three shapes of steady-state probability distributions of ComK levels. Excitable dynamics correspond to a distribution confined to low ComK molecule numbers, oscillatory dynamics correspond to a distribution mixed between low and high molecule numbers, and mono-stable dynamics correspond to a distribution centered at high molecule numbers. As described in Materials and Methods, we calculate the fraction f of the distribution in the high-molecule-number state (see the shaded regions in Fig 3B). Within our model, f represents the fraction of time a single cell spends in the competent state, or equivalently, the fraction of an isogenic population of cells found in the competent state at a given time. Importantly, f is the indicator of population heterogeneity, since unresponsive (f = 0) or fully competent (f = 1) populations are homogeneous, while mixed populations (0 < f < 1) are heterogeneous. We define the range of induction rate αk for which 0 < f < 1 as the viable response range, since unresponsive cells (f = 0) do not benefit from competence, while long-term competence (f = 1) is known to have a detrimental effect on growth rate and cell division [14, 15].

In Fig 3C, we compare the viable response range of the stochastic model with the boundaries between dynamical regimes predicted by the deterministic model. We see that for both the native and the SynEx circuit, the stochastic range extends beyond the deterministic range for both low and high induction rate αk. At low induction rate, the range expands by more in the SynEx circuit than in the native circuit (16 times vs. 8.4 times). In contrast, at high induction rate, the range expands by more in the native circuit than in SynEx circuit (20 times vs. 3.5 times). The latter effect is much stronger, such that the total expansion of the viable range is three times larger in the native circuit than in the SynEx circuit (8.4 × 20 = 168 times vs. 16 × 3.5 = 56 times). Taken together, these observations imply that noise expands the range of stress levels to which cells can respond in a dynamic way, and that this expansion depends on circuit architecture. In the next section, we elucidate the mechanisms behind this expansion.

Noise-induced oscillations underlie the expansion of the response range

Why does noise expand the viable response range at low induction levels? As shown in Fig 4A, the reason is that noise leads to repeated excitations into the competent state, which prevents the system from remaining completely unresponsive. In a completely deterministic excitable system, an excitation is caused by initializing the system away from its stable fixed point, and it occurs only once. However, in a stochastic system, noise can cause repeated perturbations away from the stable state, leading to persistent additional excitations. Indeed, in both circuits, noise at the stable state is high, because the stable state corresponds to one or more species being expressed at very low molecule number. Specifically, in the native circuit, one species is at low molecule number (ComK), and in the SynEx circuit, both species are at low molecule number (ComK and MecA); see Fig 4A. This difference is consistent with the expansion of the viable response range being larger in the SynEx circuit than in the native circuit at low induction levels (Fig 3C). Since the dynamics in both circuits are governed by Poissonian birth-death reactions, low molecule numbers correspond to high intrinsic noise (variance over the squared mean). The corresponding fluctuations are visible in S1A Fig, which shows a zoom-in of the top left panel in Fig 4A. These fluctuations then lead to excitations, driven by the the low-molecule-number species. This is seen in S1B Fig, which overlays the ComK and ComS dynamics for the native circuit in the left column of Fig 4A, and shows that maxima in ComK levels slightly precede minima in ComS levels, indicating that the low-molecule-number species, ComK, is the driver of the excitation. The net result is frequent and persistent noise-induced excitations for both circuits in the low-stress regime. This effect is consistent with the noise-induced excitations seen for these circuits in previous work [15, 16]. Here, however, we have quantified the effect of these excitations on the stochastic distribution, which describes the heterogeneous population response.

thumbnail
Download:
Fig 4. Stochastic oscillations persist outside the deterministic oscillatory regime.

The deterministic oscillatory regime is defined by for the induction rate αk. (A) At low induction rate , where the deterministic model predicts excitable dynamics, the stochastic dynamics are oscillatory. The oscillations arise from repeated noise-induced excitations. (B) At high induction rate , where the deterministic model predicts mono-stable dynamics, the stochastic dynamics are also oscillatory. The oscillations here arise because noise prevents damping to the mono-stable state (see the deterministic curves in the right panels). The effect is much stronger for the native circuit (notice that the left panel is 15 times outside the deterministically oscillatory regime) because, unlike in the SynEx circuit, one of the species, ComS, is at low copy number and therefore subject to significant intrinsic noise. The deterministic model is given in Eqs 6 and 7, while the stochastic model is given in Eqs 1–5. In A, the deterministic initial conditions are ComK molecules and ComS molecules (native), and ComK molecules and MecA molecules (SynEx). In B, the deterministic initial conditions are ComK molecules and ComS molecules (native), and ComK molecules and MecA molecules (SynEx). In the excitable regime (A), the initial conditions are chosen to demonstrate the single, transient excitation.

https://doi.org/10.1371/journal.pcbi.1004793.g004

Why does noise expand the viable response range at high induction levels? Here the mechanism is different from at low induction levels. As shown in Fig 4B, the reason is that noise prevents the damping of oscillations, which keeps the system from relaxing to the competent state. In the deterministic system, the mono-stable state is defined by a stability matrix whose eigenvalues are complex with negative real parts (S2 Fig and S1 Text section 2). This means that the solution relaxes to the mono-stable state in an oscillatory way, i.e. the oscillations are damped. The damped oscillations are clearly visible for the native circuit in S1C and S1D Fig, and for the SynEx circuit in the right panel of Fig 4B (black lines). Intrinsic noise thwarts this relaxation, continually perturbing the system away from the stable point, and preserving a finite oscillation amplitude (colored lines in Fig 4B). Similar effects have been observed in ecological and epidemic models, where they are attributed to the ability of white noise to repeatedly excite a system at its resonant frequency [19]. Here we see the effect at the molecular level in bacteria, and we find that it occurs sufficiently strongly that it supports and significantly extends a heterogeneous population response. Note that in both circuits the low- and high-stress regimes differ in both the number and the stability characteristics of their fixed points (S2 Fig and S1 Text section 2), further indicating that the noise-induced features arise via different mechanisms in the two regimes.

At high induction levels, the expansion of the viable response range is more pronounced in the native circuit than in the SynEx circuit. This effect was demonstrated at the population level in Fig 3C, and it is sufficient to make the expansion of the entire viable response range, at both low and high induction levels, three times larger in the native circuit than in the SynEx circuit. The effect is also demonstrated by the dynamics in Fig 4B: the noise-induced prevention of damping is clearly evident for the native circuit, even at the αk value shown, which is 15 times value predicted deterministically. In fact, at an αk value this large, the deterministic dynamics are no longer damped oscillations; instead they are monostable and non-oscillatory (S1 Text section 2), and yet the stochastic pseudo-oscillations persist. The reason that the effect of noise is so pronounced in the native circuit is that the mono-stable fixed point corresponds to ComS being expressed at very low molecule numbers, where the intrinsic noise is high (lower left panel). In contrast, in the SynEx circuit, the mono-stable state corresponds to both species being expressed at high molecule numbers, so the intrinsic noise is lower than in the native circuit. This difference, which stems ultimately from the difference in the architecture of the two circuits (Fig 2), was found in previous work [16] to be responsible for the increased variability in the competence durations of the native circuit compared to the SynEx circuit. Here we demonstrate that, in the context of the current model, this difference between circuits additionally leads to an increase in the expansion of the viable response range.

We have tested that the effects discussed above are robust, in that they persist when we relax the three simplifying assumptions of our minimal stochastic model (see Materials and Methods). We relax two of the assumptions by considering a non-adiabatic stochastic model in which the fast dynamics of mRNA production and enzymatic degradation that have been named as the possible source of large comK mRNA variations and linked to competence entry [18] are included explicitly, and by setting the mean molecule numbers in the tens of thousands as opposed to tens (see S1 Text). We find that all noise-induced effects persist, namely (i) repeated excitations, (ii) the prevention of damping, and (iii) the enhancement of effect ii in the native circuit over the SynEx circuit (see S6 and S7 Figs). As shown in S6 and S7 Figs, we also verify quantitatively that effects i and ii produce sufficiently oscillatory dynamics that the power spectrum is peaked, as opposed to the non-peaked power spectrum observed for purely excitable or mono-stable dynamics. Interestingly, when we raise the molecule number, but retain the adiabatic assumption, we find that the effects of noise diminish, and the stochastic model behaves like the deterministic model (see S8 and S9 Figs). This confirms that the effects we observe are rooted in the intrinsic noise arising from low molecule numbers, as expected. Importantly, however, it demonstrates that when coupled with explicit mRNA and competitive degradation dynamics, these intrinsic effects dominate the response up to a much higher molecule number regime.

Finally, we relax the third assumption by considering a three-species model for the SynExSlow circuit, in which the dynamics of ComS are accounted for explicitly (see S10 Fig). Recall that SynExSlow is a variant of SynEx that adds an additional PcomG − comS in addition to the PrpsD − comS such that the duration of competence is similar to the native strain [16]. We find that effect ii persists, while effect i does not, indicating that the expansion of the viable response regime at high induction levels is more robust than at low induction levels. Since the high-induction regime is also where the contrast between the native and SynEx circuits is greater, we focus on the high-induction regime in the next section, where we compare our model predictions with experiments.

Fluorescence microscopy confirms the predictions of the model

To test our model predictions, we use quantitative fluorescence microscopy to measure the ComK activity levels in populations of B. subtilis cells harboring either the native or the SynExSlow circuits, as described in Materials and Methods (see Fig 5A). To report the competence state, we use a PcomG − cfp reporter that was previously shown to be highly correlated with the activity of PcomK [15]. In particular, we use a previously reported strain, Hyper − αK, to measure competence in the native competence circuit, and we created an analogous SynExKSlow strain to monitor competence in the synthetic SynExSlow competence circuit. ComK expression is induced by increasing the concentration of IPTG, which corresponds to the model parameter αk. As seen in Fig 5B, in both the native and the SynEx strain, as the IPTG concentration increases, the fluorescence distribution across the population changes shape: first it is centered at low values, then it is split between low and high values, and finally it is centered at high values. This change is qualitatively reminiscent of the change seen in the stochastic model in Fig 3B. Moreover, Fig 5C also shows that the transition to a distribution centered at high values occurs at a higher IPTG concentration in cells with the native circuit than in cells with the SynExSlow circuit. This feature is also qualitatively consistent with the finding in Fig 3C that the viable response range is expanded to a greater extent in the native circuit than in the SynEx circuit at high induction levels.

thumbnail
Download:
Fig 5. Quantitative fluorescence microscopy confirms model predictions.

(A) Microscopy image of B. subtilis cells with 1.5 μM IPTG. ComK activity, measured by CFP fluorescence, is shown in red. (B) Fluorescence distributions over the imaged populations for both circuits as a function of IPTG levels (see additional distributions in S11 Fig). Note that, as in the model, the shift from a non-responsive state (low fluorescence) to a responsive state (high fluorescence) is clearly evident in the distributions. (C) Comparison of the data with the stochastic model (Eqs 1–5) in the high-induction regime. For both circuits, αk is normalized by the value of the deterministic transition from the oscillatory to the mono-stable regime (dashed line); these normalized αk values are shown for each distribution in B in the upper-right corners. Agreement between the model and data confirms both model predictions: that noise extends the viable response regime to higher stress levels than predicted deterministically, and that the effect is more pronounced in the native circuit than in the SynEx circuit.

https://doi.org/10.1371/journal.pcbi.1004793.g005

To investigate whether our experimental observations agree quantitatively with our theoretical predictions, as well as qualitatively, we fit the fluorescence distributions to the stochastic model, as described in Materials and Methods. Fig 5C shows that both of our central predictions in the high-induction regime are quantitatively confirmed by the data, namely (i) that noise extends the transition to a permanently competent state beyond the deterministically predicted induction level, and (ii) that it does so to a larger extent in the native circuit than in the SynEx circuit. Fig 5 therefore provides strong experimental support for the the notion that intrinsic noise expands the viable response range by delaying, as a function of induction level, the relaxation of cells to the competent state.

Discussion

Phenotypic heterogeneity, in which different individuals express particular genes at different levels, is an important survival strategy in uncertain environments. Here we studied dynamically maintained phenotypic heterogeneity in the competence response of B. subtilis, and how it is influenced by intrinsic fluctuations in molecule numbers. By combining theoretical modeling, stochastic simulations, and quantitative microscopy, we showed that intrinsic noise facilitates heterogeneity by expanding the range of stress levels over which heterogeneity is maintained (the viable response range). The effect manifests itself at both low and high stress levels, and the influence of noise is dramatic: in the native competence circuit, noise increases the maximal stress level at which a heterogeneous population response occurs, by 20-fold.

Our work advances previous work investigating the effects of circuit architecture on dynamic response. It was previously known that the native competence circuit exhibited higher variability in its competence duration times than a synthetic analog with different architecture (SynEx). This variability was attributed to intrinsic molecule number fluctuations and was thought to provide a fitness advantage, similar to variability in the times to commit to cell states [20]. Yet the advantages of the native design over the synthetic design were not immediately clear. Here, we showed that while the SynEx circuit is more predictable in terms of competence duration times, its dynamic response is limited to a much smaller range of stress levels, which limits its functionality.

In both circuits, the viable response range is expanded at both low and high stress levels. The mechanisms in these two cases are different. At low stress levels, intrinsic noise causes repeated, period excitations, effectively sustaining oscillations into the excitable regime. At high stress levels, noise prevents the damping of oscillations, effectively delaying, as a function of stress level, the static and indefinite entry into the competent state. Both mechanisms rely on intrinsic fluctuations and, importantly, persist even at high molecule numbers in a non-adiabatic system. The limited response range observed in the deterministic solution is recovered only in the strict limit of fast switching and high molecule numbers, suggesting that the effects of noise that we observe here are generic.

Quantitative microscopy measurements confirmed our theoretical predictions: all cells exhibited competence at high induction levels, no cells exhibited competence at low induction levels, and a bimodal population response was observed in the intermediate regime. Important differences between circuit architectures were also confirmed experimentally, namely that the SynEx circuit begins to oscillate at lower stress levels (IPTG levels) than the native circuit, and that the native circuit can withstand roughly 5-fold higher stress levels than the SynEx circuit before indefinitely entering the competent state. An independent fit of the model predictions to the experimental data showed very good agreement.

Traditionally, noise in gene expression has often been seen as a nuisance that needs to be controlled, especially in stable environments or when the reproducibility of downstream gene expression is crucial. Thus, much work has concentrated on how to ensure the reliability of gene expression and cell signaling in the presence of intrinsic noise [21–29]. However, as has been shown experimentally and theoretically in the context of antibiotic resistance [3, 30, 31], noise-induced population heterogeneity can be advantageous for adaptation to new conditions [32–35]. Functional applications of noise have also been identified in a number of settings [36] ranging from differentiation decisions to sporulate [37], apoptose [38], or allow DNA uptake, such as discussed in this paper.

In B. subtilis, a heterogeneous competence response is thought to be optimal since permanent competence curbs cell growth. The effect of phenotypic heterogeneity on the growth rate of populations has also been studied theoretically [2, 39–44], showing that while in optimal conditions fluctuations decrease the overall growth rate, in less favorable environments, diversity of gene expression increases the population fitness [44]. The effect of selection on such populations was also considered [43], and shown to influence the stability of the phenotypic states [6].

We have taken advantage of the finding that comK is the master regulator of competence in B. subtilis in that the regulation of competence primarily occurs by modifying the expression of comK[45]. However, it is not clear whether significant physiological stress would cause ComK production in the same controlled manner as we induce ComK production via IPTG. Indeed, it is likely that the regulatory input into competence via comK has been evolutionarily tuned to stay within the dynamic range of the native competence machinery. Alternatively, it may be that the competence circuit has been tuned to the accept the available range of comK regulatory input discernable from the organism’s environment. Perhaps then the topology that allows for the most dynamic response range given the available signal transduction mechanisms may be subject to natural selection in order to improve fitness.

An alternative differentiation program to competence, which is evident in the experiments (Fig 5A), is sporulation under stress, in which cells create a durable endospore. This differentiation program has been shown to occur independently of competence without crosstalk [4, 5]. This is remarkably illustrated by the existence of dual-activity cells where the mother cell of a spore continuous to undergo competence until the spore forces it to lyse [5]. Since competence has been demonstrated to occur independently of sporulation, further consideration of sporulation is beyond the scope of this work and spores are excluded from the analysis (see also S1 Text, sections 6 and 7).

We have described an example of phenotypic heterogeneity that is maintained by an oscillatory response, and we have demonstrated that intrinsic noise increases the range of stress levels for which oscillations occur. The ability of noise to facilitate oscillations has also been observed in the entrainment of NF-κB in fibroblast cells to oscillating TNF inputs [46]. There, small-molecule-number noise was shown to facilitate both oscillation and entrainment, and phenotypic variability was shown to enlarge the dynamic range of inputs for which entrainment is possible. These results, along with our findings herein, suggest that strategies that exploit the coupling between noise and phenotypic heterogeneity allow for functional population responses over a large variety of conditions.

Materials and Methods

Stochastic model

Our minimal stochastic models of the native and SynEx circuits are based on our previous modeling work [16], but employ the (stochastic) master equation instead of a (deterministic) dynamical system in order to capture the effects of intrinsic noise. The master equation describes the dynamics of the probability distribution over the numbers of the relevant molecular species inside the cell [47]. For both circuits the master equation reads (1) where pnm is the joint probability distribution over molecule numbers n and m (see Fig 2 for a diagram and explanation of all variables and parameters). In the native circuit, n is the number of ComK proteins and m is the number of ComS proteins. In the SynEx circuit, n is the number of ComK proteins and m is the number of MecA proteins. The dynamics are birth-death processes with mutual regulation: the production rates g and q increase the numbers n and m, respectively, while the degradation rates r and s decrease the numbers n and m, respectively, and the regulation is encoded in the functional dependence of the rates on n and m. The regulation functions follow from our previous work [16] and for the native circuit read (2) (3) while for the SynEx circuit they read (4) (5)

The meaning of the parameters is explained in Fig 2. The regulatory functions introduce positive and negative feedbacks (see Fig 2, top). The parameter values used in the model are as in [16] and are given in S1 Text.

The model in Eqs 1–5 makes three simplifying assumptions, all of which we later relax. First, as in [15, 16] we have assumed that mRNA dynamics and the enzymatic degradation process are substantially faster than all other biochemical reactions in the circuits, and are thus adiabatically eliminated (see S1 Text for details). This reduces each model to the two-species form in Eq 1, depicted by the cartoons in Fig 2 (top). Second, the parameters are chosen such that typical protein copy numbers are small (in the tens or hundreds per cell). Lacking information about the absolute protein numbers in the experiments, we make this assumption because we expect any effects of intrinsic noise to be most evident in the low-number regime, although as we later show, the effects we find persist out to protein numbers in the tens of thousands. Third, we model in Eqs 4 and 5 the SynEx circuit as originally constructed [16], instead of the “SynExSlow” circuit that we use in experiments (described later in this section). This reduces the model from three species to two, which is more amenable to analytic and numerical solution. Once again, however, we will see that the most important effects of noise that we elucidate are also present in a model of the SynExSlow circuit.

We solve Eq 1 in steady state in one of two ways. At low copy numbers, we use the spectral method [48, 49], a hybrid analytic-numerical technique that exploits the eigenfunctions of the birth-death process. Derivation of the spectral solution of Eq 1 is given in S1 Text. The spectral method is much more efficient than other numerical techniques [48], but we find here that it becomes numerically unstable at sufficiently high copy numbers. Therefore, at high copy numbers, we use iterative inversion of the matrix acting on pnm on the right-hand side of Eq 1 (see S1 Text). To obtain individual stochastic trajectories of the system described by Eqs 1–5, we use the Gillespie algorithm [50].

Deterministic model

The deterministic analog of Eq 1 is obtained by performing an expansion in the limit of large molecule numbers [47]. To first order one obtains (6) (7) where and are ensemble averages. Eqs 6 and 7 form a coupled dynamical system whose properties we obtain by linear stability analysis. As shown in S2 Fig, both circuits exhibit excitable, oscillatory, and mono-stable regimes, depending on the value of the control parameter αk. The transition from excitable to oscillatory is marked by the annihilation of a stable and an unstable fixed point, leaving only one unstable fixed point. The transition from oscillatory to mono-stable is marked by this unstable fixed point becoming stable. These transitions provide the dashed lines in Fig 3C.

Genetic circuit construction

For the native competence circuit, we used a variant from our previous study [13]. For the synthetic competence circuit, we reconfigured the original “SynExSlow” circuit created in [16], in order to introduce a tunable proxy for stress level. This required replacing the tunable Phyperspank − comS with an internally controlled promoter for the ribosomal gene rpsD, as well as adding in Phyperspank − comK. The result was a strain that is resistant to four antibiotics and has comS expressed from a ribosomal promoter, providing for a basal level of expression (see S1 Text for chromosomal alterations and antibiotic resistance). In both strains, ComK expression is induced by increasing the amount of IPTG in the environment. Since stress signals are usually integrated at the ComK promoter, IPTG therefore acts as a proxy for stress and triggers competence. This allows us to simulate stress directly in a controlled manner, rather than using physiological stresses that may themselves induce external variation in the responses.

Time-lapse microscopy

Cells of Bacillus subtilis were prepared by streaking from glycerol stocks onto LB agar plates containing the appropriate antibiotic for maintenance and incubated at 37°C overnight. Single colonies were then selected from the plates and grown in LB broth for three to four hours at 37°C until an OD of 1.6 to 1.8 was reached. While culturing the cells, argarose pads were made by pouring 6 mL of 0.8% w/v low-melting point agarose in re-suspension medium onto a glass coverslip. A second glass coverslip was then placed on top of the medium, and the medium was left to congeal while the culture was grown. Once the culture was ready, cells were spun down and resuspended in the resuspension medium twice to wash away the LB. To deposit cells, the top glass coverslip was removed, and then 2 μL of cells were dropped on 37°C low melting point agarose pads. The pads were then cut into squares with a 5mm edge, each containing a single drop of cells. After drying for one additional hour, the pads were flipped over and placed on a glass-bottom dish. The dish was then sealed with parafilm. Images of the cells were then obtained at 100X magnification on an Olympus IX81 system using the ImagePro software from MediaCybernetics along with customized macros.

IPTG stock solutions were dissolved in ethanol to a concentration of 100 mM. Working (1000X) stocks were diluted with Milli-Q water to 30 mM, 10 mM, 3 mM, 1.5 mM, 0.75 mM, respectively, by serial dilution and then added to the appropriate resuspension media at a ratio of 1:999 to achieve the final concentrations indicated. The resuspension media at the specified concentrations of IPTG were used in all steps above after the initial growth in LB. See S12 Fig for time-lapse images of strains containing either the native or SynExSlow competence circuits over the span of 24 hours.

Plasmid and strain construction

Template plasmids with homologous recombination arms for the Bacillus subtilis chromosomal loci were modified through restriction enzyme digest and ligation of DNA inserts (see S1 Text for loci). The inserts were created by polymerase chain reactions using primers from Integrated DNA Technologies while using genomic DNA or other plasmids as templates.

The PY79 strain of Bacillus subtilis was modified through homologous recombination using a One-Step Transformation protocol by inducing competence. 50 ng of plasmid DNA was replicated in TOP10 E. coli cells (Invitrogen, Life Sciences, Inc) and purified using a MiniPrep spin column (Sigma-Aldrich). The DNA was then mixed with culture growing in minimal salts for thirty minutes and then subsequently were rescued using 2xYT rich medium. Positive colonies were then selected on LB agar plates containing selective concentrations of antibiotics.

Image analysis

Fluorescence histograms were obtained from microscopy images using a pixel-based analysis. A mask was created on each image to identify the areas that the cells occupy (see S3 Fig). A histogram of fluorescence intensity values was then generated for pixels within that area.

Culture media

Sterlini-Mandelstram Resuspension Medium was used during time-lapse microscopy and followed the protocol as in references [51, 52]. The actual protocol used consists of making two salt solutions: A and B. Solution A consists of 0.089 g of FeCl3 ⋅ 6H2O, 0.830 g of MgCl2 ⋅ 6H2O and 1.979 g MnCl2 ⋅ 4H2O in 100 mL of filtered water. Solution A is filter sterilized (not autoclaved) and stored at 4°C. Solution B consists of 53.5 g NH4Cl, 10.6 g Na2SO4, 6.8 g KH2PO4, and 9.7 g NH4NO3. Solution B is then also filter sterilized and stored at 4°C. Sporulation salts are made by combining adding 1 mL of Solution A and 10 mL of Solution B to filtered water for a total 1 L. This solution is then autoclaved. The final Resuspension media is created by combining 93 mL of sporulation salts, 2 mL of 10% v/v L-glutamate, 1 mL of 0.1M CaCl2, and 4 mL of 1M MgSO4 on the day of the experiment.

One-step transformation media consists of 6.25 g of K2HPO4 ⋅ 3H2O, 1.5 g of KH2PO4, 0.25 g of trisodium citrate, 50 mg of MgSO4 ⋅ 7H2O, 0.5 g of Na2SO4 at pH 7.0, 125 μL of 100 mM FeCl3, 5 μL of 100 mM MnSO4, 1 g of glucose, and 0.5 g of glutamate added into filtered water for a total of 250 mL. The media is filter sterilized using 0.2 micron Millipore filters.

2xYT recovery medium consists of 16.0 g of Tryptone, 10.0 g of Yeast Extract, 5.0 g of NaCl added to filtered water to a total volume of 1L. The media is then filter sterilized using 0.2 micron Millipore filters.

Distribution analysis and comparing experiments with modeling

For the ComK distributions in the model, pn = ∑m pnm, we determine the fraction f of cells in the responsive, high ComK protein concentration state using two independent methods. First, we use a generalized method of separating the distribution’s two modes: since the distribution is often not completely bimodal (see Fig 3B, middle column), we find the average n* = (n1 + n2) / 2 of the two inflection points surrounding the putative local minimum between the two modes, and define . In the case of a bimodal pn, this method indeed well approximates the location of the actual local minimum. Second, we fit pn to a mixture of two Poisson distributions, using the Kullback-Leibler divergence as the cost function. There are three fitting parameters, the two Poisson parameters and the relative weighting between them, and the weighting provides f. We see in S4 Fig that the two methods give similar results for the dependence of f on the control parameter αk, demonstrating that our determination of f is robust to the method used.

The two-Poisson method is smoother but is numerically unstable at low αk. Furthermore, it does not capture the transition to f = 1 at high αk, since a roughly equal mixture of two Poisson distributions with similar means (f ∼ 0.5) will always provide a better fit than a single Poisson distribution (f = 1). Therefore, in Fig 3C, we use the two-Poisson method at intermediate αk, while we use the inflection-point method for low αk and to determine the transition to f = 1 at high αk. At low/intermediate αk we concatenate the curves from the two methods and use smoothing near the concatenation point. Since f = 0 is approached asymptotically, we quantify the expansion of the viable range at low stress by determining the minimum αk value for which f > ϵ = 10−2. The expansion factors at low αk are robust to the values of ϵ and the smoothing parameters.

To compare the model predictions to the experimental data, we analyze the fluorescence distributions (Fig 5B) in the same way as the model distributions. Specifically, we calculate the fraction of the comK population in the responsive state f for each experimental distribution using the two-Poisson method above. Because the mapping between IPTG concentration and the model parameter describing the level of external stress αk is unknown, we infer the most likely value of αk corresponding to each experimental distribution by fitting the theory to the data. First we use the mode of the [IPTG] = 0 distributions to subtract the background fluorescence from the remaining data. To avoid binning, we then fit the cumulative distribution instead of the probability distribution (sample fits are shown in S5 Fig). We use a maximally constrained least-squares fit, where all parameters are fixed as in S1 Text except αk and the unknown parameter X describing the conversion of pixel intensity to molecule number. Given a value of X, we find the values of αk for each distribution that minimize the sum of the squared error S in each case. X is then chosen by minimizing the sum of minimum S values over all distributions. These f and αk values inferred from the data are plotted in Fig 5C. The error bars on αk are obtained by finding the αk values where S reaches 1.25 of its minimum value (sample plots of S vs. αk are shown in S5 Fig).

Supporting Information

S1 Text. Supplementary materials and methods, including model parameter values, stability analysis of the deterministic fixed points, description of relaxing the model assumptions, spectral solution to the master equation, and chromosomal alterations and antibiotic resistance of the strains used.

https://doi.org/10.1371/journal.pcbi.1004793.s001

(PDF)

S1 Fig. Further details on simulations in Fig 4.

(A) Zoom-in of the top left panel in Fig 4A, showing fluctuations at low molecule number. (B) Overlay of the ComK and ComS dynamics in the left column of Fig 4A, demonstrating that ComK and ComS excitations are synchronized, with maxima in ComK slightly preceding minima in ComS. (C, D) As in the left column of Fig 4B, except less far outside the deterministically oscillatory regime (1.5 times the deterministic transition value of αk, instead of 15 times). In this regime, the deterministic dynamics are clearly damped oscillatory (black), while the stochastic dynamics are, as in Fig 4B, pseudo-oscillatory and not damped.

https://doi.org/10.1371/journal.pcbi.1004793.s002

(PDF)

S2 Fig. Linear stability analysis of deterministic system predicts three dynamic regimes, for both the native and SynEx circuits.

Fixed points satisfying the steady state of Eqs 6 and 7 for (A) the native circuit and (B) the SynEx circuit. Each fixed point is stable if the real parts of the eigenvalues of the Jacobian matrix evaluated at that point are negative, and unstable otherwise. The Jacobian matrix is Jij = ∂Fi/∂xj, where , , , and . For both circuits, there are three dynamic regimes. The excitable regime (low αk) is has three fixed points, one of which is stable. The oscillatory regime (intermediate αk) has one unstable fixed point. The mono-stable regime (high αk) has one stable fixed point. In the mono-stable regime, near the oscillatory regime, the eigenvalues are complex, indicating damped oscillations (see also S1 Text section 2).

https://doi.org/10.1371/journal.pcbi.1004793.s003

(PDF)

S3 Fig. Image analysis procedure.

Fluorescence histograms are generated via isolating cells from the background by identifying connected and contiguous areas, applying a binary mask, and binning the resulting pixel intensities. Only pixels that fall within cells, not the background, are included.

https://doi.org/10.1371/journal.pcbi.1004793.s004

(PDF)

S4 Fig. Determination of high-response fraction f is robust to calculation method.

Two independent methods are used to determine f: finding the inflection points, and fitting to a mixture of two Poisson distributions (see Materials and Methods). For both the native and SynEx circuit, we see that the two methods give results that correspond very closely to each other. The two-Poisson method does not capture the transition to f = 1, since a roughly equal mixture of two Poisson distributions with similar means (f ∼ 0.5) will always provide a better fit than a single Poisson distribution (f = 1).

https://doi.org/10.1371/journal.pcbi.1004793.s005

(PDF)

S5 Fig. Sample fits to experimental data and plots of sum-of-squares.

Cumulative probability distributions of fluorescence data are fit to the distributions from the stochastic model (Eqs 1–5) by minimizing the sum of squared errors, for the (A-D) native and (E-H) SynEx circuits, as described in Materials and Methods.

https://doi.org/10.1371/journal.pcbi.1004793.s006

(PDF)

S6 Fig. Effects of noise persist when assumptions are relaxed in the model of the native circuit.

Top two rows show ComK and ComS time series from Gillespie simulations of the relaxed model that includes mRNA and competitive degradation dynamics (see S1 Text, section 3, native circuit). Far from the deterministic boundaries of the control parameter, and (indicated by the dashed vertical lines), the dynamics are excitable, oscillatory, and mono-stable as predicted (columns 1, 3, and 5, respectively). However, near the boundaries, but outside the oscillatory regime, noise causes oscillations to persist, due to either repeated excitations (column 2) or prevention of damping (column 4), confirming the effects seen in the reduced model of the main text. The persistence of oscillations is verified by computing the power spectrum from the Fourier transform of the ComK time series n(t). For periodic signals, the power spectrum is peaked at a non-zero frequency ω (and in some cases its harmonics). Red line is a Gaussian fit to aid the eye.

https://doi.org/10.1371/journal.pcbi.1004793.s007

(PDF)

S7 Fig. Effects of noise persist when assumptions are relaxed in the model of the SynEx circuit.

As in S6 Fig but for the SynEx circuit (see S1 Text, section 3, SynEx circuit). Once again, oscillations persist outside the deterministic boundaries as indicated by the peaked power spectra. Note, however, that oscillations are damped at the value here, whereas in the native circuit they persist beyond this value (see S6 Fig). This confirms the effect seen in the main text that the the prevention of damping is more pronounced in the native circuit than in the SynEx circuit.

https://doi.org/10.1371/journal.pcbi.1004793.s008

(PDF)

S8 Fig. Effects of noise diminish when molecule number is raised in the adiabatic model of the native circuit.

Gillespie simulations of the adiabatically reduced model of the native circuit (Eqs 1–5, as in Fig 4), but for high molecule numbers (Γk = 25000 and Γs = 20). We see that 10% outside the deterministically oscillatory regime, the stochastic dynamics are either non-oscillatory (column 2) or weakly oscillatory (column 4). This is in contrast to the low-molecule-number regime (Fig 4), where oscillations persist in these regions and beyond. We conclude that raising molecule number in the adiabatic model of the native circuit reduces the stochastic behavior to the deterministic behavior.

https://doi.org/10.1371/journal.pcbi.1004793.s009

(PDF)

S9 Fig. Effects of noise diminish when molecule number is raised in the adiabatic model of the SynEx circuit.

As in S8 Fig, but for the adiabatically reduced model of the SynEx circuit (Eqs 1–5) at high molecule numbers (kk = 5000 and km = 2500). Comparing to Fig 4, we similarly conclude that raising molecule number in the adiabatic model of the SynEx circuit reduces the stochastic behavior to the deterministic behavior.

https://doi.org/10.1371/journal.pcbi.1004793.s010

(PDF)

S10 Fig. Effects of noise persist in a model of the SynExSlow circuit.

Top row shows the deterministic ComK time series from the SynExSlow model (see S1 Text, section 3, SynExSlow circuit), while the next three rows show the stochastic ComK, MecA, and ComS time series for the same model. Although the SynExSlow model only exhibits a damped oscillatory regime at these parameters, not a standard oscillatory regime (see S1 Text), we define a heuristic boundary /s below which oscillations are not appreciably damped within the first 24 hours (column 3), and above which they are (column 4). We see that, as in S7 Fig, noise prevents damping at large values of the control parameter, even at high molecule numbers (column 4). However, as in S9 Fig, noise does not induce repeated excitations at small values of the control parameter (column 2). We conclude that the former effect is more robust.

https://doi.org/10.1371/journal.pcbi.1004793.s011

(PDF)

S11 Fig. Supplementary fluorescence distributions.

Together with Fig 5B, these distributions provide the data analyzed for Fig 5C. As in Fig 5B, the normalized αk values are shown for each distribution in the upper-right corners.

https://doi.org/10.1371/journal.pcbi.1004793.s012

(PDF)

S12 Fig. Time-lapse fluorescence images of native and SynExSlow strains.

Composite phase contrast and fluorescence images of the Native and SynExSlow strains are shown with IPTG concentrations of 0, 0.75, 1.5, and 3 μM at 0, 6, 12, and 24 hours. Fluorescence represents ComK activity as reported by PcomG − cfp.

https://doi.org/10.1371/journal.pcbi.1004793.s013

(PDF)

Acknowledgments

We thank Tolga Cagatay for assistance in preparing Bacillus strains.

Author Contributions

Conceived and designed the experiments: AM MK CHW GMS AMW. Performed the experiments: AM MK LH JL. Analyzed the data: AM MK LH GMS AMW. Contributed reagents/materials/analysis tools: AM MK CHW GMS AMW. Wrote the paper: AM MK LH JL GMS AMW.

References

  1. 1. Fraser D, Kaern M. A chance at survival: gene expression noise and phenotypic diversification strategies. Molecular microbiology. 2009;71(6):1333–1340. pmid:19220745
  2. 2. Kussell E, Leibler S. Phenotypic diversity, population growth, and information in fluctuating environments. Science. 2005;309(5743):2075–2078. pmid:16123265
  3. 3. Balaban NQ, Merrin J, Chait R, Kowalik L, Leibler S. Bacterial persistence as a phenotypic switch. Science. 2004;305(5690):1622–1625. pmid:15308767
  4. 4. Kuchina A, Espinar L, Garcia-Ojalvo J, Suel G. Reversible and noisy progression towards a commitment point enables adaptable and reliable cellular decision-making. PLoS Comput Biol. 2011;7(11):e1002273. pmid:22102806
  5. 5. Kuchina A, Espinar L, Çağatay T, Balbin AO, Zhang F, Alvarado A, et al. Temporal competition between differentiation programs determines cell fate choice. Molecular systems biology. 2011;7(1):557. pmid:22146301
  6. 6. Mora T, Walczak AM. Effect of phenotypic selection on stochastic gene expression. The Journal of Physical chemistry B. 2013;117(42):13194–13205. pmid:23795617
  7. 7. Novick A, Weiner M. Enzyme induction as an all-or-none phenomenon. Proceedings of the National Academy of Sciences of the United States of America. 1957;43(7):553. pmid:16590055
  8. 8. Choi PJ, Cai L, Frieda K, Xie XS. A stochastic single-molecule event triggers phenotype switching of a bacterial cell. Science. 2008;322(5900):442–446. pmid:18927393
  9. 9. Savageau MA. Comparison of classical and autogenous systems of regulation in inducible operons. Nature. 1974;252:546–549. pmid:4431516
  10. 10. Kittisopikul M, Süel GM. Biological role of noise encoded in a genetic network motif. Proceedings of the National Academy of Sciences. 2010;107(30):13300–13305.
  11. 11. Koh RS, Dunlop MJ. Modeling suggests that gene circuit architecture controls phenotypic variability in a bacterial persistence network. BMC systems biology. 2012;6(1):47. pmid:22607777
  12. 12. Michel D. Kinetic approaches to lactose operon induction and bimodality. Journal of theoretical biology. 2013;325:62–75. pmid:23454080
  13. 13. Süel GM, Kulkarni RP, Dworkin J, Garcia-Ojalvo J, Elowitz MB. Tunability and noise dependence in differentiation dynamics. Science. 2007;315(5819):1716–1719. pmid:17379809
  14. 14. Haijema BJ, Hahn J, Haynes J, Dubnau D. A ComGA-dependent checkpoint limits growth during the escape from competence. Molecular microbiology. 2001;40(1):52–64. pmid:11298275
  15. 15. Süel GM, Garcia-Ojalvo J, Liberman LM, Elowitz MB. An excitable gene regulatory circuit induces transient cellular differentiation. Nature. 2006;440(7083):545–550. pmid:16554821
  16. 16. Çağatay T, Turcotte M, Elowitz MB, Garcia-Ojalvo J, Süel GM. Architecture-dependent noise discriminates functionally analogous differentiation circuits. Cell. 2009;139(3):512–522. pmid:19853288
  17. 17. Espinar L, Dies M, Çağatay T, Süel GM, Garcia-Ojalvo J. Circuit-level input integration in bacterial gene regulation. Proceedings of the National Academy of Sciences. 2013;110(17):7091–7096.
  18. 18. Maamar H, Raj A, Dubnau D. Noise in gene expression determines cell fate in Bacillus subtilis. Science. 2007;317(5837):526–529. pmid:17569828
  19. 19. Black AJ, McKane AJ. Stochastic formulation of ecological models and their applications. Trends in ecology & evolution. 2012;27(6):337–345.
  20. 20. Nachman I, Regev A, Ramanathan S. Dissecting timing variability in yeast meiosis. Cell. 2007;131(3):544–556. pmid:17981121
  21. 21. Mugler A, Tostevin F, ten Wolde PR. Spatial partitioning improves the reliability of biochemical signaling. Proceedings of the National Academy of Sciences. 2013;110(15):5927–5932.
  22. 22. Govern CC, ten Wolde PR. Optimal resource allocation in cellular sensing systems. Proceedings of the National Academy of Sciences. 2014;111(49):17486–17491.
  23. 23. Govern CC, ten Wolde PR. Energy dissipation and noise correlations in biochemical sensing. Physical review letters. 2014;113(25):258102. pmid:25554909
  24. 24. Tostevin F, Wolde PR, Howard M. Fundamental Limits to Position Determination by Concentration Gradients. PLoS Computational Biology. 2007;3(4). pmid:17465676
  25. 25. Saunders TE, Howard M. Morphogen profiles can be optimized to buffer against noise. Physical Review E. 2009;80(4):41902.
  26. 26. Tkačik G, Walczak AM, Bialek W. Optimizing information flow in small genetic networks. Physical Review E. 2009;80(3):031920.
  27. 27. Walczak AM, Tkačik G, Bialek W. Optimizing information flow in small genetic networks. II. Feed-forward interactions. Physical Review E. 2010;81(4):041905.
  28. 28. Tkačik G, Walczak AM, Bialek W. Optimizing information flow in small genetic networks. III. A self-interacting gene. Physical Review E. 2012;85(4):041903.
  29. 29. Tkačik G, Walczak AM. Information transmission in genetic regulatory networks: a review. Journal of Physics: Condensed Matter. 2011;23(15):153102. pmid:21460423
  30. 30. Rotem E, Loinger A, Ronin I, Levin-Reisman I, Gabay C, Shoresh N, et al. Regulation of phenotypic variability by a threshold-based mechanism underlies bacterial persistence. Proceedings of the National Academy of Sciences. 2010;107(28):12541–12546.
  31. 31. Gefen O, Balaban NQ. The importance of being persistent: heterogeneity of bacterial populations under antibiotic stress. FEMS microbiology reviews. 2009;33(4):704–717. pmid:19207742
  32. 32. Sato K, Ito Y, Yomo T, Kaneko K. On the relation between fluctuation and response in biological systems. Proceedings of the National Academy of Sciences. 2003;100(24):14086–14090.
  33. 33. Ito Y, Toyota H, Kaneko K, Yomo T. How selection affects phenotypic fluctuation. Mol Syst Biol. 2009;5:264. pmid:19401676
  34. 34. Kashiwagi A, Urabe I, Kaneko K, Yomo T. Adaptive response of a gene network to environmental changes by fitness-induced attractor selection. PloS one. 2006;1(1):e49. pmid:17183678
  35. 35. Shimizu Y, Tsuru S, Ito Y, Ying B, Yomo T, Kudla G. Stochastic Switching Induced Adaptation in a Starved Escherichia coli Population. PLoS ONE. 2011;6(9).
  36. 36. Eldar A, Elowitz MB. Functional roles for noise in genetic circuits. Nature. 2010;467(7312):167–173. pmid:20829787
  37. 37. Veening JW, Stewart EJ, Berngruber TW, Taddei F, Kuipers OP, Hamoen LW. Bet-hedging and epigenetic inheritance in bacterial cell development. Proceedings of the National Academy of Sciences. 2008;105(11):4393–4398.
  38. 38. Spencer SL, Gaudet S, Albeck JG, Burke JM, Sorger PK. Non-genetic origins of cell-to-cell variability in TRAIL-induced apoptosis. Nature. 2009;459(7245):428–432. pmid:19363473
  39. 39. Thattai M, Van Oudenaarden A. Stochastic gene expression in fluctuating environments. Genetics. 2004;167(1):523–530. pmid:15166174
  40. 40. Kussell E, Kishony R, Balaban NQ, Leibler S. Bacterial persistence a model of survival in changing environments. Genetics. 2005;169(4):1807–1814. pmid:15687275
  41. 41. Leibler S, Kussell E. Individual histories and selection in heterogeneous populations. Proceedings of the National Academy of Sciences. 2010;107(29):13183–13188.
  42. 42. Rivoire O, Leibler S. The value of information for populations in varying environments. Journal of Statistical Physics. 2011;142(6):1124–1166.
  43. 43. Sato K, Kaneko K. On the distribution of state values of reproducing cells. Physical biology. 2006;3(1):74. pmid:16582472
  44. 44. Tănase-Nicola S, Ten Wolde PR. Regulatory control and the costs and benefits of biochemical noise. PLoS Comput Biol. 2008;4(8):e1000125–e1000125. pmid:18716677
  45. 45. Hahn J, L K, D D. The regulation of competence transcription factor synthesis constitues a critical control in the regulation of competence in Bacillus subtilis. J Bacteriology. 1994;176(18):5753–5761.
  46. 46. Kellogg RA, Tay S. Noise facilitates transcriptional control under dynamic inputs. Cell. 2015;160(3):381–392. pmid:25635454
  47. 47. Van Kampen NG. Stochastic processes in physics and chemistry. vol. 1. Elsevier; 1992.
    • 48. Walczak AM, Mugler A, Wiggins CH. A stochastic spectral analysis of transcriptional regulatory cascades. Proceedings of the National Academy of Sciences. 2009;106(16):6529–6534.
    • 49. Mugler A, Walczak AM, Wiggins CH. Spectral solutions to stochastic models of gene expression with bursts and regulation. Physical Review E. 2009;80(4):041921.
    • 50. Gillespie DT. Exact stochastic simulation of coupled chemical reactions. The Journal of Physical Chemistry. 1977;81(25):2340–2361.
    • 51. Sterlini JM, Mandelstam J. Commitment to sporulation in Bacillus subtilis and its relationship to development of actinomycin resistance. Biochem J. 1969 Jun;113(1):29–37. pmid:4185146
    • 52. Harwood CR, Cutting SM. Molecular Biological Methods for Bacillus. West Sussex, England: John Wiley & Sons Ltd; 1990.