Noise Expands the Response Range of the Bacillus subtilis Competence Circuit

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.


II. 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 condi-tions 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][2][3][4][5][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][3][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][10][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. Un-arXiv:1508.07571v1 der 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.
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]. 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. 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 SynEx circuit. Taken together, these results constitute a fundamental In the native circuit, ComK is produced with the induction rate α k and activates its own expression with Hill function parameters β k , k k , 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 , k k , 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 δ.
example of how noise can increase the functionality of a phenotypic response.

III. 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).
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.

A. 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 Sec. V, 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. These three dynamical regimes have been confirmed in experimental studies of the native competence circuit [13].
We find that these deterministic dynamics are reflected in the stationary solutions to the minimal stochastic models. As shown in Fig. 3B, the three types of dynamics correspond to three shapes of stationary 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 Sec. V, we calculate the fraction f of the distribution in the highmolecule-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 for the induction rate α k . (A) At low induction rate α k < α (1) k , 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 α k > α (2) k , 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.
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 . Furthermore, the extension at high induction rate is significantly more pronounced for the native circuit (roughly 20 times the deterministic value) than for the SynEx circuit (roughly 3 times the deterministic value). These observations imply that noise expands the range of stress levels to which cells can respond in a dynamic way. In the next section, we elucidate the mechanisms behind this expansion.

B. 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 (ComK for the native circuit, ComK and MecA for the SynEx circuit; see Fig. 4A). Since the dynamics are governed by Poissonian birth-death reactions, low molecule numbers correspond to high intrinsic noise (variance over the squared mean), leading to frequent and persistent excitations. 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.
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 (Fig. A1). This means that the solution relaxes to the mono-stable state in an oscillatory way, i.e. the oscillations are damped (see the black lines in the right panel of Fig. 4B, for example). Intrinsic noise thwarts this relaxation, continually perturbing the system away from the stable point, and preserving a finite oscillation amplitude (see the colored lines). 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 [17]. 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.
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. It 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. The reason that the effect 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 begin expressed at higher molecule numbers, so the intrinsic noise is lower. 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 the architecture of the native circuit additionally leads to an increase in the expansion of its viable response range, which has a clear benefit for fitness.
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 Sec. V). We relax two of the assumptions by considering a non-adiabatic stochastic model in which the fast dynamics of mRNA production and enzymatic degradation are included explicitly, and by setting the mean molecule numbers in the tens of thousands as opposed to tens (see Appendix A). 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 Fig. A5 and Fig.  A6). As shown in Fig. A5 and Fig. A6, we also verified 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 Fig. A7 and Fig. A8). 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 Fig. A9). 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 this is also the more pronounced effect, we focus on the high-induction regime in the next section, where we compare our model predictions with experiments.

C. Fluorescence microscopy confirms the predictions of the model
To test our model predictions, we use quantitative fluorescence microscopy to measure the ComK expression levels in populations of B. subtilis cells harboring either the native or the SynExSlow circuits, as described in Sec. V (see Fig. 5A). 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 theoretical prediction shown in Fig. 3C.
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 Sec. V. 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. . 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.

IV. 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 8 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 [18]. 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 [19][20][21][22][23][24][25][26][27]. However, as has been shown experimentally and theoretically in the context of antibiotic resistance [3,28,29], noise-induced population heterogeneity can be advantageous for adaptation to new conditions [30][31][32][33]. Functional applications of noise have also been identified in a number of settings [34] ranging from dif-ferentiation decisions to sporulate [35], apoptose [36], 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,[37][38][39][40][41][42], showing that while in optimal conditions fluctuations decrease the overall growth rate, in less favorable environments, diversity of gene expression increases the population fitness [42]. The effect of selection on such populations was also considered [41], and shown to influence the stability of the phenotypic states [6].
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 [43]. There, small-moleculenumber 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.

A. 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 [44]. For both circuits the master equation reads dp nm dt = g n−1 p n−1,m + r n+1,m (n + 1)p n+1,m +q n p n,m−1 + s n,m+1 (m + 1)p n,m+1 −(g n + r nm n + q n + s nm m)p nm .
where p nm 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 9 functional dependence of the rates on n and m. The regulation functions follow from our previous work [16] and for the native circuit read while for the SynEx circuit they read 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 Appendix A. The model in Eqns. 1-9 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 Appendix A for details). This reduces each model to the two-species form in Eqn. 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 Eqns. 6-9 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 Eqn. 1 in steady state in one of two ways. At low copy numbers, we use the spectral method [45,46], a hybrid analytic-numerical technique that exploits the eigenfunctions of the birth-death process. Derivation of the spectral solution of Eqn. 1 is given in Appendix A. The spectral method is much more efficient than other numerical techniques [45], 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 p nm on the righthand side of Eqn. 1 (see Appendix A). To obtain individual stochastic trajectories of the system described by Eqns. 1-9, we use the Gillespie algorithm [47].

B. Deterministic model
The deterministic analog of Eqn. 1 is obtained by performing an expansion in the limit of large molecule numbers [44]. To first order one obtains wheren andm are ensemble averages. Eqns. 10 and 11 form a coupled dynamical system whose properties we obtain by linear stability analysis. As shown in Fig.  A1, both circuits exhibit excitable, oscillatory, and monostable 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.

C. 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 P hyperspank − comS with an internally controlled promoter for the ribosomal gene rpsD, as well as adding in P hyperspank − 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 Appendix A 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.

D. 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 resuspension 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 way 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 media at a ratio of 1:999 to achieve the final concentrations indicated.

E. 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 Appendix A 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.

F. 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 Fig. A2). A histogram of fluorescence intensity values was then generated for pixels within that area.

G. Culture media
Sterlini-Mandelstram Resuspension Medium was used during time-lapse microscopy and followed the protocol as in references [48,49]. The actual protocol used consists of making two salt solutions: A and B. Solution A consists of 0.089 g of FeCl 3 ·6H 2 O, 0.830 g of MgCl 2 ·6H 2 O and 1.979 g MnCl 2 ·4H 2 O 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 NH 4 Cl, 10.6 g Na 2 SO 4 , 6.8 g KH 2 PO 4 , and 9.7 g NH 4 NO 3 . 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 CaCl 2 , and 4 mL of 1M MgSO 4 on the day of the experiment.
One-step transformation media consists of 6.25 g of K 2 HPO 4 · 3H 2 O, 1.5 g of KH 2 PO 4 , 0.25 g of trisodium citrate, 50 mg of MgSO 4 · 7H 2 O, 0.5 g of Na 2 SO 4 at pH 7.0, 125 µL of 100 mM FeCl 3 , 5 µL of 100 mM MnSO 4 , 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.

H. Distribution analysis and comparing experiments with modeling
For the ComK distributions in the model, p n = m p nm , we determine the fraction f of 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 * = (n 1 + n 2 )/2 of the two inflection points surrounding the putative local minimum between the two modes, and define f = ∞ n=n * p n . In the case of a bimodal p n , this method indeed well approximates the location of the actual local minimum. Second, we fit p n 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 Fig. A3 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.
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 Fig. A4). We use a maximally constrained least-squares fit, where all parameters are fixed as in Appendix A 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 Fig. A4). dimensionless form, which is sufficient to establish their deterministic dynamical properties. To account for intrinsic noise and to compare with experiments, several physical quantities are then specified to establish molecule numbers and timescales. The conversion to the remaining model parameters, and their resulting values, are then given below.

Native circuit:
Dimensionless parameters Physical quantities Molecule numbers: Γ k = 100 Γ s = 1 Timescales: δ k = 0.001/s δ s = 0.001/s Remaining parameters Note that in [13,15,16], the higher values Γ k = 25000 and Γ s = 20 (native), and k k = 5000 and k m = 2500 (SynEx) are used to establish molecule number. We use the lower values here to elucidate the effects of intrinsic noise. This assumption is relaxed in the next section, where we return to the high-number regime. In the main text, we discuss how the effects of noise are robust to this choice.

Relaxing the model assumptions
To relax the first two simplifying assumptions made in the main text, we consider a stochastic model that has high molecular copy numbers, and that accounts explicitly for the mRNA dynamics and the dynamics of competitive degradation. The model follows from our earlier work [16], and consists of a set of coupled chemical reactions for each circuit. Due to the complexity of the model, we do not solve the master equation explicitly, but rather we simulate the dynamics using the Gillespie algorithm [47].
For the native circuit, the reactions and rates are: Additional parameters k 2 = 0.19/s k 5 = 0.0015/s k 11 = 2 × 10 −6 /s k 13 = 4.5 × 10 −6 /s k k = 5000 k s = 833 h = 2 p = 5 Ω = 1.66 µm 3 M T = 500 Here P gene denotes the promoter of the corresponding gene (constitutive or regulated), and MecA K and MecA S represent the complex of MecA bound to ComK and ComS, respectively. As in the main text, n is the number of ComK molecules per cell. Ω represents the cell volume, and M T gives the total number of MecA molecules. Upon adiabatically eliminating the faster mRNA dynamics and the MecA dynamics, one obtains the model in the previous section, except with the larger molecule numbers Γ k = 25000 and Γ s = 20 [13]. The control parameter here, k 1 , is related to the control parameter in the reduced model, α k , via k 1 = k 7 α k /k 3 [13].
For the SynEx circuit, the reactions and rates are: Once again, P gene denotes the promoter of the corresponding gene, n is the number of ComK molecules per cell, and Ω represents the cell volume. Upon adiabatically eliminating the faster mRNA dynamics, one obtains the model in the previous section, except with the larger molecule numbers k k = 5000 and k m = 2500 [16]. The control parameter here, k 1 , is related to the control parameter in the reduced model, α k , via k 1 = k 8 α k /k 5 .
To relax the third simplifying assumption made in the main text, we consider a stochastic model of the SynExSlow circuit. The model follows from our earlier work [16], and is similar to the stochastic model in the main text, except that there are three species with dynamic molecule numbers: ComK (n), MecA (m), and ComS ( ). The production rate functions are, respectively, q n = α m + β m n p k p m + n p , and the degradation rate functions are, respectively, The parameter values are [16]: Parameters α k (varied) β s = 0.5/s α m = 0.075/s Γ k = 25000 α s = 0.5/s Γ s = 20 δ k = δ s = 2 × 10 −6 /s k k = 5000 λ k = λ m = λ s = 10 −4 /s k m = 2500 β k = 7.5/s k s = 500 β m = 2.5/s h = p = 2 Note that these values correspond to the high-molecule-number regime of the native and SynEx models of the main text, and thus this model of the SynExSlow circuit also relaxes the first simplifying assumption of low molecule number.
Finally, we note that at these parameter values, the deterministic analog of this model predicts that as a function of the control parameter α k , the excitable regime transitions directly into a damped oscillatory (mono-stable) regime, where two of the three eigenvalues of the Jacobian matrix are complex, and all have negative real part. That is, there is no standard oscillatory regime. Therefore, we define a heuristic boundary α (2) k = 0.15/s after which the damping is clearly evident within the first 24 hours of the deterministic dynamics (see Fig. S7). The fact that the stochastic dynamics exhibit sustained oscillations in the absence of a deterministically oscillatory regime, even beyond this heuristic boundary (Fig. S7), indicates that the ability of noise to induce oscillations in the SynExSlow circuit is especially strong.

Spectral solution to the master equation
We write the master equation (Eqn. 1) as dp nm dt = − L n [g n , r nm ] +L m [q n , s nm ] p nm .
Here −L is the linear birth-death operator, whose action on the probability distribution p is described by −L n [g n , r n ]p n = g n−1 p n−1 + r n+1 (n + 1)p n+1 − (g n + r n n)p n , where in general on an operator (here,L) we use a subscript to denote the sector (n or m) on which it acts. Explicitly, then, the master equation reads dp nm dt = g n−1 p n−1,m + r n+1,m (n + 1)p n+1,m − (g n + r nm n)p nm +q n p n,m−1 + s n,m+1 (m + 1)p n,m+1 − (q n + s nm m)p nm .
We will now derive the spectral decomposition of the master equation by introducing the generating function. For intuition, we will first introduce the generating function in the context of the one-dimensional system described by Eqn. A4, then extend our results to the full master equation.
The generating function is an expansion in a complete set of states, indexed by molecule number, for which the probabilities provide the expansion coefficients. Denoting the states abstractly as |n , the generating function is defined With Eqn. A34, we have turned the original master equation (Eqn. A17) into a linear algebraic equation, involving only matrix and tensor multiplication. In steady state Eqn. A34 can be solved perturbatively by treating the first term on the right-hand side as larger than the others and iterating to convergence. Algorithmically, then, the procedure prescribed by the spectral method is the following. First we compute the deviation matrices and tensors from the given rate functions g n , r nm , q n , and s nm , and our chosen 'gauges'ḡ,r,q, ands. Then we solve for the steady state of Eqn. A34 by iteration to obtain the expansion coefficients G jk . Finally, we compute the probability distribution by taking the inverse transform: The probability distribution provides the complete stochastic description of the steady-state process. Note that both the deviation matrices/tensors and the probability distribution are computed via multiplication against the eigenmodes n|j and m|k or their conjugates j|n and k|m ; these are efficiently precomputed via recursive update rules. For the SynEx circuit, the dynamics of the spectral expansion coefficients are simplified because we have r nm → r m and s nm →s:Ġ where We see that the sparser coupling of the degradation terms (compared to the native circuit) has left us with a linear algebraic equation involving only matrices, not tensors. The solution follows the steps outlined above for the native circuit.
The spectral method is more efficient than other solution techniques by many orders of magnitude [45]. However, we find that here it becomes numerically unstable at sufficiently high numbers. Therefore, at high copy numbers, we solve for the steady state of the original master equation (Eqn. A3) by iteration. This is less efficient but more numerically stable. Specifically, we tile the columns of p nm into one vector P i , and write the master equation as where M ii is the reshaped operator in Eqn. A3, and −D i and N ii are its diagonal and non-diagonal components, respectively. In steady state, we have which we solve iteratively until convergence.

Chromosomal alterations and antibiotic resistance
The native strain was a variant from a previous study [13] and had the following chromosomal alterations and antibiotic resistance:

Locus Construct
Antibiotic Resistance AmyE P hyperspank comK Spectinomycin SacA P comG cfp, P comS yfp Chloramphenicol The SynExSlow strain was modified from [16] as described in the main text and had the following chromosomal alterations and antibiotic resistance:  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 gives smoother results, but 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). Therefore, in Fig. 3C and Fig. 5C, we use the two-Poisson method for f < 1 and use the inflection-point method to determine the transition to f = 1.

Native SynEx
ComK molecule # (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 P (ω) = |ñ(ω)| 2 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. here, whereas in the native circuit they persist beyond this value (see Fig. A5). This confirms the effect seen in the main text that the prevention of damping is more pronounced in the native circuit than in the SynEx circuit. . We see that, as in Fig. A6, noise prevents damping at large values of the control parameter, even at high molecule numbers (column 4). However, as in Fig. A8, noise does not induce repeated excitations at small values of the control parameter (column 2). We conclude that the former effect is more robust.