Skip to main content
Advertisement
  • Loading metrics

Role of Autoregulation and Relative Synthesis of Operon Partners in Alternative Sigma Factor Networks

Abstract

Despite the central role of alternative sigma factors in bacterial stress response and virulence their regulation remains incompletely understood. Here we investigate one of the best-studied examples of alternative sigma factors: the σB network that controls the general stress response of Bacillus subtilis to uncover widely relevant general design principles that describe the structure-function relationship of alternative sigma factor regulatory networks. We show that the relative stoichiometry of the synthesis rates of σB, its anti-sigma factor RsbW and the anti-anti-sigma factor RsbV plays a critical role in shaping the network behavior by forcing the σB network to function as an ultrasensitive negative feedback loop. We further demonstrate how this negative feedback regulation insulates alternative sigma factor activity from competition with the housekeeping sigma factor for RNA polymerase and allows multiple stress sigma factors to function simultaneously with little competitive interference.

Author Summary

Understanding the regulation of bacterial stress response holds the key to tackling the problems of emerging resistance to anti-bacteria’s and antibiotics. To this end, here we study one of the longest serving model systems of bacterial stress response: the σB pathway of Bacillus subtilis. The sigma factor σB controls the general stress response of Bacillus subtilis to a variety of stress conditions including starvation, antibiotics and harmful environmental perturbations. Recent studies have demonstrated that an increase in stress triggers pulsatile activation of σB. Using mathematical modeling we identify the core structural design feature of the network that are responsible for its pulsatile response. We further demonstrate how the same core design features are common to a variety of stress response pathways. As a result of these features, cells can respond to multiple simultaneous stresses without interference or competition between the different pathways.

Introduction

Bacteria survive in stressful environmental conditions by inducing dramatic changes in their gene expression patterns [1,2]. For a variety of stresses, these global changes in gene expression are brought about by the activation of alternative σ-factors that bind the RNA polymerase core enzyme and direct it towards the appropriate stress response regulons [3]. Consequently, to ensure that these σ-factors are only active under specific environmental conditions, bacteria have evolved regulatory systems to control their production, activity and availability [3,4]. These regulatory networks can be highly complex but frequently share features such as anti-σ-factors, partner switching mechanisms and proteolytic activation [4]. The complexity of these networks has impeded a clear mechanistic understanding of the resulting dynamical properties. In this study, we focus on one of the best studied examples of alternative σ-factors, the general stress-response regulating σB in Bacillus subtilis [5] to understand how the structure of the σ-factor regulatory networks is related to their functional response.

The σB-mediated response is triggered by diverse energy and environmental stress signals and activates expression of a broad array of genes needed for cell survival in these conditions [5]. Activity of σB is tightly regulated by a partner-switching network (Fig 1A and 1B) comprising σB, its antagonist anti-σ-factor RsbW, and anti-anti-σ-factor RsbV. In the absence of stress, RsbW dimer (RsbW2) binds to σB and prevents its association with RNA polymerase thereby keeping the σB regulon OFF. Under these conditions most of RsbV is kept in the phosphorylated form (RsbV~P) by the kinase activity of RsbW2. RsbV~P has a low affinity for RsbW2 and cannot interact with it effectively [6]. However, in the presence of stress, RsbV~P is dephosphorylated by one or both of the dedicated phosphatase complexes: RsbQP for energy stress and RsbTU for environmental stress [710]. Dephosphorylated RsbV attacks the σB-RsbW2 complex to induce σB release, thereby turning the σB regulon ON [11]. Notably, the genes encoding σB and its regulators lie within a σB-controlled operon [12], thereby resulting in positive and negative feedback loops.

thumbnail
Fig 1. σB general stress response network.

A. Network diagram of the σB general stress response. The network has two modules: a transcriptional module that inputs the free σB level and outputs the total concentrations of operon proteins RsbV (RsbVT), RsbW (RsbWT) and σB (BT); and a post-translational module that uses RsbVT, RsbWT and BT and the stress phosphatase levels as inputs to output the level of free σB. In the post-translational module, energy and environmental stresses activate the stress-sensing phosphatases RsbQP (QP) and RsbTU (TU) which dephosphorylate RsbV which in turn activates σB by releasing it from the σB-RsbW2 complex. Note only the monomeric forms of RsbW and RsbV have been shown for simplicity. B. Simplified view of the σB network. The σB network works as a feedback loop wherein free level controls RsbVT, RsbWT and BT levels via operon transcription and the three operon components together with the stress level determine the free σB level. The feedback loop sign is can be either positive or negative depending on whether increase in operon component levels impacts free σB level either positively or negatively. C-E. Dynamics of free σB in response to a step-increase in phosphatase concentration for different combinations of the relative synthesis rates of σB operon partners (λW = RsbWT/BT, λV = RsbVT/BT).

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

Recently, it was shown that under energy stress σB is activated in a stochastic series of transient pulses and increasing stress resulted in higher pulse frequencies [13]. It has also been shown that increase in environmental stressor such as ethanol leads to a single σB pulse with an amplitude that is sensitive to the rate of stressor increase [14]. While it is clear that the pulsatile activation of σB is rooted in the complex architecture of its regulatory network (Fig 1A and 1B) its mechanism is not fully understood. Previous mathematical models of the σB network either did not produce the pulsatile response [15] or made simplifications to the network [13] that are somewhat inconsistent with experimentally observed details. As a result, it remains unclear which design features of the σB network enable its functional properties.

To address these issues we develop a detailed mathematical model of the σB network and examine its dynamics to understand the mechanistic principles underlying the pulsatile response. By decoupling the post-translational and transcriptional components of the network we show that an ultrasensitive negative feedback between the two is the basis for σB pulsing. Moreover we find that the relative synthesis rates of σB and its operon partners RsbW and RsbV, plays a critical role in determining the nature of the σB response. We also use our model, together with previously published experimental data from [13,14], to explain how the σB network is able to encode the rate of stress increase and the size of stochastic bursts of stress phosphatase into the amplitudes of σB pulses.

We further develop this model to investigate how the network functions in the context of other σ-factors. As in many other bacteria, σB is one of the many σ-factors that complex with RNA-polymerase core that is present in limited amounts [3,16]. Therefore, when induced these alternative σ-factors compete with one another and the housekeeping σ-factor σA for RNA polymerase. We use our model to investigate how the design of this network enables it to function even in the presence of competition from σA which has a significantly higher affinity for RNA polymerase [17]. Lastly, we investigate how multiple alternative σ-factors compete when cells are exposed to multiple stresses simultaneously. Using our model we identify design features that are ubiquitous in stress σ-factor regulation and critical to bacterial survival under diverse types of stresses.

Results

Biochemically accurate model of σB pulsing

In a recent study, Locke et al. [13] demonstrated that a step-increase in energy stress results in pulsatile activation of σB. The study also proposed a minimal mathematical model of the network which reproduced pulsing in σB. However, this model included several assumptions inconsistent with experimentally observed details: (i) Phosphorylation and dephosphorylation reactions were assumed to follow Michaelis-Menten kinetics despite the fact that kinase (RsbW) and phosphatase concentrations are known to be comparable to substrate (RsbV) concentrations [18] so the approximation breaks down [19], (ii) σB and RsbV are represented as a single lumped variable rather than separate species and, (iii) partner-switching, and the formation and dissociation of various RsbW2 complexes were not included explicitly. Though this minimal model produces pulses resembling their experimental observations, it does not depict a biochemically accurate picture of the σB network. Consequently it cannot be used to uncover the design features that enable σB pulsing.

To understand the σB network response we built on our earlier study [15] to develop a detailed mathematical model that explicitly includes all known molecular interactions in the network. Note that we made one significant change to the model discussed in [15]. The model in [15] assumed that the synthesis rates for σB and its operon partners (RsbW and RsbV) follow the stoichiometry of their binding ratios (i.e. RsbWT/BT = 2 and RsbWT/RsbVT = 1; where BT, RsbWT and RsbVT represent total σB, RsbW and RsbV concentrations respectively). However experimental measurements have shown that σB, RsbW and RsbV are produced in non-stoichiometric ratios [18]. The exact mechanism underlying these non-stoichiometric ratios is currently only incompletely understood. However, analysis of the open-reading frames in the operon showed that rsbV and rsbW may be translationally coupled due to overlapping termination and initiation codons [20] which may ensure that they are expressed in similar amounts. The same analysis also showed that the rsbW and sigB reading frames overlapped and that this overlap was preceded by a region of dyad symmetry which may form a stem-loop structure [20]. These features may interfere with sigB translation and lead to lower expression of σB than its binding partners RsbV and RsbW. To account for these features, in contrast to our earlier study, we assumed σB, RsbW and RsbV can be produced in non-stoichiometric ratios and studied how changes in relative synthesis rates of σB operon partners affect the response of the σB network to step-increases in energy stress phosphatase levels. We note that RsbX, a negative regulator of RsbTU phosphatase [21], is not included in our model. RsbX was excluded for simplicity since it is not essential for the pulsatile response of the σB network [14].

Simulations of this detailed model showed that different combinations of RsbW:σB (λW) and RsbV:σB (λV) relative synthesis rates lead to qualitatively different dynamical responses of the σB network. For operon partner synthesis ratios similar to those estimated in [18], (i.e. RsbWT > 2BT and RsbWT ≈ RsbVT) our model responded to a step-up increase of the phosphatase with a pulsatile σB response (Fig 1C) that resembled the experimentally observed behavior [13]. In contrast, when RsbW:σB and RsbV:σB relative synthesis rates follow the stoichiometry of their binding ratios pulsing is not observed and the σB activity monotonically increases over time (Fig 1D). Pulsing also disappears when RsbW synthesis is high enough to neutralize both its binding partners (Fig 1E).

Pulsing originates from emergent negative feedback in the network

To understand why the pulsatile response is only observed for certain operon partner synthesis rates, we investigated our mathematical model by decoupling the network’s transcriptional and post-translational responses (as shown in Fig 1A). By varying the σB operon transcription rate, while keeping the relative synthesis rates of RsbW:σB (λW) and RsbV:σB (λV) fixed, we were able to calculate the post-translational response (Fig 2A, blue curve) of the σB network: [σB] = Fp ([BT], [PT]). This function describes how the free σB concentration varies as a function of BT (total concentration of σB) and PT (total phosphatase concentration). Note that although we refer only to BT for brevity, RsbWT and RsbVT are always assumed to increase in proportion to the BT for this post-translational response. This post-translational function is analogous to an in vitro assay wherein various combinations of total σB (BT−and proportional amounts of RsbWT and RsbVT) and total phosphatase (PT) are mixed together and then the resulting free σB concentration is measured. In parallel, we calculated the transcriptional response (Fig 2A, black curve) [BT] = FT ([σB]) which analogous to a transcriptional reporter construct in vivo, describes how changes in the free σB concentration affect total σB concentrations (and RsbWT and RsbVT concentrations which are always proportional to BT). In this analysis framework, the steady state of the complete closed loop network can be determined by simultaneously solving the post-translational and transcriptional equations, [σB] = FP ([BT], [PT]) and [BT] = FT ([σB]) at each phosphatase concentration PT. Graphing both functions provided the steady-state solution as their intersection point (Fig 2A, red circle).

thumbnail
Fig 2. Negative feedback drives the pulsatile response of the σB network.

A. Decoupled post-translational (blue curve) and transcriptional (black curve) responses of the σB network for λW = RsbWT/BT = 4, λV = RsbVT/BT = 4.5. σB and BT represent the concentrations of free and total σB. Red circle marks the steady states of the full system. Note that RsbWT and RsbVT are assumed to always increase in proportion to the BT for both post-translational and transcriptional responses. B. Sensitivity of the post-translational response (LGP) to changes in total σB concentration (operon production). C. Representation of the σB pulsatile trajectory in the σB-BT phase plane (green curve). Blue and cyan curves are the decoupled post-translational responses at high and low phosphatase concentrations. Black curve is the transcriptional response. D. (λW, λV) relative synthesis parameter space is divided into regions with positive (Region I), negative (Region II) and zero (Region III) post-translational sensitivity that respectively correspond to an effective positive, negative and no feedback in the σB network. Red and black lines represent the analytically calculated region boundaries λW = 2 + λV and λW = 2(1 + λVkdeg / kk).

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

This decoupling approximation allows us to quantify the sign and strength of feedback in the full model. The effective sign of the feedback in the σB network is given by the sign of the product of the sensitivities of two response functions, i.e. sign ((∂FT / ∂[σB])·(∂FP / ∂[BT])). Since σ-factors function as activators of transcription, FT ([σB]) is a monotonically increasing function of σB (i.e. ∂FT / ∂[σB] > 0). Consequently, the sign of the feedback in the σB network is given by the sign of the sensitivity of the post-translational response to RsbBT (i.e. ∂FP / ∂[BT]). In other words, if increase in the operon production leads to an increase in free σB then the feedback is positive, whereas if increase in the operon production leads to a decrease in free σB then the feedback is negative. Our results show that for the parameters chosen in Fig 1C FP is a non-monotonic function of BT (Fig 2A, blue curve). At low RsbBT, free σB increases as a function of BT because RsbW is sequestered in the W2V2 complex. However at higher BT, the kinase flux dominates the phosphatase flux resulting in an increased RsbV~P and the freeing of RsbW2 from RsbV. Freed RsbW2 sequesters σB in the W2σB complex. Furthermore, in the total σB concentration range where ∂FP / ∂[BT] < 0 in Fig 2B, the post-translational response is quite steep (Fig 2A), i.e. small changes in BT lead to significant decreases in free σB. This ultrasensitivity can be quantified by calculating the slope in logarithmic space, i.e.

This dimensionless quantity characterizes the ratio of relative changes in σB and BT at steady state (Fig 2B). The sign of LGP defines the effective sign of the feedback loop and if the magnitude of |LGP| > 1 defines an ultrasensitive response. For the σB network, in the region around the steady state LGP < −1 indicating that the σB network operates in an ultrasensitive negative feedback regime. Two types of post-translational reactions that are known to produce ultrasensitivity play a role here (S1A and S1B Fig): (1) Zero-order ultrasensitivity due to competition between RsbW kinase and RsbQP/RsbTU phosphatases for RsbV and (2) molecular titration due to sequestration of σB by RsbW. Notably around the steady state, whereas both the fraction of unphosphorylated RsbV and the fraction of free σB decrease ultrasensitively as a function of increase in operon expression (proportional to BT) the latter is far more sensitive (S1C Fig). This indicates that molecular titration between σB and its binding partners may contribute more to the ultrasensitivity of the post-translational response than the zero-order competition between RsbW and stress phosphatases. Irrespective of their relative contributions however, our results show that both mechanisms combine to ensure that near the steady state the σB network operates in an ultrasensitive negative feedback regime.

Notably, negative feedback is one of the few network motifs capable of producing adaption-like pulsatile responses [22]. Moreover, ultrasensitivity of the feedback ensures homeostatic behavior—making the steady state robust to variations of parameters [22]. This explains why in Fig 1C a step-increase in the phosphatase concentration in our model leads to a σB pulse followed by return to nearly the same steady state. Plotting the trajectory of the σB pulse (green curve, Fig 2C) on the ([σB], [BT]) plane and over the post-translational and transcriptional responses (Fig 2C) illustrates the mechanism driving this pulsatile response. Starting at the initial steady state (red circle), an increase in phosphatase shifts the ultrasensitive post-translational response (cyan to blue curve) so that free σB is rapidly released from the RsbW2B complex whereas total σB levels remain relatively unchanged. The increase in σB operon transcription eventually causes accumulation of total σB and the anti-σ-factor RsbW. This in turn forces the σB level to decrease, following the post-translational response curve, to the new steady state (gray circle) which has very little free σB thereby completing the σB pulse.

The same analysis can be applied for different values of relative synthesis rates, i.e. those that correspond to Fig 1D and 1E. As shown in S2 Fig these parameter values do not produce an ultrasensitive non-monotonic post-translational response. Consequently they do not lead to the emergence of overall negative feedback explaining their non-pulsing dynamics. To determine if the presence or absence of negative feedback more generally explains the different dynamical responses in Fig 1C–1E, we sampled different combinations of relative synthesis rates ([RsbWT] / [BT] = λW and [RsbVT] / [BT] = λV) and calculated the post-translational sensitivities. Our calculations showed that based on the sign of post-translational sensitivity (LGP) the relative synthesis parameter space can be divided into three regions (Fig 2D). For (λW, λV) combinations in Region I the sensitivity is always positive. Increase in λW leads the system into an ultrasensitive negative regime (LGP < 0 and |LGP| ≫ 1) in Region II. A further increase in λW or a decrease in λV transitions the system into a non-responsive (LGP ∼ 0) state in Region III. Dynamic simulations for sampled (λW, λV) combinations confirm that pulsatile responses to step-up in phosphatase concentration are restricted to Region II where the effective feedback is negative (S2 Fig).

To understand the boundaries between the three regions and how the level of the phosphatase affects the network, we developed a simplified analytical model that is based on the observation that RsbW and RsbV bind strongly to each other [18] (see S1 Text for details). This approximation allowed us to determine the boundaries in Fig 2D (black and red lines) and resulted in a clear biological interpretation of the three regions. In Region I the amount of RsbW, irrespective of phosphatase level, is insufficient to bind all of its partners and consequently some fraction of σB always remains free or unbound to RsbW. In contrast in Region II, the amount of phosphatase determines how much RsbV is in its inactive phosphorylated form RsbV~P and therefore whether the amount of RsbW is sufficient to bind all of its partners depends on the levels of RsbV~P. As a result, for this region, the ratio of kinase and phosphatase (PT) fluxes determines the post-translational response. Lastly, Region III is the opposite of Region I in that the amount of RsbW is more than sufficient to bind all of its partners, even when all RsbV is unphosphorylated. As a result, irrespective of phosphatase levels, very little σB is free and its level is nearly insensitive to changes in total σB. Thus negative feedback and consequently pulsing are only possible in Region II where changes in phosphatase can shift the balance between the prevalent partner complexes.

The role of negative feedback in producing a pulsatile response also explains why pulsing does not occur in strains where σB operon is transcribed constitutively [13]. In this case, the σB network lacks the negative feedback necessary to produce a pulsatile response. A step-increase in phosphatase still leads to an increase in free σB due to the change in the post-translational response; however, this not followed by an increase in total σB levels (S2C Fig). Consequently, an increase in phosphatase results in a monotonic increase in free σB rather than a pulse (S2F Fig).

The only actual measurements of λW and λV were made by Delumeau et al. [18] using a quantitative western blot assay. Interestingly they report that λW = 2.9, λV = 1.7 in the absence of stress and λW = 2.4, λV = 2.65 in the presence of stress. These measurements suggest that the ratios might change depending on whether cells are under stress. Although the mechanism underlying this change is unclear, our model predicts that both measured ratio pairs lie within the negative feedback regime shown in Fig 2D. Accordingly our simulations show that the network responds to step-increases in phosphatase levels with a pulsatile response for both pre-stress and post-stress (λW, λV) values (S3 Fig). However, due to reduced ultrasensitivity of the system for these parameters, concertation of free σB following increase in the stress (phosphatase) does not perfectly adapt to the pre-stress value (S3 Fig). In an attempt to match the near-perfect adaptation reported in Refs. [13,14] we’ve chosen to do further analysis with λW = 4 and λV = 4.5.

Notably, our simulations also showed that it is not essential for the phosphatase level to remain fixed after a stress-induced step-increase. In fact, we found that a dilution mediated decline in phosphatase level post-step-increase has little impact on the pulse amplitude (S4 Fig). This observation can be explained by the relatively rapid dynamics of the post-translational response as compared to the gradual nature of dilution and suggests that pulsatile dynamics are relevant even for experimental conditions where phosphatase levels do not remain fixed in stressful conditions [14,23].

Further our decoupling method also sheds light on another experimental observation by Locke et al. [13]: the dependence of σB pulse amplitude on the phosphatase level. Specifically, we found that σB pulse amplitude is a threshold-linear function of the phosphatase concentration (S5 Fig). Our decoupling method shows that this threshold-linear behavior arises because the σB network only operates in a negative feedback regime for phosphatase concentrations higher than a threshold. Below the phosphatase threshold, the post-translational response [σB] = FP ([BT], [PT]) ∼ 0 and is insensitive to RsbBT (S5B and S5C Fig). Thus, the full system lacks the negative feedback and as a result σB does not pulse. Using our analytical approximation we found that this phosphatase threshold is proportional to the basal level of RsbW kinase synthesis rate and the ratio of the kinase and phosphatase catalytic rate constants (S5D and S5E Fig). Increase in the basal σB operon expression rate increases the phosphatase threshold. Further, an increase in the relative synthesis rate of RsbW (λW = [RsbWT] / [BT]) makes the phosphatase threshold more sensitive to the σB operon expression rate, whereas a decrease in ratio of the kinase and phosphatase catalytic rate constants makes it less sensitive (S5D and S5E Fig). This shows that the phosphatase threshold represents the concentration at which the phosphatase is able to match the basal kinase flux.

Altogether these results show how the ultrasensitive negative feedback plays a critical role in determining many properties of the σB network pulsatile response and how the decoupling method can facilitate the identification of essential design features that enable the existence of this negative feedback.

Under energy stress conditions σB network encodes phosphatase burst size into pulse amplitudes

In the preceding sections we have shown how the σB network responds to a step-increase in RsbQP or RsbTU phosphatases by producing a single pulse of activity. However, Locke et al. [13] have shown that an increase in energy stress leads to a sustained response with a series of stochastic pulses in σB activity. This study further showed that this sustained pulsing response is driven by noisy fluctuations in level of energy-stress-sensing phosphatase RsbQP. While the mean level of RsbQP is regulated transcriptionally by energy stress, its concentration in single cells can fluctuate due to the stochasticity of gene expression [8,13]. To determine if our model could explain this response to stochastic fluctuations in RsbQP, we modified it to include fluctuations in the concentration of this phosphatase.

Based on previous theoretical [24,25] and experimental [26] studies we assume that fluctuating phosphatase level follows a gamma distribution which is described by two parameters—burst size (b, average number of molecules produced per burst) and burst frequency (a, number of bursts per cell cycle). The mean phosphatase in this case is the product of burst size and burst frequency (〈PT〉 = ab). Thus, energy stress can increase mean phosphatase by changing burst size or burst frequency or both. In other words, stress conditions can increase phosphatase levels by either producing more phosphatase molecules per transcription-translation event or by making these events more frequent. While the results of [13] cannot exclude either mechanism, we can use our model to uncover which mechanisms is dominant.

First, we performed stochastic simulations in which mean phosphatase concentration was varied by changing burst size. These simulations reproduced all the experimentally-observed features of the σB pulsatile response. Specifically our results show that stochastic bursts in stress phosphatase levels lead to pulses of σB activity (Fig 3A). Moreover, consistent with the experimental observations of [13], our model showed that the amplitude of σB pulses increases linearly with the stress phosphatase level (Fig 3A and 3B). Finally, we found that stress-mediated increases in phosphatase concentration lead to an ultrasensitive (effective Hill coefficient ~5.6) increase in the frequency of σB pulsing (Fig 3C) and an ultrasensitive (effective Hill coefficient ~2) increase in the level of σB target expression (Fig 3D).

thumbnail
Fig 3. Pulsatile response of the σB network to stochastic phosphatase bursts during energy stress.

Model simulations for σB network response where energy stress leads to an increase in stress-sensing phosphatase RsbQP burst size (A-D) or RsbQP burst frequency (E-H). A,E. Simulations show stochastic bursts in levels of RsbQP lead to pulses of σB target promoter activity. Light and dark green curves are sample trajectory from stochastic simulation at high and low stress respectively. Note that σB target promoter activity pulse amplitude increases significantly with increasing stress for burst size modulation (A) but not for burst frequency modulation (E). B,F. Mean σB pulse amplitude increases linearly as a function of mean phosphatase level for burst size modulation (B) but is insensitive to mean phosphatase level for burst frequency modulation (F). Green circles and errorbars show means and standard deviations calculated from stochastic simulations. Black line is a linear fit. C,G. With increasing mean phosphatase level, mean σB pulse frequency increases ultrasensitively for burst size modulation (C) and linearly for burst frequency modulation (G). Green circles and errorbars show means and standard deviations calculated from stochastic simulations. Black curves are a Hill-equation fit with nHill = 5.6 in (C) and a linear fit in (G) respectively. D,H. Mean σB target expression increases ultrasensitively as a function of mean phosphatase level for both burst size (D) and burst frequency (H) modulation. Green circles are the mean σB target expression calculated from stochastic simulations. Black curve is a Hill-equation fit with nHill = 2 in (D) and in nHill = 1.2 (H).

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

Next, we compared these results with stochastic simulations in which burst frequency was modulated (Fig 3E–3H). These simulations also led to an increase in σB pulsing (Fig 3E) and a non-linear increase in the level of σB target expression as mean phosphatase level was increased with more frequent bursts (Fig 3H). However, we found that σB pulse amplitude remains constant for burst frequency modulation (Fig 3E and 3F) unlike the ~5-fold increase for burst-size modulation (Fig 3B). Moreover, the frequency of σB pulses increase linearly with phosphatase level unlike the non-linear increase observed with burst-size-increase simulations (compare Fig 3C and 3G). Notably the experimental observations reported in [13] show that σB pulse amplitude does increase (~3-fold) with an increase in energy stress thus suggesting that increase in phosphatase concentration at high stress is primarily the result of increase in burst size.

To further reinforce the role of mean burst-size modulation in controlling the σB pulsatile response we next examined the cumulative histograms of pulse amplitudes at different phosphatase concentrations. These histograms carry different signatures for burst-size or burst-frequency encoding. The distribution of pulse amplitudes is unchanged with increase in burst frequency (S6A Fig) because σB pulse amplitude is determined by phosphatase burst size and not burst frequency. In contrast, if phosphatase levels are controlled by changing mean burst size then the distribution of pulse amplitudes changes accordingly. Consequently, the normalized cumulative histograms of pulse amplitudes overlap for burst-frequency encoding (S6A Fig) but not burst-size encoding (S6B Fig). Applying this test to the data from [13], we found that the normalized cumulative pulse amplitudes histograms do not overlap (S6C Fig). These results predict that stress affects the σB network via burst-size modulation of phosphatase production which is then encoded into σB pulse amplitudes. While the molecular mechanism that introduces energy stress to the network is still not fully understood, our prediction places an important constraint on it.

σB network encodes rate of environmental stress increase into pulse amplitudes

Our model can also be used to study the response of σB network to environmental stress. Unlike the energy stress phosphatase, the environmental stress phosphatase RsbU is regulated post-translationally by binding of RsbT [2729]. RsbT is trapped by its negative regulators under unstressed conditions but is released upon stress. Consequently, the concentration of RsbTU complex is tightly controlled at the post-translational level and is therefore expected to be relatively insensitive to gene expression fluctuations but sensitive to the level of environmental stress. As a result, step-up increases in environmental stress agents like ethanol produce rapid increases in RsbTU and result in only a single pulse of σB activity [14]. However it has been shown that for gradual increases in stress, σB pulse amplitude depends on the rate of stress increase [14]. To explain this response, we modeled gradual stress with ramped increase in RsbTU complex concentration (Fig 4A). Our simulations showed that the detailed model of σB network is indeed able to capture the effect of rate of stress increase on σB pulse amplitudes. Specifically for a fixed increase in RsbTU complex, the pulse amplitude decreases non-linearly as a function of the duration of phosphatase ramp (Fig 4B and 4E).

thumbnail
Fig 4. Rate sensitivity of the σB pulsatile response to environmental stress.

A. Ramped increases in RsbTU complex concentration were used as model inputs to simulate different rates of stress increase in σB network. B. σB pulse amplitudes in the wildtype model (kdeg = 0.72 hr-1 is the degradation rate of σB operon proteins) resulting from the ramped increases in phosphatase concentration shown in (A). C,D. σB pulse amplitudes resulting from the ramped increase in phosphatase concentration shown in (C) for various degradation/dilution rates (D). E. Non-linear dependence σB pulse amplitude on phosphatase ramp duration for various degradation/dilution rates. Circles and solid curves represent simulation results and Hill-equation fits respectively. Colors represent different kdeg values as in (D). F. Kramp, the half-maximal constant of the non-linear dependence of amplitude on ramp duration, as a function of kdeg.

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

We hypothesized that this ramp rate encoding is the result of the timescale separation between the fast post-translational and the slow transcriptional responses of the σB network. During the pulsed σB activation, post-translational response is rate-limited by the phosphatase ramp. In contrast, the transcriptional response is slow and its rate is set by the degradation rate of σB operon proteins. Following a step-increase in phosphatase, the fast post-translational response ensures that σB reaches its post-translational steady state before the slow increase in RsbW sequesters σB and turns off the pulse (Fig 4A and 4B). However, for a ramped increase in phosphatase the post-translational increase in σB is limited by the rate of phosphatase ramp. This allows RsbW to catch up and terminate the σB pulse earlier, thereby decreasing the pulse amplitude. To test this, we varied the degradation rate of σB operon proteins and proportionally changed the operon transcription rate to ensure that the total concentrations of σB, RsbW and RsbV are kept fixed. We found that indeed pulse amplitude decreases with increase in degradation/dilution rate (Fig 4C and 4D). Our simulations showed that Kramp, the half-maximal constant for the dependence of pulse amplitude on ramp duration, was indeed sensitive to the degradation rate (Fig 4E and 4F). This suggests that the timescale separation between the post-translational and transcriptional responses is the basis of ramp rate encoding into pulse amplitude.

The design of the σB network enables it to compete with σA for RNA polymerase

The results thus far indicate that σB network functions in the effectively negative feedback regime where increase in the operon expression decreases σB activity. Negative feedback loops have been shown to increase the robustness of the system to perturbations. We therefore decided to investigate how the σB network design affects its performance when it faces competition for RNA polymerase from other σ-factors, e.g. from the housekeeping σ-factor σA [16,30,31]. Since σA has a much higher affinity for RNA polymerase [17], a small increase in σA can dramatically increase the amount of σB necessary to activate the transcription of the σB regulon. Thus, changes in σA can alter the input-output relationship of a stress-response σ-factor like σB (S7A and S7B Fig) and thereby adversely affect the survival of cells under stress.

To understand how the σB network handles competition for RNA polymerase, we expanded our model to explicitly include σA, RNA polymerase (RNApol) and its complexes with both σ-factors. The presence of σA will affect transcriptional activity of σB but not post-translational interactions between σB operon partners (Fig 5A, left panel). Therefore, post-translational response [σB] = FP ([BT], [PT]) is not affected by σA. In contrast, in the transcription response, an increase in σA decreased the ‘effective affinity’ of σB for RNApol and consequently higher levels of free σB are necessary to achieve the same production rate for σB target genes.

thumbnail
Fig 5. Negative feedback insulates the σB response from competition with houskeeping σ-factor σA.

A. Simplified network diagrams of stress σ-factor σB competing with housekeeping σ-factor σA for RNA polymerase. In all cases, a σB phosphatase controls the stress-signal driven activation of σB. (B,C). Trajectories of free σB (B) and σB target promoter activity (C) in response to stochastic phosphatase input for both networks at two different levels of σAA = 9μM—low competition-regime and σA = 12μM—high-competition regime for RNA polymerase). D-E. Mean free σB concentration (D) and mean σB target promoter activity (E) as a function of total σA concentration (AT) for both networks in (A) at fixed mean phosphatase (mean PT = 0.5 μM). Gray vertical line shows the total RNA polymerase level which was fixed at 10 μM.

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

Using our model, we examined how changes in σA level affect the network response to energy stress signal, i.e. under stochastically fluctuating RsbQP phosphatase levels. Our simulations showed that phosphatase bursts lead to pulses of free σB and pulsatile transcription of σB-controlled promoters (Fig 5B and 5C) as the presence of σA does not affect the effective feedback sign. Notably our results also showed that the amplitudes of σB target promoter pulses are hardly affected by a ~30% increase in σA (Fig 5C, left panel). This surprising insensitivity of the phosphatase-σB target dose-response to RNApol competition is the result of the ultrasensitive negative feedback between free σB and total σB. Due to the ultrasensitivity of this feedback, a small decrease in total σB levels resulting from the increase in σA causes a large increase in σB pulse amplitude (Fig 5B left panel, Fig 5D green line). This increased amplitude compensates for the increased competition for RNApol and insulates the network from perturbations (Fig 5D and 5E, green curves).

To further illustrate the importance of the negative feedback in insulating the network, we compared the response of the wildtype network to an “in silico” mutant network wherein the σB operon is constitutive rather than σB dependent (Fig 5A). Consequently this network lacks any feedback between free σB and total σB. Our simulations (Fig 5B, right panel) show that the free σB concentration of the no-feedback-network does not show adaptive pulsing and therefore σB concentration fluctuates along with the phosphatase levels. Increase in σA did not affect this response. This is expected since in the absence of feedback σA only affects the expression of σB targets in this network (Fig 5A, right panel). Without an increase in free σB (Fig 5D), the increased competition for RNApol at higher σA reduced the σB target promoter activity (Fig 5C and 5E). Similarly a positive feedback network design is also incapable of increasing free σB in response to an increase in σA (S5C, S5D and S5E Fig). Thus fluctuations in σA can interfere with the σB stress-response of these alternative network designs. In contrast, the wildtype σB network with its ultrasensitive negative feedback design can compensate for competition effects (Fig 5D and 5E).

Negative feedback designs of stress-response σ-factor networks minimizes interference

The emergent negative feedback design of the network discussed here is not unique to σB. Transcription of many alternative σ-factors in B. subtilis as well in other bacteria is often positively auto-regulated but sigma-factor operons often include post-translational negative regulators [3,12,3235]. For example σW, a σ-factor in B. subtilis that controls the response to alkaline shock [36] is co-transcribed with its anti-σ-factor RsiW. In the absence of stress, RsiW sequesters σW in an inactive complex. σW is activated by stress signals which trigger the cleavage and degradation of RsiW thereby releasing and activating σW target expression [37]. Although it is unknown whether the σW network functions in a negative feedback regime similar to σB or if it pulses, it is possible for this network to exhibit these design properties. If RsiW is expressed in stoichiometric excess of its binding partner σW from the σW-regulated operon which they share [38], then similar to the σB network, σW would operate in a negative feedback regime.

To determine if negative feedback control offers any advantages when multiple stress σ-factors are active, we built a new model that includes three σ-factors: σB, σW and σA (Fig 6A). In this model we use σB and σW to denote two generic stress σ–factors and accordingly anti-σ-factors RsbW (RsiW) and other details of post-translational regulation were excluded for simplicity. Thus our results apply to any combination of alternative σ-factors with ultrasensitive negative feedback control. In this general model, regulation of free σB and σW was modeled with simplified identical versions of the negative feedback design of the σB network (S7A Fig). Under this simplification, free σB and free σW are non-monotonic functions of their respective total concentrations BT and WT. These non-monotonic functions are qualitatively similar to the post-translational response function shown in Fig 2B and depend on a signaling proteins PB (for σB) and PW (for σW). Following the previous section, this model explicitly includes σA, RNApol and its complexes with σ-factors. As a result, transcriptional activity of both σB and σW depend on σA and RNApol concentrations (see S1 Text). Concentrations of RNApol and σA were chosen to ensure that amount of RNApol is insufficient to bind to all σ-factors at the same time. All other parameters of the simplified model were chosen to approximately match the full σB network model and ensure that both σB and σW operate in the negative feedback regime. Consequently for the chosen parameters this simplified model acts like our detailed model and responds to step increases in the stress signaling protein PB (or PW) by producing a pulse of σB (or σW) activity (S7C and S7D Fig). To enable a comparison of the competition between σ-factors for different types of feedback we hereafter focus on only steady state response, however our conclusions are also valid for the averaged pulsatile dynamical responses that could be characteristic of the negative feedback σ-factor networks.

thumbnail
Fig 6. Negative feedback minimizes competition between stress σ factors for RNA polymerase.

A,B. Simplified network diagrams of stress σ-factors σB and σW and housekeeping σ-factor σA competing with each other for RNA polymerase. σB and σW activities are regulated by negative and positive feedbacks in (A) and (B) respectively. In both cases, signaling proteins PB and PW control the stress-signal driven activation of σB and σW respectively. C, D. Dependence of free σB and σW levels on PB at fixed PW (= 2μM). In the wildtype negative feedback system (C), increase in σB phosphatase leads to an increase in both free σB (green curve) and free σW (red curve). In the positive feedback system (D), increase in σB phosphatase leads to an increase in free σB (green curve) and a decrease in free σW (red curve). E, F. σB and σW target promoter activities as a function of PB at fixed PW in the wildtype negative feedback system (E), and the positive feedback system (F). G, H. RNA polymerase bound σB (Rpol-σB) as a function of PB at fixed PW in the wildtype negative feedback system (G) and the positive feedback system (H). Increase in σB phosphatase (PB) leads to an increase in Rpol-σB (green curve) and corresponding decreases ΔRpol-σW (core complex with σW, red area) and ΔRpol-σA (complex with σA, blue area).

https://doi.org/10.1371/journal.pcbi.1005267.g006

We used this simple model to study the response when cells are simultaneously exposed to multiple stresses creating competition for RNApol. For these simulations we fixed σA levels and studied how activation signals for one alternative σ-factor affects the activity of another. As before (S7A and S7B Fig), increased availability of one stress σ-factor leads to a competition for RNA polymerase and as a result reduces the activity of another stress σ-factor (S8E and S8F Fig). However, when negative feedback loops are present, surprisingly, increasing the stress signal for one σ-factor did not lead to any significant change in the activity of another σ-factor. For example, increasing stress signaling protein PB while keeping PW fixed leads to an increase in free σB but also results in a small increase in free σW (Fig 6C). This response can be explained by the ultrasensitive negative feedback loops controlling the two stress σ-factors. An increase in free σB by stress signaling protein PB leads to increased competition for RNApol resulting in a decrease in the production of RsbW. But since σW is regulated by a negative feedback, a decrease in total RsbW concentration actually frees up more σW thereby insulating σW target activity from the effects of RNApol (Fig 6E). Similarly the dynamic response of the stress σ-factors is also insulated from competition and an increase in fixed PW levels increases the pulse amplitude of σB in response to step changes in stress signaling protein PB (S8A–S8D Fig). This compensation of changes in RNA polymerase availability comes about because both σB and σW are regulated by ultrasensitive negative feedbacks in our model. As a result of this negative feedback, both σ-factor networks function as homeostatic modules. Homeostatic resistance to changes in signals is an intrinsic property of ultrasensitive negative feedback motifs.

Thus the two stress σ-factors are able to function simultaneously despite the scarcity of RNApol. The mechanism minimizing competition between stress σ-factors becomes clearer when we track the changes in σ–RNApol complexes as a function of the stress signaling protein PB. As PB increases, more free-σB becomes available and binds to RNApol (Fig 6G). However this RNApol must be accounted for by the RNApol lost by the other operating σW and σA factors. Comparing the contributions of each σ-factor shows that despite the fact that σA has a much higher affinity for RNApol, most of the RNApol in the σB-RNApol complex is drawn from the σA-RNApol pool rather than σW-RNApol pool (Fig 6G). Thus the negative feedback design allows stress σ-factors to minimize their competition with each other at the expense of the housekeeping factor σA.

The role of the negative feedback in producing this response becomes clear when we compare the response of an “in silico” mutant network with positive feedback loops between σB and BT and σW and WT (Fig 6B). These positive feedback loops are expected to display no homeostatic properties and as a result, in this network activation of σB should significantly decrease σW activity. Indeed, our simulation for the positive feedback network (Fig 6D) demonstrates that with increase in stress signaling protein PB and the resulting increase in free σB, the free σW concentration decreases. As a result of the increased competition for RNApol and the decreased free σW, σW target promoter activity in this network decreases as a function of PB (Fig 6F). Moreover comparing changes in σ–RNApol complexes as a function of stress signaling protein PB we find that most of the RNApol in the σB-RNApol complex is drawn from the σW-RNApol pool rather than σA-RNApol pool (Fig 6H). Thus the negative feedback designs are essential for stress σ-factors not only to tolerate competition from σA, but also to avoid competing with each other when the cell is simultaneously exposed to multiple types of stresses.

Discussion

Taken together, our results show how the design of the σB network includes an implicit ultrasensitive negative feedback that plays multiple functional roles. This design enables pulsatile activation of σB in response to energy stress and rate-sensitivity to increases in environmental stress. Moreover, our model predicts that the same design feature allows the network to effectively compete with house-keeping and other alternative σ-factors for RNA polymerase core.

Prompted by recent observations of the highly dynamic pulsatile response of the σB network [13,14], we have developed a mathematical model that reproduces all reported features of the response including pulsatile activation in response to stress. Our model avoids making ad hoc simplifications and instead captures all the known molecular details of the network. By decoupling the post-translational and transcriptional responses in our model we were able to derive a simplified view of the network that illustrates how the pulsatile response is mechanistically based on the ultrasensitive negative feedback in the network. Using this method we identified the ratios of σB, RsbW and RsbV synthesis rates as the most critical design property, which by controlling the post-translational response determines the sign of the feedback in the network as well as all qualitative features of the network response. This highlights how ignoring non-transcriptional interactions and focusing on transcriptional regulatory interactions alone can be misleading when trying to identify or characterize network motifs. Notably, recent analyses of networks like bacterial two-component systems [39] and the sporulation phosphorelay [40] have similarly shown how the effective sign of feedback in these networks depends critically on their post-translational interactions.

Interestingly previously reported measurements of the relative amounts of σB, RsbW and RsbV using a quantitative western blot assay showed that the ratios between these components may change depending upon stress [18]. Our model shows that for both pre- and post-stress values of these synthesis ratios, the σB network lies within the negative feedback regime and responds to step-increases in phosphatase levels with a pulsatile response. Thus these stress induced expression changes do not affect our conclusions about the function of the σB network. Nevertheless their mechanistic basis, whether transcriptional or post-translational, remains unclear and may add an additional layer of complexity to regulation. Together with our modeling results these observations highlight the need for more quantitative experimental methods to determine the relative synthesis rates of σB operon components.

The decoupling of the post-translational and transcriptional response greatly facilitated the identification of critical design features despite the complexity of the network. This separation greatly reduces the dimensionality of the dynamical system by enabling an independent input–output analysis for the two modules. Similar methods have also been applied to deduce core functional properties in other bacterial networks comprising two-component systems and alternative σ-factors [4143]. Interestingly our analysis revealed that the post-translational and transcriptional module structures of the σB network and the phosphorelay controlling B. subtilis sporulation are remarkably similar [40]. Despite the differences in molecular details, in both networks increase in total transcription factor levels produces a non-monotonic response in the active transcription factor. Combining this response with the transcriptional feedback produces an ultrasensitive negative feedback in both networks. The relevance of these similarities is evidenced by the fact that both networks produce dynamically similar pulsatile responses even though they are activated by entirely different stimuli.

We further used our model to establish that energy stress controls σB pulses frequency by modulating the size of stochastic bursts of energy stress phosphatase rather than burst frequency. We reached this conclusion by showing that only burst size modulation can explain the experimentally observed changes in mean pulse amplitude and pulse amplitude distribution with increasing energy stress level. It should be noted that although our results provide clear evidence in favor of phosphatase burst-size modulation, direct confirmation of this mechanism necessitates the use of single-molecule techniques such as RNA-FISH [44,45] that can be used to estimate stochastic properties of gene expression. This result raises the question whether pulsatile σB response can achieve proportional expression of downstream genes, as was previously suggested [13,46]. This proportional control requires the distribution of pulse amplitudes to remain fixed even as stress levels increase. However under the burst-size encoding strategy, pulse amplitude distributions change as stress levels increase thereby negating the efficacy of a pulsed response in producing proportional expression of downstream genes.

The functional significance of pulsatile response may instead lie in its ability to encode the rate of environmental stress increase. Our model showed that this rate encoding follows from the timescale separation between the fast post-translational and the slow transcriptional responses in the network. As a result cells are able to encode the rate of stress increase into σB pulses. This rate responsiveness is only possible with adaptive pulsatile responses and thus may explain the need for σB pulsing to control the general stress response.

We also used our model to understand the response when placed in the larger context of other σ-factor networks and competition for RNA polymerase. Our results show how the network design is uniquely suited to insulating its response from RNA polymerase competition from the housekeeping σ-factor. Finally we demonstrated how ultrasensitive negative feedback, a ubiquitous feature of stress σ-factor regulation enables different stress σ-factors to operate simultaneously without inhibiting each other. These results are relevant not only for understanding the stress response of bacteria but also increasingly for the design of synthetic circuits. The movement towards the construction of larger genetic circuits has produced numerous recent designs that include multiple independent modules that rely on shared resources or actuators to function [4749]. Our results highlight how competition between modules for shared resources can significantly affect the performance of these synthetic circuits. Further, inspired by the design of naturally occurring stress σ-factor network we provide new design rules that can improve the performance and robustness of the synthetic networks.

Methods

The details of all biochemical reactions in the model and the corresponding differential equations are described in the S1 Text.

Mathematical model of σB stress-response network

Our mathematical model of σB network is based on a previous model proposed in [15]. This ODE-based model explicitly includes all known molecular species, post-translational reactions and the transcriptional regulation of the σB operon by σB. Below we formulate the set of reactions and associated differential equations.

Model reactions

The events shown in Fig 1A can be described by the following set of biochemical reactions:

  • Dimerization of anti-σ-factor RsbW (1)
  • Reversible binding of the anti-anti-σ-factor RsbV to anti-σ-factor dimer RsbW2 to form the complexes RsbW2-RsbV and RsbW2-RsbV2 (2) (3)
  • Phosphorylation of the anti-anti-σ-factor RsbV by RsbW2 (4) (5)
  • Reversible binding of σB to RsbW2 to form the complex RsbW2B (6)
  • Reversible displacement of σB by RsbV in the complex RsbW2 (7)
  • Dephosphorylation of phosphorylated anti-anti-σ-factor RsbV~P (8)
  • Protein degradation/dilution due to cell growth (9) where X is any protein or complex in the σB network. For simplicity the degradation/ dilution rate of all proteins and complexes are assumed to be equal.
  • Production of σB, RsbW and RsbV (10) σB, RsbW and RsbV were assumed to be synthesized proportionally as all three are part of the same operon. λW and λV are respectively the ratio of synthesis rate of RsbW and RsbV to σB. Because of σB autoregulation, synthesis was modeled as a hyperbolically increasing function of σB concentration, [σB]

Here ν0 is the basal synthesis rate, f is the fold change in protein synthesis due to positive autoregulation and K is the equilibrium dissociation constant for the binding of σB to the promoter DNA.

The stress signals were assumed to control the concentrations of stress phosphatases RsbTU and RsbQP. For RsbQP, energy stress was assumed to regulate the transcription rate of the phosphatase and the phosphatase concentration was assumed to be subject to stochastic fluctuations resulting from gene expression noise. In contrast, RsbTU concentration is regulated by environmental stress post-translationally [2729], consequently RsbTU concentration was assumed to be stress-dependent but not subject to stochastic fluctuations.

Model equations

We assume mass-action kinetics for all the above reactions (Eqs 110) to obtain the following set of equations that describe network dynamics:

Here [σB] is the concentration of free σB; [W2] is the concentrations of dimeric RsbW; [V] and [VP] are the concentrations of unphosphorylated and phosphorylated RsbV; [W2σB], [W2V], [W2V2] and [VPP] are the concentrations of the corresponding protein complexes. [BT], [RsbWT], [RsbVT] and [PT] are the concentrations of total σB, RsbW, RsbV and phosphatase: (11)

All model parameters are summarized in Table 1.

thumbnail
Table 1. List of parameters values used in the model for σB network.

https://doi.org/10.1371/journal.pcbi.1005267.t001

Additional reactions for the model of competition between σB and σA

To model the competition for RNA polymerase between σB and the housekeeping σ-factor σA (Figs 5 and S5), we extended the model described above and supplemented reactions (1–9) with the following reactions for σA, RNA polymerase (RNApol) and σ–RNApol binding:

  • Reversible binding of σ-factors and RNA polymerase
  • Reversible binding of RNApol-σB complexes to target promoters
  • Production of σB, RsbW and RsbV Where v0 is the basal synthesis rate and νB = ν0f / [pB]tot is the maximal rate. f is the fold change in protein synthesis due to positive autoregulation and [pB]tot is the total concentration of the σB promoter.

To investigate the competition between σB, σW and σA, we used a phenomenological non-monotonic function to model the post-translational regulation of stress σ-factors (σB and σW; see S1 Text for details).

Calculation of steady state post-translational and transcriptional responses

The decoupled transcriptional and post-translational responses of the network at steady state were calculated using the bifurcation package MATCONT[50]. The post-translational response [σB] = FP ([BT],[PT]), was calculated by varying the rate of operon transcription while keeping the component synthesis rates (λW, λV) and the total phosphatase concentration (PT) fixed. Similarly, the transcriptional response [BT] = FT ([σB]), was calculated by varying the free σB concentration as an independent variable to calculate the total concentrations of σB, RsbW and RsbV.

Simulations

The parameter values for reversible binding and phosphorylation reactions were taken from [15] or were analysis driven to obtain pulsing in σB. All the parameters used in the model are summarized in Table 1. In the deterministic set-up (Figs 1, 2, 4 and 6, S1S5 and S8) the system of differential equations was solved using standard ode15s solver in MATLAB(The Mathworks Inc., Natick, MA). For stochastic simulations in Figs 3, 5, S6 and S7, the time-varying total phosphatase level [PT] (= [P] + [VPP]) was pre-computed using a gamma distributed Ornstein-Uhlenbeck process as in [13]. This gamma distributed Ornstein-Uhlenbeck process permits independent modulation of mean burst size (b) and frequency (a) [51]. For each phosphatase level, 50 simulations were performed each lasting 10 hours. Pulses were detected by examining local maxima and minima of the simulated trajectories, and subsequently this information was used to compute statistics for pulse amplitude and frequency.

For the simulations of the effect of competition for RNA polymerase (Figs 5 and S7), the total housekeeping σ-factor concentration was varied between 5 and 15 μM. In these simulations we used λW = 4, λV = 4.5 and λW = 2, λV = 2 to simulate the wildtype (negative feedback) and positive feedback networks respectively. For the simulations of the no feedback network we used (λW = 4, λV = 4.5) and f = 0 and v0 = 8.64 μMhr-1 to model the σB–independent constitutive production of operon components.

For the simulations of the competition between σB, σW and σA (Figs 6 and S8), the total housekeeping σ-factor concentration was kept fixed at 12 μM. The post-translational response of stress σ-factors σB (and σW) in this model was described using the following phenomenological function: σBfree = BT / (1 + ([BT] / KB)nb / [PB]mb). Here the constants nb and mb describe the non-linear dependence of free σ on total-σ. We used (nb = 7, mb = 5) and (nb = 0, mb = 3) to simulate the wildtype (negative feedback) and positive feedback networks respectively. KB and KW were fixed at 5μM for simulations of both networks.

Supporting Information

S1 Text. Supplementary Methods and Derivations.

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

(PDF)

S1 Fig. Ultrasensitive negative feedback in the σB network.

A. Decoupled post-translational (blue curve) and transcriptional (black curve) responses of the σB network for λW = RsbWT / BT = 4, λV = RsbVT / BT = 4.5. σB and BT represent the concentrations of free and total σB. Gray circle marks the steady states of the full system. Red and blue lines represent the piecewise analytical approximations of the post-translational response. B. Decrease in the fraction of phosphorylated RsbV (VP + VPP–orange curve) and unbound σB (green curve) as a function total operon expression level according to the post-translational response. C. Sensitivity of the post-translational response of phosphorylated RsbV (Vp+VpP—orange curve) and unbound σB (green curve) to changes in total operon expression level (BT). At the shown steady state (gray circles) both responses have LG<-1.

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

(PDF)

S2 Fig. σB does not pulse for networks that lack negative feedback.

A-C. Decoupled post-translational and transcriptional responses of σB networks that lack negative feedback. (A) λW = 2, λV = 2 (Region I in Fig 2D)—positive feedback system (LG>0); (B) λW = 8, λV = 4.5 (Region III in Fig 2D) a non-responsive system (LG~0); (C) λW = 4, λV = 4.5 with no transcriptional feedback—no feedback system (LG = 0). In each panel cyan and blue curves show the post-translational response at low and high phosphatase concentrations, and black curve shows the transcriptional response. Gray and black circles mark the steady states of the full system. Step-increase in phosphatase causes a shift in the post-translational response from low phosphatase-cyan to high phosphatase-blue and leads to an increase in σB (green curve) in all three systems. D-F. Time-course representations of green trajectories described in A-C. Note that σB does not pulse in any of the three systems.

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

(PDF)

S3 Fig. σB pulse response for λW, λV calculated from the results of Delumeau et. al. [18].

Dynamics of free σB in response to a step-increase in phosphatase concentration for the two different ratios measured in Delumeau et. al. [18]. (A) Pre-stress: λW = 2.9, λV = 1.7 (B) Post-stress: λW = 2.4, λV = 2.65.

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

(PDF)

S4 Fig. Phosphatase level decrease after step-up does not affect σB pulsing.

A. Dynamics of free σB in response to a step-increase in phosphatase concentration. B. Dynamics of free σB in response to a step-increase in phosphatase concentration which is followed by decay due to dilution.

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

(PDF)

S5 Fig. Dependence of σB pulse amplitude on phosphatase concentrations and post-translational parameters.

A. Time-course representations of σB pulse trajectories for small (0.1μM-orange curve) and large (0.4μM-blue curve) step-increases in phosphatase. λW = 4, λV = 4.5 for both trajectories. B,C. Representation of the σB pulse trajectories and decoupled post-translational and transcriptional responses of σB network for small (B) and large (C) step-increases in phosphatase. Cyan and green curves show the post-translational responses at initial and final phosphatase levels. Black curves show the transcriptional response. Black and gray circles mark the steady states of the full system. Note that at the initial phosphatase level the σB~0 and BT is at the basal level of σB operon transcription. The small step-increase in phosphatase does not significantly shift the post-translational response around the initial steady state leading to minor, transient increase in σB (orange curve in B). The large step-increase in phosphatase (C) does significantly shift the post-translational response around the initial steady state leading to prominent pulse in σB (blue curve in C). D. σB pulse amplitudes show a threshold linear response to increase in phosphatase level. The threshold phosphatase level increases with increasing basal level of σB operon transcription (Basal BT). E,F. Phosphatase threshold for pulsing as a function of Basal BT for different values of the (E) RsbW relative synthesis rate (λW) and (F) the ratio of phosphatase to kinase rates (kp / kk). The circles represent threshold levels calculated from simulations. The black lines represent the analytical approximation: PT = ν0kk (λW / 2–1 –λWkdeg / kk) / (kp / kdeg), where v0 is the basal rate of σB operon transcription and kdeg is protein degradation/dilution rate. Basal BT = ν0 / kdeg.

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

(PDF)

S6 Fig. Pulsatile response of the σB network encodes phosphatase burst size not burst frequency.

A-C. Normalized pulse amplitude cumulative histograms for stochastic simulations with (A) burst frequency modulation, (B) burst-size modulation and (C) experimental data taken from [13]. Different colors represent varying levels of mean phosphatase (PT) in the model or mycophenolic acid (MPA, energy stress) in experiments.

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

(PDF)

S7 Fig. Sensitivity of the σB target expression to σA and competition for RNA polymerase.

A. Steady-state dependence of σB target expression on free σB for different total levels of the housekeeping σ-factor (AT). B. KsigB, the half-maximal constant of the dependence of σB target expression, as a function of the total levels of the housekeeping σ-factor (AT). C. Simplified network diagrams of a positive feedback regulated stress σ-factor σB competing with housekeeping σ-factor σA for RNA polymerase. D-E. Trajectories of free σB (D) and σB target promoter activity (E) in response to stochastic phosphatase input at two different levels of total σA (AT = 9μM-low competition for RNA polymerase; AT = 12μM-high competition for RNA polymerase).

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

(PDF)

S8 Fig. Pulsatile response and RNApol competition in the multiple stress σ-factor model.

A,B. Decoupled σB post-translational and transcriptional components in the simplified model for the competition of stress σ-factors. Cyan and blue curves show the post-translational response at low and high concentration of σB stress signaling protein PB. Black curve shows the transcriptional responses. A step-increase in PB causes a shift in the steady-state post-translational response (from low phosphatase-cyan to high phosphatase-blue) and leads to a pulsatile σB response trajectory (green curve). Concentration of σW stress signaling protein PW was kept fixed at 0.1 μM (A) and 2 μM (B). C,D. Time-course representations of the green trajectories in (A,B) showing σB (C) and σB promoter activity (D) respectively. E. Steady state dependence of σB target promoter activity on the level of free σB for different levels of the σ-factor σW. F. KsigB, the half-maximal constant of the dependence of target expression on σB as a function of the concentration of the stress σ-factor σW for different levels of the housekeeping σ-factor σA. Total RNA polymerase core concentration was kept fixed at 10μM for all simulations.

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

(PDF)

Acknowledgments

The authors are grateful to James Locke for sharing raw data from Ref. [13] and Chet Price for feedback on the manuscript and useful discussion.

Author Contributions

  1. Conceptualization: JN OAI.
  2. Formal analysis: JN AT.
  3. Funding acquisition: OAI.
  4. Methodology: JN AT OAI.
  5. Software: JN AT.
  6. Supervision: OAI.
  7. Writing – original draft: JN AT OAI.

References

  1. 1. Nichols RJ, Sen S, Choo YJ, Beltrao P, Zietek M, et al. (2011) Phenotypic landscape of a bacterial cell. Cell 144: 143–156. pmid:21185072
  2. 2. Nicolas P, Mader U, Dervyn E, Rochat T, Leduc A, et al. (2012) Condition-dependent transcriptome reveals high-level regulatory architecture in Bacillus subtilis. Science 335: 1103–1106. pmid:22383849
  3. 3. Helmann JD (2002) The extracytoplasmic function (ECF) sigma factors. Adv Microb Physiol 46: 47–110. pmid:12073657
  4. 4. Osterberg S, del Peso-Santos T, Shingler V (2011) Regulation of alternative sigma factor use. Annu Rev Microbiol 65: 37–55. pmid:21639785
  5. 5. Hecker M, Pane-Farre J, Volker U (2007) SigB-dependent general stress response in Bacillus subtilis and related gram-positive bacteria. Annu Rev Microbiol 61: 215–236. pmid:18035607
  6. 6. Dufour A, Haldenwang WG (1994) Interactions between a Bacillus subtilis anti-sigma factor (RsbW) and its antagonist (RsbV). J Bacteriol 176: 1813–1820. pmid:8144446
  7. 7. Kang CM, Brody MS, Akbar S, Yang X, Price CW (1996) Homologous pairs of regulatory proteins control activity of Bacillus subtilis transcription factor sigma(b) in response to environmental stress. J Bacteriol 178: 3846–3853. pmid:8682789
  8. 8. Vijay K, Brody MS, Fredlund E, Price CW (2000) A PP2C phosphatase containing a PAS domain is required to convey signals of energy stress to the sigmaB transcription factor of Bacillus subtilis. Mol Microbiol 35: 180–188. pmid:10632888
  9. 9. Voelker U, Voelker A, Haldenwang WG (1996) Reactivation of the Bacillus subtilis anti-sigma B antagonist, RsbV, by stress- or starvation-induced phosphatase activities. J Bacteriol 178: 5456–5463. pmid:8808936
  10. 10. Yang X, Kang CM, Brody MS, Price CW (1996) Opposing pairs of serine protein kinases and phosphatases transmit signals of environmental stress to activate a bacterial transcription factor. Genes Dev 10: 2265–2275. pmid:8824586
  11. 11. Alper S, Dufour A, Garsin DA, Duncan L, Losick R (1996) Role of adenosine nucleotides in the regulation of a stress-response transcription factor in Bacillus subtilis. J Mol Biol 260: 165–177. pmid:8764398
  12. 12. Haldenwang WG (1995) The sigma factors of Bacillus subtilis. Microbiol Rev 59: 1–30. pmid:7708009
  13. 13. Locke JC, Young JW, Fontes M, Hernandez Jimenez MJ, Elowitz MB (2011) Stochastic pulse regulation in bacterial stress response. Science 334: 366–369. pmid:21979936
  14. 14. Young JW, Locke JC, Elowitz MB (2013) Rate of environmental change determines stress response specificity. Proc Natl Acad Sci U S A 110: 4140–4145. pmid:23407164
  15. 15. Igoshin OA, Brody MS, Price CW, Savageau MA (2007) Distinctive topologies of partner-switching signaling networks correlate with their physiological roles. J Mol Biol 369: 1333–1352. pmid:17498739
  16. 16. Grigorova IL, Phleger NJ, Mutalik VK, Gross CA (2006) Insights into transcriptional regulation and sigma competition from an equilibrium model of RNA polymerase binding to DNA. Proc Natl Acad Sci U S A 103: 5332–5337. pmid:16567622
  17. 17. Rollenhagen C, Antelmann H, Kirstein J, Delumeau O, Hecker M, et al. (2003) Binding of sigma(A) and sigma(B) to core RNA polymerase after environmental stress in Bacillus subtilis. J Bacteriol 185: 35–40. pmid:12486038
  18. 18. Delumeau O, Lewis RJ, Yudkin MD (2002) Protein-protein interactions that regulate the energy stress activation of sigma(B) in Bacillus subtilis. J Bacteriol 184: 5583–5589. pmid:12270815
  19. 19. Rami Tzafriri A, Edelman ER (2007) Quasi-steady-state kinetics at enzyme and substrate concentrations in excess of the Michaelis-Menten constant. J Theor Biol 245: 737–748. pmid:17234216
  20. 20. Kalman S, Duncan ML, Thomas SM, Price CW (1990) Similar organization of the sigB and spoIIA operons encoding alternate sigma factors of Bacillus subtilis RNA polymerase. J Bacteriol 172: 5575–5585. pmid:2170324
  21. 21. Voelker U, Dufour A, Haldenwang WG (1995) The Bacillus subtilis rsbU gene product is necessary for RsbX-dependent regulation of sigma B. J Bacteriol 177: 114–122. pmid:8002609
  22. 22. Ma W, Trusina A, El-Samad H, Lim WA, Tang C (2009) Defining network topologies that can achieve biochemical adaptation. Cell 138: 760–773. pmid:19703401
  23. 23. Brody MS, Vijay K, Price CW (2001) Catalytic function of an alpha/beta hydrolase is required for energy stress activation of the sigma(B) transcription factor in Bacillus subtilis. J Bacteriol 183: 6422–6428. pmid:11591687
  24. 24. Cai L, Friedman N, Xie XS (2006) Stochastic protein expression in individual cells at the single molecule level. Nature 440: 358–362. pmid:16541077
  25. 25. Friedman N, Cai L, Xie XS (2006) Linking stochastic dynamics to population distribution: an analytical framework of gene expression. Phys Rev Lett 97: 168302. pmid:17155441
  26. 26. Taniguchi Y, Choi PJ, Li GW, Chen H, Babu M, et al. (2010) Quantifying E. coli proteome and transcriptome with single-molecule sensitivity in single cells. Science 329: 533–538. pmid:20671182
  27. 27. Dufour A, Voelker U, Voelker A, Haldenwang WG (1996) Relative levels and fractionation properties of Bacillus subtilis sigma(B) and its regulators during balanced growth and stress. J Bacteriol 178: 3701–3709 sigma. pmid:8682769
  28. 28. Kang CM, Vijay K, Price CW (1998) Serine kinase activity of a Bacillus subtilis switch protein is required to transduce environmental stress signals but not to activate its target PP2C phosphatase. Mol Microbiol 30: 189–196. pmid:9786195
  29. 29. Marles-Wright J, Grant T, Delumeau O, van Duinen G, Firbank SJ, et al. (2008) Molecular architecture of the "stressosome," a signal integration and transduction hub. Science 322: 92–96. pmid:18832644
  30. 30. Gruber TM, Gross CA (2003) Multiple sigma subunits and the partitioning of bacterial transcription space. Annu Rev Microbiol 57: 441–466. pmid:14527287
  31. 31. Mauri M, Klumpp S (2014) A model for sigma factor competition in bacterial cells. PLoS Comput Biol 10: e1003845. pmid:25299042
  32. 32. Brooks BE, Buchanan SK (2008) Signaling mechanisms for activation of extracytoplasmic function (ECF) sigma factors. Biochim Biophys Acta 1778: 1930–1945. pmid:17673165
  33. 33. Kingston AW, Liao X, Helmann JD (2013) Contributions of the sigma(W), sigma(M) and sigma(X) regulons to the lantibiotic resistome of Bacillus subtilis. Mol Microbiol 90: 502–518. pmid:23980836
  34. 34. Mascher T (2013) Signaling diversity and evolution of extracytoplasmic function (ECF) sigma factors. Curr Opin Microbiol 16: 148–155. pmid:23466210
  35. 35. Yoshimura M, Asai K, Sadaie Y, Yoshikawa H (2004) Interaction of Bacillus subtilis extracytoplasmic function (ECF) sigma factors with the N-terminal regions of their potential anti-sigma factors. Microbiology 150: 591–599. pmid:14993308
  36. 36. Wiegert T, Homuth G, Versteeg S, Schumann W (2001) Alkaline shock induces the Bacillus subtilis sigma(W) regulon. Mol Microbiol 41: 59–71. pmid:11454200
  37. 37. Heinrich J, Hein K, Wiegert T (2009) Two proteolytic modules are involved in regulated intramembrane proteolysis of Bacillus subtilis RsiW. Mol Microbiol 74: 1412–1426. pmid:19889088
  38. 38. Huang X, Gaballa A, Cao M, Helmann JD (1999) Identification of target promoters for the Bacillus subtilis extracytoplasmic function sigma factor, sigma W. Mol Microbiol 31: 361–371. pmid:9987136
  39. 39. Ray JC, Igoshin OA (2010) Adaptable functionality of transcriptional feedback in bacterial two-component systems. PLoS Comput Biol 6: e1000676. pmid:20168997
  40. 40. Narula J, Kuchina A, Lee DY, Fujita M, Suel GM, et al. (2015) Chromosomal Arrangement of Phosphorelay Genes Couples Sporulation and DNA Replication. Cell 162: 328–337. pmid:26165942
  41. 41. Miyashiro T, Goulian M (2008) High stimulus unmasks positive feedback in an autoregulated bacterial signaling circuit. Proc Natl Acad Sci U S A 105: 17457–17462. pmid:18987315
  42. 42. Tiwari A, Balazsi G, Gennaro ML, Igoshin OA (2010) The interplay of multiple feedback loops with post-translational kinetics results in bistability of mycobacterial stress response. Phys Biol 7: 036005. pmid:20733247
  43. 43. Tiwari A, Ray JC, Narula J, Igoshin OA (2011) Bistable responses in bacterial genetic networks: designs and dynamical consequences. Math Biosci 231: 76–89. pmid:21385588
  44. 44. Chong S, Chen C, Ge H, Xie XS (2014) Mechanism of transcriptional bursting in bacteria. Cell 158: 314–326. pmid:25036631
  45. 45. Corrigan AM, Tunnacliffe E, Cannon D, Chubb JR (2016) A continuum model of transcriptional bursting. Elife 5.
  46. 46. Levine JH, Lin Y, Elowitz MB (2013) Functional roles of pulsing in genetic circuits. Science 342: 1193–1200. pmid:24311681
  47. 47. Cookson NA, Mather WH, Danino T, Mondragon-Palomino O, Williams RJ, et al. (2011) Queueing up for enzymatic processing: correlated signaling through coupled degradation. Mol Syst Biol 7: 561. pmid:22186735
  48. 48. Nielsen AA, Voigt CA (2014) Multi-input CRISPR/Cas genetic circuits that interface host regulatory networks. Mol Syst Biol 10: 763. pmid:25422271
  49. 49. Segall-Shapiro TH, Meyer AJ, Ellington AD, Sontag ED, Voigt CA (2014) A 'resource allocator' for transcription based on a highly fragmented T7 RNA polymerase. Mol Syst Biol 10: 742. pmid:25080493
  50. 50. Dhooge A, Govaerts W, Kuznetsov YA (2003) MATCONT: A MATLAB package for numerical bifurcation analysis of ODEs. Acm Transactions on Mathematical Software 29: 141–164.
  51. 51. Barndorff-Nielsen OE, Shephard N (2001) Non-Gaussian Ornstein-Uhlenbeck-based models and some of their uses in financial economics. Journal of the Royal Statistical Society Series B-Statistical Methodology 63: 167–207.
  52. 52. Arrieta-Ortiz ML, Hafemeister C, Bate AR, Chu T, Greenfield A, et al. (2015) An experimentally supported model of the Bacillus subtilis global transcriptional regulatory network. Mol Syst Biol 11: 839. pmid:26577401
  53. 53. Nannapaneni P, Hertwig F, Depke M, Hecker M, Mader U, et al. (2012) Defining the structure of the general stress regulon of Bacillus subtilis using targeted microarray analysis and random forest classification. Microbiology 158: 696–707. pmid:22174379