https://orcid.org/0000-0003-0670-1943
Figures
Abstract
Single-cell approaches are revealing a high degree of heterogeneity, or noise, in gene expression in isogenic bacteria. How gene circuits modulate this noise in gene expression to generate robust output dynamics is unclear. Here we use the Bacillus subtilis alternative sigma factor σB as a model system for understanding the role of noise in generating circuit output dynamics. σB controls the general stress response in B. subtilis and is activated by a range of energy and environmental stresses. Recent single-cell studies have revealed that the circuit can generate two distinct outputs, stochastic pulsing and a single pulse response, but the conditions under which each response is generated are under debate. We implement a stochastic mathematical model of the σB circuit to investigate this and find that the system’s core circuit can generate both response types. This is despite one response (stochastic pulsing) being stochastic in nature, and the other (single response pulse) being deterministic. We demonstrate that the main determinant for whichever response is generated is the degree with which the input pathway activates the core circuit, although the noise properties of the input pathway also biases the system towards one or the other type of output. Thus, our work shows how stochastic modelling can reveal the mechanisms behind non-intuitive gene circuit output dynamics.
Author summary
Experimental advances have enabled the measurement of the dynamics of gene regulatory networks at the single-cell level. This has revealed surprising heterogeneity, or noise, in gene expression between genetically identical cells. This noise can be beneficial, for example by allowing a bacterial population to ‘bet-hedge’ against future environmental change by having a few cells randomly enter a stress prepared state. Here, we use mathematical modelling to investigate a noisy gene regulatory circuit, the σB mediated general stress response of the bacterium Bacillus subtilis. By creating a stochastic model of the σB network, we can replicate the two different response behaviours the system has previously been shown to produce. Interestingly, the first of these behaviours (a single response pulse) is non-stochastic in nature, while the other (stochastic pulsing) is distinctly stochastic. We scan system parameters (properties) to determine how these affect which behaviour the system produces. We show that relatively minor perturbations can push the system from a single pulse response to a stochastic pulsing regime, helping explain previous contradictory experimental results. Our work furthers understanding of how noise in gene expression can enable novel gene circuit dynamics.
Citation: Loman TE, Locke JCW (2023) The σB alternative sigma factor circuit modulates noise to generate different types of pulsing dynamics. PLoS Comput Biol 19(8): e1011265. https://doi.org/10.1371/journal.pcbi.1011265
Editor: Tobias Bollenbach, Universitat zu Koln, GERMANY
Received: October 10, 2022; Accepted: June 12, 2023; Published: August 4, 2023
Copyright: © 2023 Loman, Locke. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: Scripts for generating all of the figures in this article, as well as the simulations on which they are based, can be found at https://gitlab.developers.cam.ac.uk/slcu/teamjl/loman_locke_2023.
Funding: The work in the Locke laboratory is supported by a fellowship from the Gatsby Foundation (GAT3272/GLC to JCWL). TEL has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No. 721456. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
1 Introduction
As we gain more detailed knowledge of the composition of biochemical systems, the use of mathematical chemical reaction network (CRN) models is not only becoming increasingly feasible, but also an important tool to discern how these networks process information. CRN models have been used for such wide applications as the study of animal development [1, 2], plant circadian rhythms [3, 4], and bacterial stress response [5, 6]. In addition, they are an important tool in synthetic biology, where they can be used to design circuits with desired regulatory properties [7]. Traditionally such models have been deterministic in nature. However, stochastic modelling techniques are now often used to capture noisy system dynamics [8]. Such noise can both be caused by low molecule numbers of system components (intrinsic noise) and external quasi-random effects such as cell cycle stage and variability in cellular energy supply (extrinsic noise) [9–11]. Noise in gene expression has been shown to be important to the dynamics of many systems, enabling new and sometimes even beneficial biological phenotypes [12–16]. In this article, we use the general stress response alternative sigma factor of Bacillus subtilis, σB, as a model system to study the influence of noise on the output of a cellular network.
Sigma factors are an essential component of bacterial RNA polymerase, allowing it to recognise and bind to genomic promoter regions. Under standard environmental conditions, a so-called housekeeping sigma factor is used. However, in response to environmental cues, an alternative sigma factor can be activated, replacing the housekeeping sigma factor in RNA polymerase and redirecting the bacterium’s transcriptional program [17–21]. σB of B. subtilis is one of the most well-studied alternative sigma factors. It responds to two classes of stress: energy stress (caused by ATP depletion) and environmental stress (caused by factors such as ethanol, heat, or salt). In response to these stresses, it activates approximately 200 genes involved in the general stress response [22–24]. σB activity is controlled by a core circuit, which can be activated by two distinct upstream pathways (which are triggered by energy and environmental stress, respectively) (Fig 1) [22, 24, 25].
Under non-stress conditions, σB is bound, and held inactive by, its anti-sigma factor (RsbW, W in the figure). An anti-anti-sigma factor (RsbV, V in the figure) is inactivated by a phosphate group. The upstream pathways are triggered by two distinct types of stress (environmental and energy stress), each activating their respective phosphatase (RsbTU and RsbP, respectively). These will dephosphorylate, and thus activate, RsbV. Once activated, RsbV binds RsbW, which simultaneously releases σB in a partner switching mechanism. This permits σB to activate the general stress response of B. subtilis. RsbW is a kinase that re-phosphorylates RsbV, allowing RsbW to re-bind σB and shut the activation off. In addition, σB activates the production of itself, RsbW, and RsbV, creating a mixed positive/negative feedback loop. The environmental stress response phosphatase complex, RsbTU, consists of the co-factor RsbT and the phosphatase RsbU (T and U, respectively, in the figure). The availability of RsbT is controlled by a multi-protein complex called the stressosome and the phosphatase RsbX (both not depicted in the figure for simplicity). The energy stress response phosphatase, RsbP (P in the figure), also depends on a second protein, RsbQ (not depicted in this figure).
The core σB circuit consists of σB, its anti-sigma factor (RsbW), and its anti-anti-sigma factor (RsbV). In the absence of stress, σB is bound to, and kept inactive by, RsbW. Meanwhile, RsbV is kept inactivated by a phosphate group. The energy and environmental stress sensing pathways each activate a phosphatase (RsbP and RsbTU, respectively) [26, 27]. These dephosphorylate RsbV, activating it. Active RsbV binds RsbW in a partner switching mechanism that also releases (and thus activates) σB [28–31]. Active σB will not only activate the general stress response, but also activate an operon containing itself, RsbW, and RsbV, creating a mixed positive/negative feedback loop [32, 33]. RsbW also acts as a kinase, rephosphorylating RsbV, thus resetting the system in the absence of stress [30].
Advances in single-cell fluorescent microscopy have enabled the study of dynamic gene expression at the single-cell level [34, 35]. These techniques have been deployed on alternative sigma factor systems, with reporters for alternative sigma factor activity often displaying heterogeneous activation dynamics (even across isogenic populations experiencing homogeneous inputs) [36–42]. It has been suggested that this heterogeneity can be generated by the amplification of intrinsic cellular noise, and that it may create beneficial, variable, phenotypes through bet-hedging behaviours [43]. Investigations of σB using single-cell fluorescent microscopy have shown that the two stress types (energy and environmental stress) can produce two distinct response behaviours. In response to energy stress, σB displays a persistent stochastic pulsing behaviour. It has been observed that σB is activated in short (∼ 1 hour) activity pulses, randomly and independently distributed across the population and time [36, 38]. The environmental stress response is instead a single response pulse, where the population displays a synchronised pulse (∼ 1 hour long) at the time of stress onset, but the system thereafter remains OFF as stress persists [44]. A model of a simplified σB network was implemented in [36]. Later, a model based on the full network CRN was implemented in [45]. Using this model, the authors showed how the relative stoichiometery of the synthesis rates of σB, RsbW, and RsbV can cause the network to function as an ultrasensitive negative feedback loop, generating pulses. The model can explain the two pathways’ different responses by assuming that environmental stress produces a tightly controlled supply of phosphatase (triggering the core circuit on stress onset only), while energy stress instead produces a naturally fluctuating phosphatase supply (repeatedly triggering the core circuit).
Further experiments, however, have demonstrated that the two response behaviours are not unique to their respective pathways. In [37] it was shown that energy stress could produce a single response pulse, without following stochastic pulses. In addition, it was shown that the environmental stress response, in backgrounds where the stressosome (a multi-protein complex that controls the availability of the RsbT component of RsbTU) was mutated, could also produce a stochastic pulsing-type behaviour. The system’s ability to generate both behaviours from both pathways suggests that it is the core σB circuit that is able to produce both a single response pulse and stochastic pulsing. It is possible that minor perturbations in the two stress sensing pathways could bias the core circuit towards either behaviour.
In this article, we develop a stochastic model of the σB network, based on a previous deterministic model of σB regulation [45] (hereafter referred to as the Narula model). By considering noise in all system components we confirm that the core circuit can generate both stochastic pulsing and single pulse dynamics without requiring different assumptions about the noise properties of upstream energy and environmental stress pathways. Next, we find a minimal set of parameters that tune the system’s response behaviour, and investigate how the tuning of these parameters modulates which behaviour is displayed. We show how a single system parameter (the effective dephosphorylation rate of RsbV) is the primary determinant of which behaviour is produced. Furthermore, we show how the system transitions from a single pulse response, to stochastic pulsing, to oscillations as this parameter is increased. This provides an explanation for recent data showing that the energy stress pathway and the environmental stress pathway can display both stochastic pulsing and single pulse response dynamics. Finally, we demonstrate how properties of noise in the system’s upstream pathways may still bias it towards either response behaviour.
2 Results
2.1 The core σB circuit can produce both the stochastic pulsing and the single response pulse behaviour
The Narula model recreates the energy and environmental stress response behaviours (stochastic pulsing and a single response pulse, respectively) by assuming that the two upstream pathways present qualitatively different input processes to the core σB circuit [45]. Since then, it has been shown experimentally that both pathways can generate both responses under specific environmental or genetic perturbations [37]. This suggests that the core circuit can generate both behaviours. The Narula model relies on a stochastic input process (to a deterministic model) to generate noisy behaviours. By instead using a stochastic CRN interpretation of this model, we investigated the core circuit’s ability to generate both response behaviours from constant stress inputs.
To implement a stochastic CRN, we used the chemical Langevin equation (CLE) [46]. The CLE, which is an established technique for modelling noise in gene regulatory networks [47–50], is an approximation of the Gillespie algorithm [51, 52]. The Gillespie algorithm permits accurate simulations of the actual reaction events of a CRN, correctly taking into account the inherent randomness of the reactions to create stochastic simulations. Gillespie simulations, however, have long simulation runtimes. By approximating the system as a stochastic differential equation (SDE), the CLE permits fast stochastic simulations. While the exactness of this approximation is reduced with increasing noise levels, it is more accurate than other approximations such as the linear noise approximation [53]. For our application, the faster simulation times the CLE permits were important for performing the parameter sweeps carried out later in the paper. The CLE also allows independent tuning of the system’s noise amplitude through the introduction of a system size parameter Ω. Using an equivalent formulation, our noise term in the CLE contains a noise scaling parameter (where η = 0 means no noise and η = 1 means unscaled noise), and we use this to investigate the noise amplitude’s effect on the system. To enforce non-negativity in our CLE simulations, we used absolute values of any negative numbers in the noise terms [54]. To confirm that our conclusions are not due to the assumptions of the CLE, throughout the paper we verify that key results can be reproduced using the Gillespie algorithm [51, 52]. We note that the CLE and Gillespie algorithm both model intrinsic noise only, and not extrinsic noise. While sources of extrinsic noise (such as cell cycle variations) could be added to the model [55], we omitted these as current experimental evidence points to transcriptional noise driving alternative sigma factor pulse initiation and not extrinsic factors such as the cell cycle [36, 38].
We first simulated our CLE adaptation of the Narula model under a range of conditions to test what output dynamics are generated. By subjecting the model, at various noise levels, to various degrees of stress, we noted that modulation of η (the degree of noise) tunes the amplitude of the transient response pulse (Fig 2A and 2B). However, while such modulations created increased fluctuations in the asymptotic state, they were unable to produce stochastic pulsing (S1 Fig). It is known that pulses are related to oscillations [56] (a phenomenon that can also be further emphasised by noise [57, 58]). By screening bifurcation diagrams of all parameters we identified kK2 (the re-phosphorylation rate of the anti-anti-sigma factor RsbV) as the only parameter that can induce oscillations under non-absurd tunings (S2 and S3 Figs). Tuning of a combination of parameters can also induce pulsing. However, we wished to minimise the number of model parameter values we modified. By modifying a smaller set of parameters we reduce the dimensionality of the parameter space we need to explore, reducing computational complexity and enabling more complete scans. We thus selected kK2 only as our proxy for the system’s proneness to pulsing.
(A,B) For various noise amplitudes (η), the CLE adaptation of the Narula model is exposed to a stress step (at t = 0, red dashed line). (A) Mean [σB] in response to the input (average over n = 150 simulations). The amplitude of the single response pulse increases with noise. (B) For each value of η, four different simulations are shown. There is little variation between the individual trajectories. (C,D) Illustration of our measures for the degree with which the system exhibits the single response pulse (C) and stochastic pulsing (D) behaviours. A simulation is divided into a transient phase (t ∈ [0.0, 5.0]) and an asymptotic phase (t ∈ [5.0, 200.0], but we note that in this figure the x-axis is cut to t = 50 to better display both phases). Next, we find the maximum activity in the transient phase, the maximum activity in the asymptotic phase, and the mean activity in the asymptotic phase. The single response pulse measure (C) is defined as the maximum transient activity (orange line) divided by the maximum asymptotic activity (magenta line). The stochastic pulsing measure (D) is defined as the maximum asymptotic activity (orange line) divided by the mean asymptotic activity (magenta line). In practice, a mean measure over several (n > 50) simulations is always used. (E,F) The parameters η (noise amplitude) and kK2 (a proxy for the system’s proneness to oscillations) are varied. For each parameter combination, the maximum magnitude of the single response pulse (E) and stochastic pulsing (F) behaviours that can be achieved by varying the parameter pstress (20.0 μM< pstress < 200.0 μM) is found and plotted. (G) The single response pulse behaviour is maximised at (kK2, η) = (15.0hr−1, 0.04). For these values, four simulations are shown. (H) The stochastic pulsing behaviour is maximised at (kK2, η) = (9.0hr−1, 0.06). For these values, four simulations are shown. (G,H) These simulations demonstrate that the CRN of the Narula model can generate both behaviours while exposed to intrinsic noise only. (I,J) It is possible to recreate both response behaviours, single pulse response (I) and stochastic pulsing (J), using the Gillespie algorithm. This demonstrates that the responses are not dependent on the modelling approach used. Parameter values and other details on simulation conditions for this figure are described in S1–S3 Tables.
We wished to determine to what degree the system can generate the single response pulse and stochastic pulsing behaviours. To do this, we developed an automatic measure that takes a parameter set as input, and returns the magnitude of both behaviours (Fig 2C and 2D). An extensive description of the measure can be found in Section 4.5.1, while this paragraph contains a briefer description. First; we utilise that σB pulse durations rarely exceed 1 hour to divide σB stress response simulations (with stress added at time t = 0) into two different time phases: a transient phase (0 − 5 hours after stress addition) and an asymptotic phase (> 5 hours after stress addition). Next, we find the maximum σB activity in the transient and asymptotic phases, as well as the mean σB activity in the asymptotic phase. The degree of the single response pulse behaviour, in a single simulation, is measured as the maximum σB activity in the transient phase divided by the maximum activity in the asymptotic phase. Similarly, for stochastic pulsing, the definition is the maximum activity in the asymptotic state divided by the mean activity in the asymptotic state. In both cases, to gain a more precise measure for a given parameter set, we take each behaviour’s mean magnitude across a large number of simulations (between 50 and 200).
Using these measures, we evaluated the system’s proneness to both behaviours across η − kK2 space (Fig 2E and 2F). This revealed how, in addition to the single pulse response behaviour, as intrinsic noise is added to the ODE model (by increasing η), the stochastic pulsing behaviour emerges (but only when kK2 is small enough). Through sample simulations, we demonstrated that the core circuit can generate both response behaviours, and that this holds even when the phosphatase concentration in the upstream pathways remains constant (Fig 2G and 2H, and S4 Fig). We also confirmed that these behaviours can be recreated using the Gillespie algorithm (Fig 2I and 2J). For the Gillespie algorithm simulations, as we could not tune a noise parameter η, instead we modified a larger set of parameters to generate the different dynamics (S3 Table). Finally, we investigated the initiation of the stochastic pulses, noting that these typically are preceded by a drop in RsbW concentration (S5 Fig).
2.2 The primary determinant of the system’s response behaviour is the total RsbV dephosphorylation rate
We have shown that through modulation of η and kK2, the system can reproduce both stochastic pulsing and single response pulse behaviours. However, η and kK2 are properties of the core circuit, and these are identical for the energy and environmental stress responses. Thus, the system must be able to generate both a single response pulse and stochastic pulsing while keeping η and kK2 fixed and varying the properties of the upstream pathway only. In the Narula model, this pathway is determined by the parameters pstress, kB5, kD5, and kP. These denote the total amount of phosphatase (pstress), and the rates at which it binds (kB5), dissociates from (kD5), and dephosphorylates RsbV (kP), according to these reactions:
Here, kB5, kD5, and kP, are fixed throughout a simulation, while pstress denotes the post-stress amount of phosphatase (with the phosphatase amount changed from pinit to pstress at the time of stress onset).
We first wanted to determine for which set of parameters the core circuit can generate both response behaviours as the upstream parameters alone are varied. To do this, we defined a measure of the system’s ability to robustly generate both response behaviours (stochastic pulsing and single response pulse) while varying only the upstream parameters (S6 Fig). By scanning this measure across η − kK2 space, we found that (η, kK2) = (0.025, 7.0hr−1) optimises this measure (S7 Fig), meaning that this parameter combination allows the core circuit to generate both behaviours, with the behaviour that is generated being dependent on the values of upstream parameters (S8 Fig). We next fixed η and kK2 at this optimum, and then proceeded to investigate how the parameters pstress, kB5, kD5, and kP modulate the system’s response.
To determine which parameters governing the upstream pathways are most important for determining the response behaviours, we next scanned the two behaviours’ magnitudes across pstress-kB5-kD5-kP space (S9 Fig). Our scans show that both behaviours occur in regions along the curve pstress ⋅ kP = C, for various values of C (Fig 3A and 3B). This suggests that the total RsbV dephosphorylation efficiency (the amount of phosphatase, pstress, times its dephosphorylation efficiency, kP) is an important determinant for the system’s response. We next performed parameter substitutions pprod = pstress ⋅ kP (pprod = total RsbV dephosphorylation efficiency) and pfrac = pstress/kP. By redoing the magnitude scans using this substitution, we confirmed that the total dephosphorylation efficiency, pprod, is important for determining the response behaviour (Fig 3C and 3D, and S10–S15 Figs). Finally, we evaluated the behavioural magnitudes’ sensitivity to change in the four parameters pprod, pfrac, kB5, and kD5 (Fig 3E and 3F). This showed that both behaviours are robust as the values of pfrac, kB5, and kD5 change, but are sensitive to the value of pprod. This suggested that the dynamics encoded by the upstream pathway are primarily defined by the value of pprod. Thus, for further analysis, we set kB5 = 3600μM−1hr−1, kD5 = 18hr−1 (their original values) and pfrac = 100Mhr1 (an intermediate value), and let the upstream pathway be defined by the value of pprod only (greatly reducing the dimensionality of the parameter space we need to consider for further analysis).
(A,B) The magnitude of the single pulse response (A) and stochastic pulsing (B) behaviours across kP-pstress-space. The regions where either behaviour occurs are similar to that of the kP ⋅ pstress = C curve (show for various values, green dashed lines). (C,D) Parameter substitutions generate a parameter pprod = pstress ⋅ kP, to which both behaviours are sensitive, and a parameter pfrac = pstress/kP, to which both behaviours are insensitive. While the regions corresponding to either behaviour are adjacent, they do not overlap. The stochastic pulsing behaviour exists for slightly larger values of pprod, as compared to the single response pulse behaviour. (E,F) For 4096 different parameter sets, we characterise both behaviours’ sensitivity to change (, E, and
, F) in the parameters pprod, pfrac, kB5, and kD5 (Section 4.5.4). We do this by evaluating
and
for the four different parameters across all 4096 parameter sets. We then put each set of 4096 evaluations in ascending order and plot them in E and F. For a few parameter sets, the behaviours show some sensitivity to kB5. However, pprod has the far largest effect on either behaviour. In both cases, changes to pfrac and kD5 have little effect on the system. Hence, these lines coincide (both following the x-axis closely). Parameter values and other details on simulation conditions for this figure are described in S2 and S5 Tables.
Next, we investigated how the system’s response is modulated by the upstream pathway’s critical parameter pprod (the total RsbV dephosphorylation efficiency). We simulated the system for pprod values ranging from small (where the response is absent) to large (where the response saturates) (as pprod = pstress⋅kP, this increase corresponds to an increase in stress input) (Fig 4). For small values of pprod the system exhibits no response and σB concentrations are low as stress is added. As pprod is increased, the system exhibits the single response pulse behaviour. Here, the amplitude of the pulse increases with the stress (Fig 4D and 4E). This is in agreement with previous experiments claiming the σB environmental stress response to be amplitude modulated [44]. For larger values of pprod, the system displays stochastic pulsing (Fig 4F and 4G). The model predicts that the pulse frequency increases with the stress, again similar to experimental data for the energy stress response [36]. Next, the system transitions into a region of oscillation (Fig 4H–4J). For strong stresses, the system exhibits a single pulse, followed by an elevated level of σB activity. The activity in this asymptotic state increases with the magnitude of the stress, however, this increase saturates for large enough stresses (Fig 4K–4M). These last two behaviours have not yet been observed experimentally. A bifurcation-stability analysis of the system revealed that oscillatory behaviours appear for lower values of pprod in stochastic simulations compared to deterministic ones (S16 Fig). This further demonstrates how noise can influence σB dynamics. Finally, we demonstrate that a similar transition between single pulse, stochastic pulsing, and oscillatory dynamics can be recreated using the Gillespie algorithm (S17 Fig). Since Gillespie simulations cannot be performed using the substituted parameter pprod, we instead used pstress (however, since pprod scales linearly with pstress their transitions are equivalent).
(A) The magnitude of the two behaviours for our selected parameter set (kK2, η, pfrac, kB5, kD5) = (7.0 h-1, 0.025, 100.0 μM h-1, 3600.0 μM-1h-1, 18.0 h-1). 12 different selected values of pprod (used in B-M) are marked with grey lines. (B-M) For 12 different values of pprod a single simulation is displayed (stress added at t = 0, red dashed line). (B) For pprod small, the system does not respond. (C,D) As pprod is increased, the system exhibits a single response pulse. The amplitude increase with pprod. (E,F) For larger values of pprod, stochastic pulsing is exhibited. The frequency of the pulses increases with pprod. (G-I) As the stress is increased further, the system enters a limit cycle. (J-L) For large pprod, the system exhibits a single response pulse, and then enters a persistent state of elevated σB activity. The activity in this state increases with pprod. (M) For large enough stresses, the system saturates at some maximum activity. An expanded version of this figure, including bifurcation analysis of system steady state properties, can be found in S16 Fig. A similar transition, but generated through Gillespie algorithm simulations, can be found in S17 Fig. Parameter values and other details on simulation conditions for this figure are described in S4 and S5 Tables.
2.3 Properties of noise in the upstream pathway can bias the network towards either response behaviour
Next, we investigated how noise in the upstream pathway affects the system’s response. We noted that relatively minor perturbations in the stress magnitude can push the system between the two behaviours (Fig 4E and 4F). However, experiments have shown that the two response types are robust under large variations in stress levels. It is possible that the structure of the upstream pathways makes their respective response behaviours more stable. Previously we have assumed an absence of noise in the availability of the phosphatase. By letting the phosphatase switch between an active and an inactive state, we added intrinsic upstream noise through the CLE. This introduced the parameters ηamp and ηfreq, denoting the amplitude and the frequency of the upstream noise, respectively (S18 Fig and Section 4.2). We next proceeded to investigate how the values of ηamp and ηfreq affect the system’s ability to robustly generate the two behaviours as pprod is varied.
We developed a measure for the system’s ability to generate one behaviour distinctly under perturbations to the upstream pathway (with these perturbations previously simplified to modifying the parameter pprod) (Section 4.5.3). The measure tests, as pprod is varied, if the system generates one behaviour uniquely while being unable to produce the other one. We scanned this distinctness-of-behaviour measure (for both behaviours) across ηamp-ηfreq space (Fig 5A and 5B). Our scan suggested that the single response pulse behaviour is distinct for small values of ηamp. Meanwhile, stochastic pulsing requires low ηfreq and high ηamp. This suggested that while upstream noise is not required to produce either response behaviour, the presence (or lack) thereof, makes either behaviour more robust (Fig 5C–5F).
(A,B) The amplitude (ηamp) and frequency (ηfreq) of the upstream pathway’s noise is varied. Plots show, for each parameter combination, the distinctness of the single response pulse (A) and stochastic pulsing (B) behaviours. The distinctness (of either behaviour) designates the system’s ability to uniquely generate that behaviour (and not the other behaviour) as a parameter is varied (here pprod). This measure is described in detail in Section 4.5.3 (there designated Dsrp or Dsp). Light green dots mark parameter sets optimising either behaviour’s distinctness. (C,D) For the parameter sets that maximise the distinctness of the single response pulse (C) and stochastic pulsing (D) behaviours, the magnitudes of the two behaviours are shown as functions of pprod. For each parameter set, 7 selected values of pprod are marked with grey lines. (E,F) For the 7 parameter sets marked in C (E), and D (F), a single simulation is displayed (stress added at t = 0, red dashed line). As pprod is varied, the two behaviours are generated much more robustly than what they were for the parameter set in Fig 4. Parameter values and other details on simulation conditions for this figure are described in S6 and S7 Tables.
Finally, our separation of the upstream noise from that of the core circuit allowed us to investigate the relative effect of the two types of noise on the system’s behaviour. We measured the relative distinctness of the two behaviours (Section 4.5.3) across ηamp-η space, where η denotes the amplitude of the core circuit’s noise (to reduce combinatorial complexity, we set ηfreq = 1). Here, the two stress sensing pathways will have the same core noise amplitude, but different upstream noise amplitudes (ηamp). We note that a network where the core noise is low, and the noise amplitude is low in one pathway and high in the other, should be able to robustly generate the two types of behaviours (Fig 6).
(A,B) For each combination of core circuit noise amplitude (η) and upstream pathway noise amplitude (ηamp) we vary the parameter pprod (over the interval 20 < pprod < 200) and calculate each behaviour’s distinctness (Section 4.5.3). The total amount of noise in the system, rather than how it is distributed between the core and upstream pathway, is an important determinant for both behaviours’ occurrence. (A) The single pulse response behaviour is distinct when both types of noise are low. (B) The stochastic pulsing behaviours require some amount of noise (in either pathway), but do diminish if the total amount of noise becomes too large. It is especially prominent when core circuit noise amplitude (η) is small and upstream pathway noise amplitude (ηamp) is intermediately valued. Parameter values and other details on simulation conditions for this figure are described in S7 Table.
3 Discussion
In recent years, the importance of cellular noise to generate biological phenotypes has been demonstrated in a range of systems. These advances have been aided by the development of techniques for dynamic single-cell measurements of cellular components. Such techniques have revealed that the σB network responds to two types of inputs (energy and environmental stress) through two distinct response behaviours (a single response pulse and stochastic pulsing, respectively). Here, we use the CLE to implement noise in the Narula model (which was previously based on the deterministic reaction rate equation). This has allowed us to study how noise modulates the network’s two responses. By modelling the intrinsic noise of the system’s biochemical reactions, we demonstrated how the core circuit can reproduce both response behaviours, without making assumptions about the distinguishing properties of the two upstream pathways. This is in agreement with recent experiments that suggest that both stress types can generate both responses, as experimental conditions are varied [36–38, 44].
To determine how the two upstream pathways generate distinct responses while feeding into a common core circuit, we characterised the parameters defining the pathways. We showed that the response was mainly determined by the rate of RsbV dephosphorylation (pprod). The two responses can thus be explained by the energy stress sensing pathway dephosphorylating RsbV at a higher rate, as compared to the environmental one. Furthermore, we showed how the system transitioned through a range of response behaviours, including single pulse and stochastic pulsing dynamics, as this crucial property was varied.
We also investigated how noise in the upstream pathways affected which behaviour was favoured. We found that the single response pulse behaviour was favoured by low amplitude noise, while the stochastic pulsing behaviour required larger (but not too large) amplitude noise. Furthermore, we noted that the stochastic pulsing behaviour was deterred by high-frequency noise (when fluctuations in the upstream pathway grow larger than the time scale on which the core circuit operates, the core is likely to only see a deterministic mean input). By modelling the amplitude of noise in the core and the pathway separately, we showed how a system with low noise in the core could optimally generate the two response behaviours from two upstream pathways. Here, the environmental stress sensing pathway would be distinctly less noisy than the energy stress sensing pathway. Our unbiased approach thus helps justify assumptions in previous models, where the energy stress pathway was simulated with a gamma distributed Ornstein-Uhlenbeck process whilst the environmental stress was modelled without noise [36, 45], but goes further as it allows us to quantify the effects of different noise levels on each behaviour. Our results suggest that the components of the environmental stress sensing pathway may regulate cellular noise. Indeed, this pathway contains the stressosome protein complex, which controls the availability of phosphatase, and mutations of which permit environmental stress to also generate a stochastic pulsing type response [37, 42]. This suggests that the stressosome complex may reduce cellular noise.
Previous experimental work disagrees on the type of dynamics generated by the energy and environmental stress pathways [36–38, 44]. By showing that the two behaviours are highly related, we can help explain these differing results. We have shown how relatively minor changes in core circuit activation strength, or pathway noise levels, can cause a switch between the two behaviours. This suggests that differences in the stressors used, or possibly the experimental set-up, could affect which response was generated in each set of experiments.
Finally, our results further demonstrate how noise may have non-intuitive effects on dynamic systems [56–58]. This includes our result that the amplitude of the initial response pulse strictly increases with the noise amplitude (Fig 2A and 2B). It also includes our bifurcation-stability analysis of the system (S16 Fig), which demonstrates that oscillatory behaviours appear for lower values of pprod in the stochastic system (as compared to the deterministic one).
To make them tractable for analysis, models must simplify the real system. In this work, we used the chemical Langevin approach (CLE) to model system noise. This is an approximation of the Gillespie algorithm that, due to its simulation speed, allowed us to simulate over a large range of parameters and noise values. However, it is satisfying that we could confirm, using Gillespie algorithm simulations, the key results of our work. This is that the core circuit can generate both single pulse and stochastic pulsing dynamics, and that the system transitions from a single pulse, to stochastic pulsing, to oscillations as the rate of dephosphorylation of RbsV is increased. The fact that simulations using both the CLE and the Gillespie algorithm give qualitatively similar dynamics increases the robustness of our conclusions.
In the future, it will be interesting to examine further, both experimentally and in models, the source of the noise in the σBcircuit. Although technically challenging due to the small size of the proteins involved and difficulties with fluorescent tags affecting function, it would be powerful to examine the levels of RsbW, RsbV and σBsimultaneously in single cells. This would allow us to use the experimentally measured noise level of each component to directly inform our models. On the modelling side, a more detailed model could include noise due to transcriptional and translational bursts [59, 60], whereas here transcription and translation are combined for simplicity. In addition, we assume uncorrelated noise for each reaction channel, but there is interesting theoretical work that suggests that the coupled expression of genes in operons could lead to correlations in noise that can affect the system dynamics [61]. Further simplifications include that dilution terms are modelled as independent, while in reality they are correlated across all species, as well as the fact that we choose to model intrinsic noise only. Although focusing on intrinsic noise is motivated by our experimental understanding of the system, it would be interesting to examine how extrinsic noise might also affect output dynamics. For example, noise due to the partitioning of molecules at division events is not considered, meaning we could be underestimating the noise due to division [62]. It will be important to address how such additional complexity can affect the σB circuit in future models. More generally, stochastic pulsing dynamics have been observed in multiple alternative sigma circuits [38], as well as other stress pathways in microbes [63], and stochastic modelling could allow us to further address whether there are any general principles behind these dynamics. Finally, in Narula (2016) it was shown that the dynamics of the σB circuit are robust to competition from other sigma factors [45]. It would be interesting to investigate the effect of intrinsic noise on such competition dynamics.
As methods for stochastic CRN modelling have become more widespread [64, 65], the effect of noise on cellular systems has increasingly been characterised [66–68]. This has not only been important to understand the design principles of native circuits, but also when designing synthetic ones. Indeed, being able to quantify the effect of intrinsic noise is required to optimally design synthetic regulator networks. This will be especially true in the context of sigma factors, as their ability to create orthogonal regulatory systems has gathered interest as synthetic regulators [69–71]. In this article, we have shown that intrinsic cellular noise can generate novel stochastic behaviours in a genetic circuit, and even affect the properties of what is otherwise a purely deterministic behaviour. As the noise properties of more systems are characterised (both experimentally and through models), it will increasingly become clear whether stochastic dynamics, and not deterministic ones, constitute the rule rather than the exception.
4 Methods
4.1 Model implementation
As a starting point, we use the CRN model implementation of the core σB circuit model as introduced in Narula (2016) [45]. Our only modification is that the phosphatase concentration before and after stress have been formalised as model parameters. The model describes the concentration of σB, RsbW, and RsbV, as well as the various complexes and dimers they form. It also models the input phosphatase, P, as a species of the system.
The strength of the stress is determined by the concentration of the species P. In the original model, environmental stress is modelled as a step, where [P] is increased from a low to a higher level at the time of stress (generating a single response pulse). Energy stress is modelled by letting the phosphate concentration adopt a (pre-simulated) gamma-distributed Ornstein-Uhlenbeck process, and then performing reaction rate equation based ODE simulations with this noisy input (creating a stochastic pulsing output). Since we will use the CLE to model system noise, we do not need to add specific noise to the concentration of P. Hence, we can allow all inputs to be step increases in [P]. Formalising this, we let the parameter pinit be the concentration of P before stress, and pstress be the concentration of P after stress.
4.1.1 Reactions.
The model consists of 27 reactions, which are provided in Table 1.
4.1.2 Parameters.
The model consists of 24 parameters, all of which stay constant throughout a simulation. The parameter values, as they were initially set in [45], are given in Table 2. Throughout the article, the parameters are modified as we investigate how this changes the system’s behaviour (exact details of which parameter values are used for each set of simulations can be found in S1–S7 Tables). Finally, we have also introduced an additional 3 parameters to the model: pinit is the pre-stress concentration of the phosphatase, pstress is the concentration of the phosphatase in the stressed state, and η is the amplitue of noise in the system (using linear scaling of the noise in the CLE).
4.1.3 SDE and ODE systems.
We will primarily use the CLE to simulate the model. From a set of chemical reactions, it describes a system of SDEs according to
(1)
where xi denotes the concentration of the ith species, vi,j the net stoichiometric change in the ith species as result of the jth reaction, and
is the propensity of the jth reaction. dWj denotes a Wiener process with 0 mean and unit variance. Here, the noise scaling parameter is marked as η.
Due to the additional noise terms, this system becomes too large to yield any useful information when displayed in SDE form. We will instead only write out the ODE system (as generated by the RRE), noting that the SDE system can be unambiguously generated by the CRN’s reactions (Table 1) according to the previous equation.
4.2 Modifications to the Narula model
In the CLE adaptation of the Narula model, we can investigate the σB response to a step increase in phosphatase input. Here, the upstream pathway activating either RsbTU (environmental stress) or RsbP (energy stress) produces a constant, non-fluctuating, level of phosphatase in response to the stress (while the amount of phosphatase is constant, we note that the relative fraction between phosphatase that is free, or bound to RsbV, fluctuate). In this section, we expand the CLE adaptation of the Narula model, adding fluctuations to the upstream phosphatase input. This allows us to investigate how upstream noise affects the system’s response behaviour.
Upstream noise is implemented by assuming that the phosphatase has an active and an inactive state, and that it switches between these according to the law of mass action. This allows us to model the upstream noise through the CLE. By introducing a second CLE noise-scaling parameter for the phosphatase, we can scale the upstream noise separately from that of the core circuit. Through the way we write the new parameters, we can introduce one, ηamp, which scales the amplitude of the noise in the phosphatase, and one, ηfreq, which scales the frequency of the noise in the phosphatase (S18 Fig). Our modified model introduces 1 new species, 3 new reactions, and 2 new parameters.
Similarly to the original Narula model, stress is modelled by increasing the total amount of phosphatase in the system. At the time of stress addition, [P] and [PI] are both increased by pstress − pinit (setting the total amount of phosphatase in the system to 2pstress, with on average half being in the active state).
4.2.1 Reactions.
The modified model introduces one new species, PI (inactive phosphatase), and three reactions. These are the activation and the deactivation of the phosphatase. In addition, phosphatase bound to RsbV can also deactivate (thus dissociating from RsbV) (Table 3).
4.2.2 Parameters.
The modified model introduces 2 new parameters (Table 4). The new parameter ηamp scales the noise amplitude in the reactions in Table 3. Meanwhile, the old parameter η still scales the noise of all the reactions in Table 1. The new parameter ηfreq scales the rate of phosphatase activation and deactivation. This has no effect on the system’s deterministic properties, but will scale the frequency of the fluctuations in the levels of active phosphatase. The parameters pinit and pstress, while in principle similar to in the previous model, have their properties slightly altered. Instead of representing (before and after stress, respectively) the total amount of phosphatase ([PI] + [P] + [P-VP]), they represent the average level of active phosphatase (before and after stress, respectively). The total amount of phosphatase is 2pinit (before stress) and 2pstress (after stress).
4.2.3 ODE System.
While we will primarily use the SDE system as generated by the CLE, due to the addition of the noise terms, the system becomes too large to yield any useful information when displayed in SDE form. We will instead only write out the ODE system (as generated by the RRE). The SDE system can be unambiguously generated by the CRN’s reactions (Tables 1 and 3) through the CLE.
4.3 Simulations
The models were implemented in the Julia programming language using the Catalyst.jl modelling package [72]. Simulations were carried out using the DifferentialEquations.jl package [73].
For Gillespie simulations, we use the Direct method (Gillespie’s direct method), which is the recommended method for small problems [51, 52, 74]. These simulations also require a stepping algorithm (which is used internally to manage the simulation). Since no reactions were time-dependent, we used the recommended SSAStepper. All Gillespie simulations were performed using the reactions in Table 1, and the parameters described in S3 Table.
For the SDE simulations, we used the implicit Euler-Maruyama method, using fixed time steps [75]. Due to numeric error, differential equations occasionally produce negative species concentration. To prevent these from crashing simulations (due to negative numbers occurring in the CLE’s square root noise term), every term within a square root was replaced with its absolute value (e.g. would become
, where dW20 is a Wiener process of the CLE SDE). This has no effect as long as species concentrations stay positive. However, if numeric errors and the noise term push a species into a negative concentration, the absolute value prevents the solver from causing an error due to attempting to take the square root of a negative number [54]. For full details of the simulation, as well as the opportunity to reproduce them, please see the provided scripts (see Section 4.6 for a link to the code).
4.4 Bifurcation analysis
Bifurcation analysis was carried out using the BifurcationKit package (which tracks steady states using pseudo-arclength continuation) [76].
4.5 Automated behaviour measures
To aid our analysis, we develop four different measures of the system’s properties. Of these, the first three are related (measuring the system’s ability to generate either behaviour under various circumstances). The last one is a sensitivity measure of whenever either behaviour is sensitive to change in a specified parameter.
4.5.1 Measure of behavioural magnitude.
Throughout the paper, we wish to determine the magnitude with which either behaviour occurs at a given parameter set . To do so, we simulate the system n times, where, n = 50, 100, or 200 (values of n vary from figure to figure and are described in S2, S5 and S7 Tables). Each simulation last 200 hours after the addition of stress (and begins 10 hours before stress). For each simulation we measure:
- The amplitude of the transient pulse. This is measured as the maximum activity of the simulation in the transient phase (t ∈ [0.0, 5.0]).
- The maximum asymptotic pulse amplitude. This is measured as the maximum activity of the simulation in the asymptotic phase (t ∈ [5.0, 200.0]).
- The mean asymptotic activity. This is measured as the mean activity of the simulation in the asymptotic phase (t ∈ [5.0, 200.0]).
Using these, we define the magnitude of the single response pulse behaviour () and the stochastic pulsing behaviour (
) as:
- Single response pulse: The amplitude of the transient pulse divided by the maximum asymptotic pulse amplitude.
- Stochastic pulsing: The maximum asymptotic pulse amplitude divided by the mean asymptotic activity.
By, for the single response pulse, dividing by the maximum asymptotic pulse amplitude we penalise the single response pulse magnitude for the occurrence of stochastic pulses. This ensures that for a single parameter set, the measures cannot both score highly. These measures are illustrated in Fig 2C and 2D.
The five-hour window for the transient phase was chosen to ensure the first (asymptotic) pulse was fully captured. We note that it is possible that this interval might include further pulses, but we still capture this possible sustained pulsing behaviour using our asymptotic phase. To ensure that 5 hours captured the initial pulse, we made an automatic scan of a large number (10, 000) of random parameter sets of our model. None of these exhibited an initial pulse longer than 5 hours.
Finally, to provide additional intuition for these measures, we describe how they will work in two different examples. First, consider a simulation with an oscillatory output. Here, the mean asymptotic activity will typically be half of the maximum asymptotic activity, yielding a stochastic pulsing measure of 2. This value is low compared to a simulation that exhibits stochastic pulsing, where the stochastic pulsing measure often reaches 50. Next, consider a simulation where a single stochastic pulse is exhibited in the asymptotic phase (t ∈ [5.0, 200.0]). This can maximise the magnitude of the stochastic pulsing behaviour, as further pulses will increase the mean asymptotic activity, but only if the further pulses are of lower amplitude, as otherwise the maximum would also increase. We note though that this is only true for the individual simulation. In practice, our measures are always evaluated as the mean over a large (n > 50) number of simulations. If one simulation-realisation of a parameter set exhibits a single asymptotic pulse, it is likely that other simulations will exhibit none. For these, the stochastic pulsing measure will equal 1, reducing the overall measure for the parameter set (taken as a mean over all the simulations).
4.5.2 Measure of the system’s ability to robustly generate both behaviours.
Next, we wish to measure to what extent the system can generate both behaviours, as a single parameter is tuned. This is used in S7 Fig to find fixed values for the core parameters so that both behaviours can be generated by modulating the upstream parameters only. The behaviours should be distinct, that is, the system should clearly exhibit one behaviour to a higher magnitude than the other (and be able to do so for both behaviours). With all other parameters fixed, define Msrp(p) and Msp(p) as the magnitude of the two behaviours (both functions of a single parameter p). We define
(2)
with
defined similarly. These are the degrees to which either behaviour surpasses the other. We then define our measure (of the degree with which the system can, as a parameter p is varied, distinctly generate both behaviours):
(3)
where (pstart, pend) is the interval over which we sample the parameter p to calculate the integrals (In practice sums over a dense grid of parameter values are used to estimate the integrals). This measure is illustrated in S6 Fig. Finally, we note that 0 < Dsrp,sp(p) < 1.
4.5.3 Measure of the system’s ability to distinctly generate one behaviour.
Next, we wish to measure to what extent the system can generate one behaviour distinctly. That is, to what extent, for a given parameter region, only one behaviour is generated. By optimising this measure, one gets a system with a preference for one of the two behaviours. We do this by modifying the previously defined Dsrp,sp(p) measure. Again, we will define our measure as a function of only a single parameter, p (in practice it will be carried out using the parameter pprod). Our new measures Dsrp(p) and Dsp(p) can be defined as
(4)
with Dsp(p) defined similarly. Here, (pstart, pend) is the interval over which we sample the parameter p. Using the terminology in S6 Fig, we have
. By squaring the numerator, we create a preference for regions in parameter space where the behaviour is prominent (not only prominent in comparison to the other behaviour).
4.5.4 Measure of the behaviours’ sensitivity to the phosphatase parameters.
Finally, we wish to measure to what extent the magnitude of a behaviour changes as we tune a single parameter, p. This will be used in Fig 3E and 3F to determine which upstream parameters are important for determining whichever behaviour is generated. Our parameters are sampled at discrete values on a grid, . We define our measures
and
(where p is the values of the model’s remaining parameters, which are kept fixed) as:
(5)
with
defined similarly. The factor (n − 1) is added to remove any bias introduced by the density with which the grid is sampled (fewer grid samples mean longer distances between the individual values, these are emphasised by the square, this factor compensates for this). By squaring each value we give preference for sudden dramatic changes, as opposed to gradual change with the parameter.
For simplicity, the and
measures shown in Fig 3E and 3F has been normalised by the maximum value of the y-axis.
4.6 Code availability
Scripts for generating all of the figures in this article, as well as the simulations on which they are based, can be found at https://gitlab.developers.cam.ac.uk/slcu/teamjl/loman_locke_2023. All scripts are written in the Julia programming language. This enables the definition of a Project.toml file, defining the exact package versions used. This will enable anyone to replicate the exact conditions under which all figures were generated.
Supporting information
S1 Fig. Stochastic pulsing is not achievable in the model by tuning η alone.
For each combination of stress magnitude (pstress) and noise amplitude (η) four simulations are shown (stress added at red dashed line, t = 0). Parameter values and other details on simulation conditions for this figure are described in S1 Table.
https://doi.org/10.1371/journal.pcbi.1011265.s001
(PDF)
S2 Fig. Only through the tuning of two parameters can oscillations be achieved.
(A) Bifurcation diagrams for the various parameters, each plot shows three diagrams (for the varying stress levels pstress = 0.05 μM, 0.20 μM, and 0.80 μM). The stars mark the parameter value for the original Narula model. Each x-axis is log10 scaled, and if the parameter’s original value is p0, it is varied over the range (p0/10, 10p0), corresponding to a tenfold decrease and increase in the target parameter, respectively. Only by tuning kK2 or λW can instability be achieved. For some parameters (kK2, kP, and F), the curve for pstress = 0.80 μM reaches much larger values compared to the other two curves, making these hard to distinguish. To avoid figure crowding, periodic orbits are not displayed in these diagrams, however, they are instead shown in S3 Fig. (B) Bifurcation diagram for the parameter pstress (the magnitude of the stress) over the interval (0.1 μM,10.0 μM), with the x-axis log10 scaled. Instability cannot be produced by tuning pstress only. Parameter values and other details on simulation conditions for this figure are described in S1 Table.
https://doi.org/10.1371/journal.pcbi.1011265.s002
(PDF)
S3 Fig. Bifurcation diagrams with periodic orbits shown.
Four of the bifurcation diagrams in S2 Fig display periodic orbits, these are shown in more details here. The stars mark the parameters’ values for the original Narula model and all x-axes are log10 scaled. (A) Bifurcation diagram with respect to the parameter kK2. (B-D) Bifurcation diagram with respect to the parameter λV, for three different values of pstress. Note that for these the y-axes are log10 scaled (to help show the periodic orbits more clearly), unlikely in S2 Fig where they are linearly scaled. Parameter values and other details on simulation conditions for this figure are described in S1 Table.
https://doi.org/10.1371/journal.pcbi.1011265.s003
(PDF)
S4 Fig. Both behaviours can be reproduced even as we remove all noise from upstream reactions.
Here we test whether removing all noise from phosphatase reactions can still allow both response dynamics. In Fig 2 we demonstrate that both system behaviours can be recreated even as phosphatase levels remain constant (unlike in [45] which assumed noisy concentrations of P). To do this, we use the CLE. The CLE adds noise to reaction channels, rather than species concentrations (with the former affecting the latter). This means that there’s still noise in P due to noise in the reactions involving P. (A,B) Here we remove noise from all reactions involving P (and components containing P). To do this, we use our modified Narula model that allows us to tune upstream noise (Section 4.2) and set ηamp = 0.0. Next, we recreate both the single response pulse (A, pprod = 50.0, η = 0.01) and stochastic pulsing (B, pprod = 25.0, η = 0.09) behaviours. This demonstrates that both behaviours can be generated by the core circuit, and are not dependent on any form of upstream noise. Stress is added at red dashed lines (t = 0) and each plot shows four simulations. Full parameter sets for this figure are described in S6 Table.
https://doi.org/10.1371/journal.pcbi.1011265.s004
(PDF)
S5 Fig. The stochastic pulses in σB are preceded by a reduction in total RsbW concentration.
(A) In the time leading up to a stochastic pulse, the total amount of σB (blue, scale marked at the left) and RsbW (blue, scale marked at the right) is plotted. Four simulations are shown. (B) For the same four simulations, the system is shown in σB-RsbW phase space. The state at the beginning of the simulation is marked with a blue dot. Just before the pulse is initiated (σB levels start to rise), the total RsbW concentration dips. Although harder to see, the dip can also be distinguished in (A). Parameter values and other details on simulation conditions for this figure are described in S1 Table.
https://doi.org/10.1371/journal.pcbi.1011265.s005
(PDF)
S6 Fig. A measure of the system’s ability to generate both behaviours distinctly as a single parameter is tuned.
We find the two functions Msrp(p) and Msp(p) of our target parameter (this example used pstress) (Section 4.5.2). We define three areas: Asrp is the area which is beneath the Msrp(p) curve but above Msp(p), with Asp defined similarly. We also define Asrp,sp as the area which is beneath both curves. Finally, our measure is defined as . Parameter values and other details on simulation conditions for this figure are described in S2 Table.
https://doi.org/10.1371/journal.pcbi.1011265.s006
(PDF)
S7 Fig. Heatmap showing which combinations of kK2 and η enables the system to optimally generate both behaviours.
The function measures the system’s ability to distinctly generate both behaviours, while varying the parameter pstress only (Section 4.5.2). The optimal value, (kK2, η) = (7.0hr−1, 0.025), is found close to the bottom left corner (light green dot). We chose the upper limit of η (0.15) as beyond this point behaviours start to become obscured by noise levels. Parameter values and other details on simulation conditions for this figure are described in S2 Table.
https://doi.org/10.1371/journal.pcbi.1011265.s007
(PDF)
S8 Fig. An optimised parameter set is able to generate both behaviours by varying pstress only.
(A,B) The system’s response to stress (red dashed line, t = 0), for (kK2, η) = (7.0hr−1, 0.025) and pstress = 0.24 μM (A) and pstress = 0.28 μM (B). Each plot shows four simulations. Parameter values and other details on simulation conditions for this figure are described in S1 Table.
https://doi.org/10.1371/journal.pcbi.1011265.s008
(PDF)
S9 Fig. The shape of the two behaviours’ regions of occurrence is stable as kB5 and kD5 are changed.
Heatmaps showing the two behaviours occurrences in kB5-kD5-space. In all plots, the behaviours region of occurrence is similar to a kP ⋅ pstress = C curve. Parameter values and other details on simulation conditions for this figure are described in S2 Table.
https://doi.org/10.1371/journal.pcbi.1011265.s009
(PDF)
S10 Fig. Heatmaps describing the magnitude of the single response pulse behaviour for various values of kB5 and kD5.
Each heatmap describes the behaviour’s magnitude as the parameters pprod (x-axis) and pfrac (y-axis) are varied. A total of 36 heatmaps are plotted and placed in a 6x6 grid for a range of values of kB5 and kD5. There is a distinct spike in magnitude as pprod is varied. Changes to pfrac, pprod, and kD5 all have some effect on the magnitude, but not as distinct as changes to pprod. Parameter values and other details on simulation conditions for this figure are described in S5 Table.
https://doi.org/10.1371/journal.pcbi.1011265.s010
(PDF)
S11 Fig. Heatmaps describing the magnitude of the stochastic pulsing response behaviour for various values of kB5 and kD5.
Each heatmap describes the behaviour’s magnitude as the parameters pprod (x-axis) and pfrac (y-axis) are varied. A total of 36 heatmaps are plotted and placed in a 6x6 grid for a range of values of kB5 and kD5. There is a distinct spike in magnitude as pprod is varied. Changes to pfrac, pprod, and kD5 all have some effect on the magnitude, but not as distinct as changes to pprod. Parameter values and other details on simulation conditions for this figure are described in S5 Table.
https://doi.org/10.1371/journal.pcbi.1011265.s011
(PDF)
S12 Fig. Heatmaps describing the magnitude of the single response pulse behaviour for various values of pprod and pfrac.
Each heatmap describes the behaviour’s magnitude as the parameters kD5 (x-axis) and kB5 (y-axis) are varied. A total of 36 heatmaps are plotted and placed in a 6x6 grid for a range of values of pprod and pfrac. Only for one value of pprod do changes in the other parameters have a major effect on the behaviour’s magnitude. Parameter values and other details on simulation conditions for this figure are described in S5 Table.
https://doi.org/10.1371/journal.pcbi.1011265.s012
(PDF)
S13 Fig. Heatmaps describing the magnitude of the stochastic pulsing response behaviour for various values of pprod and pfrac.
Each heatmap describes the behaviour’s magnitude as the parameters kD5 (x-axis) and kB5 (y-axis) are varied. A total of 36 heatmaps are plotted and placed in a 6x6 grid for a range of values of pprod and pfrac. Only for one value of pprod do changes in the other parameters have a major effect on the behaviour’s magnitude. Parameter values and other details on simulation conditions for this figure are described in S5 Table.
https://doi.org/10.1371/journal.pcbi.1011265.s013
(PDF)
S14 Fig. Heatmaps describing the magnitude of the single response pulse behaviour for various values of pfrac and kD5.
Each heatmap describes the behaviour’s magnitude as the parameters pprod (x-axis) and kB5 (y-axis) are varied. A total of 36 heatmaps are plotted and placed in a 6x6 grid for a range of values of pfrac and kD5. There is a distinct spike in magnitude as pprod is varied. Changes to pfrac, pprod, and kD5 have little effect on the behaviour’s magnitude. Parameter values and other details on simulation conditions for this figure are described in S5 Table.
https://doi.org/10.1371/journal.pcbi.1011265.s014
(PDF)
S15 Fig. Heatmaps describing the magnitude of the stochastic pulsing response behaviour for various values of pfrac and kD5.
Each heatmap describes the behaviour’s magnitude as the parameters pprod (x-axis) and kB5 (y-axis) are varied. A total of 36 heatmaps are plotted and placed in a 6x6 grid for a range of values of pfrac and kD5. There is a distinct spike in magnitude as pprod is varied. Changes to pfrac, pprod, and kD5 have little effect on the behaviour’s magnitude. Parameter values and other details on simulation conditions for this figure are described in S5 Table.
https://doi.org/10.1371/journal.pcbi.1011265.s015
(PDF)
S16 Fig. Fig 4 with a system bifurcation diagram added.
Subplot A is identical to subplot A of Fig 4. Subplots C-N are identical to subplots B-M of Fig 4. Please see Fig 4 for legends of these subplots. (B) Bifurcation diagram, showing the system’s steady state as a function of pprod. After a region of inactivity, the steady state becomes unstable (implying a limit cycle). As pprod increases further, the steady state becomes stable, with an increasing concentration that eventually saturates. The transition in the bifurcation diagram corresponds to the transition in C-N (stability at an inactive state, a limit cycle, stability at an active state). We note that the bifurcation diagram is computed from the deterministic (ODE) system, while subfigures A, and C-N all use the stochastic (SDE) system. Adding noise to a non-linear system affects the parameter values at which stability occurs. This explains that the region of instability in (B) occurs for different pprod values as compared to A. While the deterministic bifurcation analysis in B cannot be directly translated to our stochastic system, it is still worth noting this similarity between the two cases. Parameter values and other details on simulation conditions for this figure are described in S5 Table.
https://doi.org/10.1371/journal.pcbi.1011265.s016
(PDF)
S17 Fig. The behavioural transition, as pprod is varied, in Fig 4 can be recreated using the Gillespie algorithm.
In Fig 4 we observed that the system undergoes a behavioural transition as pprod is increased from small to high. For small values of pprod the system is inactive. As pprod is increased the system exhibits, in order, single response pulse, stochastic pulsing, oscillating, and persistent activity, behaviours. Here, we recreate the same transition using Gillespie simulations. While we in Fig 4 vary the parameter pprod, we never introduce this parameter substitution for the Gillespie approach. We instead vary the parameter pstress. However, since pprod = pstress ⋅ kP, the two transitions should be equivalent. (A-L) Gillespie simulation of the Narula model for various values of pstress. Each frame contains 4 simulations, and stress is added at the red dashed lines (t = 0). The behavioural transition from Fig 4 is recreated. Parameter values and other details on simulation conditions for this figure are described in S3 Table.
https://doi.org/10.1371/journal.pcbi.1011265.s017
(PDF)
S18 Fig. In the modified Narula model, two parameters (ηamp and ηfreq), allow us to scale the amplitude and frequency of the noise in the input phosphatase.
By splitting the phosphatase into an active state (P) and an inactive state (Pi), with only the active state being able to form the P-VP complex, we can introduce fluctuations in phosphatase levels. Before the addition of stress, the total amount of phosphatase ([Pi] + [P] + [P-VP]) is 2pinit (with, on average, at every time point, half being in the active state). At the time of stress onset (here at the red dashed line at t = 0), the amount of phosphatase is increased to 2pstress (again, with on average half being active). The frequency of the fluctuations scales with the parameter ηfreq. The amplitude of the fluctuations scales with the parameter ηamp. Parameter values and other details on simulation conditions for this figure are described in S7 Table.
https://doi.org/10.1371/journal.pcbi.1011265.s018
(PDF)
S1 Table. Parameter values for simulation of the CLE adaptation of the Narula model.
For each figure where the model is simulated, the parameter values used for the simulations are marked. If not marked, the following parameter values are used: kBw = 3600 μM-1hr-1, kDw = 18 hr-1, kB1 = 3600 μM-1hr-1, kB2 = 3600 μM-1hr-1, kB3 = 3600 μM-1hr-1, kB4 = 1800 μM-1hr-1, kB5 = 3600 μM-1hr-1, kD1 = 18 hr-1, kD2 = 18 hr-1, kD3 = 18 hr-1, kD4 = 1800 μM-1hr-1, kD5 = 18 hr-1, kK1 = 36 hr-1, kK2 = 36 hr-1, kP = 180 hr-1 hr-1, kDeg = 0.7 hr-1, v0 = 0.4 μM-1hr-1, F = 30, K = 0.2 μM, λW = 4, λV = 4.5, η = 0.05, pinit = 0.001 μM, pstress = 0.4 μM. Finally, in certain figures, some parameter values are varied as marked on the figures. Which parameters are varied across each is marked in the last column.
https://doi.org/10.1371/journal.pcbi.1011265.s019
(PDF)
S2 Table. Parameter values for parameter scans of the CLE adaptation of the Narula model.
For each figure where we perform parameter scans of the model’s behaviour, the parameter values used are marked. If not marked, the following parameter values are used: kBw = 3600 μM-1hr-1, kDw = 18 hr-1, kB1 = 3600 μM-1hr-1, kB2 = 3600 μM-1hr-1, kB3 = 3600 μM-1hr-1, kB4 = 1800 μM-1hr-1, kB5 = 3600 μM-1hr-1, kD1 = 18 hr-1, kD2 = 18 hr-1, kD3 = 18 hr-1, kD4 = 1800 μM-1hr-1, kD5 = 18 hr-1, kK1 = 36 hr-1, kK2 = 36 hr-1, kP = 180 hr-1, kDeg = 0.7 hr-1, v0 = 0.4 μM-1hr-1, F = 30, K = 0.2 μM, λW = 4, λV = 4.5, η = 0.05, and pinit = 0.001 μM. The third column denotes how many simulations (n) are performed for each parameter combination. Finally, in certain figures, some parameter values are varied as marked on the figures. Which parameters are varied across each is marked in the last column.
https://doi.org/10.1371/journal.pcbi.1011265.s020
(PDF)
S3 Table. Parameter values for Gillespie simulation of the Narula model.
Here, Ms stand for Molecules. For each figure where the Narula model is simulated using the Gillespie algorithm, the parameter values used for the simulations are marked. If not marked, the following parameter values are used: kBw = 3600 Molecules-1hr-1, kDw = 18 hr-1, kB1 = 3600 Molecules-1hr-1, kB2 = 3600 Molecules-1hr-1, kB3 = 3600 Molecules-1hr-1, kB4 = 1800 Molecules-1hr-1, kB5 = 3600 Molecules-1hr-1, kD1 = 18 Molecules-1, kD2 = 18 hr-1, kD3 = 18 hr-1, kD4 = 1800 Molecules-1hr-1, kD5 = 18 hr-1, kK1 = 36 hr-1, kP = 180 hr-1, λW = 4, λV = 4.5, and pinit = 0 Molecules. Finally, in S18 Fig, the value of pstress is varied as marked in the figure. This is designated by, in this table, putting pstress in the “Varied parameters” column.
https://doi.org/10.1371/journal.pcbi.1011265.s021
(PDF)
S4 Table. Parameter values for simulation of the CLE adaptation of the Narula model with parameter substitution.
For each figure where the model is simulated, the parameter values used for the simulations are marked. If not marked, the following parameter values are used: kBw = 3600 μM-1hr-1, kDw = 18 hr-1, kB1 = 3600 μM-1hr-1, kB2 = 3600 μM-1hr-1, kB3 = 3600 μM-1hr-1, kB4 = 1800 μM-1hr-1, kB5 = 3600 μM-1hr-1, kD1 = 18 hr-1, kD2 = 18 hr-1, kD3 = 18 hr-1, kD4 = 1800 μM-1hr-1, kD5 = 18 hr-1, kK1 = 36 hr-1, kK2 = 36 hr-1, kDeg = 0.7 hr-1, v0 = 0.4 μM-1hr-1, F = 30, K = 0.2 μM, λW = 4, λV = 4.5, η = 0.05, pinit = 0.001 μM, and pfrac = 100 μM hr1. Finally, in certain figures, some parameter values are varied as marked on the figures. Which parameters are varied across each is marked in the last column.
https://doi.org/10.1371/journal.pcbi.1011265.s022
(PDF)
S5 Table. Parameter values for parameter scans of the CLE adaptation of the Narula model with parameter substitution.
For each figure where we perform parameter scans of the model’s behaviour, the parameter values used are marked. If not marked, the following parameter values are used: kBw = 3600 μM-1hr-1, kDw = 18 hr-1, kB1 = 3600 μM-1hr-1, kB2 = 3600 μM-1hr-1, kB3 = 3600 μM-1hr-1, kB4 = 1800 μM-1hr-1, kB5 = 3600 μM-1hr-1, kD1 = 18 hr-1, kD2 = 18 hr-1, kD3 = 18 hr-1, kD4 = 1800 μM-1hr-1, kD5 = 18 hr-1, kK1 = 36 hr-1, kK2 = 36 hr-1, kDeg = 0.7 hr-1, v0 = 0.4 μM-1hr-1, F = 30, K = 0.2 μM, λW = 4, λV = 4.5, η = 0.05, pinit = 0.001 μM, and pfrac = 100 μM hr1. The third column denotes how many simulations (n) are performed for each parameter combination. These simulations are then used to determine behavioural magnitudes (larger n yields smoother plots). Finally, in certain figures, some parameter values are varied as marked on the figures. Which parameters are varied across each is marked in the last column.
https://doi.org/10.1371/journal.pcbi.1011265.s023
(PDF)
S6 Table. Parameter values for simulation of the modified Narula model.
For each figure where the modified Narula model is simulated, the parameter values used for the simulations are marked. If not marked, the following parameter values are used: kBw = 3600 μM-1hr-1, kDw = 18 hr-1, kB1 = 3600 μM-1hr-1, kB2 = 3600 μM-1hr-1, kB3 = 3600 μM-1hr-1, kB4 = 1800 μM-1hr-1, kB5 = 3600 μM-1hr-1, kD1 = 18 hr-1, kD2 = 18 hr-1, kD3 = 18 hr-1, kD4 = 1800 μM-1hr-1, kD5 = 18 hr-1, kK1 = 36 hr-1, kK2 = 36 hr-1, kDeg = 0.7 hr-1, v0 = 0.4 μM-1hr-1, F = 30, K = 0.2 μM, λW = 4, λV = 4.5, η = 0.05, pinit = 0.001 μM, pfrac = 100 μM hr1, ηamp = 0.05, and ηfreq = 1. Finally, in certain figures, some parameter values are varied as marked on the figures. Which parameters are varied across each is marked in the last column.
https://doi.org/10.1371/journal.pcbi.1011265.s024
(PDF)
S7 Table. Parameter values for parameter scans of the modified Narula model.
For each figure where we perform parameter scans of the modified Narula model’s behaviour, the parameter values used are marked. If not marked, the following parameter values are used: kBw = 3600 μM-1hr-1, kDw = 18 hr-1, kB1 = 3600 μM-1hr-1, kB2 = 3600 μM-1hr-1, kB3 = 3600 μM-1hr-1, kB4 = 1800 μM-1hr-1, kB5 = 3600 μM-1hr-1, kD1 = 18 hr-1, kD2 = 18 hr-1, kD3 = 18 hr-1, kD4 = 1800 μM-1hr-1, kD5 = 18 hr-1, kK1 = 36 hr-1, kK2 = 36 hr-1, kDeg = 0.7 hr-1, v0 = 0.4 μM-1hr-1, F = 30, K = 0.2 μM, λW = 4, λV = 4.5, η = 0.05, pinit = 0.001 μM, pfrac = 100 μM hr1, ηamp = 0.05, and ηfreq = 1. The third column denotes how many simulations (n) are performed for each parameter combination. These simulations are then used to determine behavioural magnitudes (larger n yields smoother plots). Finally, in certain figures, some parameter values are varied as marked on the figures. Which parameters are varied across each is marked in the last column.
https://doi.org/10.1371/journal.pcbi.1011265.s025
(PDF)
Acknowledgments
We are grateful for the high performance computing resources provided by the Cambridge Service for Data Driven Discovery. We thank Dr. Katie Abley and Dr Christian Schwall (University of Cambridge) for critical reading of the manuscript.
References
- 1. Eldar A, Dorfman R, Weiss D, Ashe H, Shilo DZ, Barkal N. Robustness of the BMP morphogen gradient in Drosophila embryonic patterning. Nature. 2002;419(6904):304–308. pmid:12239569
- 2. Exelby K, Herrera-Delgado E, Perez LG, Perez-Carrasco R, Sagner A, Metzis V, et al. Precision of tissue patterning is controlled by dynamical properties of gene regulatory networks. Development (Cambridge, England). 2021;148(4):1–14. pmid:33547135
- 3. Zeilinger MN, Farré EM, Taylor SR, Kay SA, Doyle FJ. A novel computational model of the circadian clock in Arabidopsis that incorporates PRR7 and PRR9. Molecular Systems Biology. 2006;2. pmid:17102803
- 4. Foo M, Bates DG, Akman OE. A simplified modelling framework facilitates more complex representations of plant circadian clocks. PLOS Computational Biology. 2020;16(3):1–34. pmid:32176683
- 5. Narula J, Kuchina A, Zhang F, Fujita M, Süel G, Igoshin O. Slowdown of growth controls cellular differentiation. Molecular Systems Biology. 2016;12:871–871. pmid:27216630
- 6. Zorzan I, Del Favero S, Giaretta A, Manganelli R, Di Camillo B, Schenato L. Mathematical modelling of SigE regulatory network reveals new insights into bistability of mycobacterial stress response. BMC bioinformatics. 2021;22(1):558. pmid:34798803
- 7. Chandran D, Copeland WB, Sleight SC, Sauro HM. Mathematical modeling and synthetic biology. Drug Discovery Today: Disease Models. 2008;5(4):299–309. pmid:27840651
- 8. Hou Z, Xin H. Internal noise stochastic resonance in a circadian clock system. Journal of Chemical Physics. 2003;119(22):11508–11512.
- 9. Mcadams HH, Arkin A. Stochastic mechanisms in gene expression. Proceedings of the National Academy of Sciences of the United States of America. 1997;94(3):814–819. pmid:9023339
- 10. Elowitz MB, Levine AJ, Siggia ED, Swain PS. Stochastic gene expression in a single cell. Science. 2002;297(5584):1183–1186. pmid:12183631
- 11. Perez-Carrasco R, Beentjes C, Grima R. Effects of cell cycle variability on lineage and population measurements of messenger RNA abundance. Journal of the Royal Society Interface. 2020;(17). pmid:32634365
- 12. Ackermann M, Stecher B, Freed NE, Songhet P, Hardt WD, Doebeli M. Self-destructive cooperation mediated by phenotypic noise. Nature. 2008;454(7207):987–990. pmid:18719588
- 13. Cerulus B, New AM, Pougach K, Verstrepen KJ. Noise and Epigenetic Inheritance of Single-Cell Division Times Influence Population Fitness. Current Biology. 2016;26(9):1138–1147. pmid:27068419
- 14. Ietswaart R, Rosa S, Wu Z, Dean C, Martin H. Cell-Size-Dependent Transcription of FLC and Its Antisense Long Non-coding RNA COOLAIR Explain Cell-to-Cell Expression Variation. Cell Systems. 2017; p. 622–635.e9. pmid:28624615
- 15. Benzinger D, Khammash M. Pulsatile inputs achieve tunable attenuation of gene expression variability and graded multi-gene regulation. Nature Communications. 2018;9(1). pmid:30166548
- 16. Johnston IG, Bassel GW. Identification of a bet-hedging network motif generating noise in hormone concentrations and germination propensity in Arabidopsis. Journal of the Royal Society Interface. 2018;15(141):20180042. pmid:29643226
- 17. Wösten MMSM. Eubacterial sigma-factors. FEMS Microbiology Reviews. 1998;22(3):127–150. pmid:9818380
- 18. Gruber TM, Gross CA. Multiple Sigma Subunits and the Partitioning of Bacterial Transcription Space. Annual Review of Microbiology. 2003;57:441–466. pmid:14527287
- 19. Feklístov A, Sharon BD, Darst SA, Gross CA. Bacterial sigma factors: A historical, structural, and genomic perspective. Annual Review of Microbiology. 2014;68:357–376. pmid:25002089
- 20. Davis MC, Kesthely CA, Franklin EA, MacLellan SR. The essential activities of the bacterial sigma factor. Canadian Journal of Microbiology. 2017;63(2):89–99. pmid:28117604
- 21. Helmann JD. Where to begin? Sigma factors and the selectivity of transcription initiation in bacteria. Molecular Microbiology. 2019;112(2):335–347. pmid:31119812
- 22. Hecker M, Pané-Farré J, Völker U. SigB-dependent general stress response in Bacillus subtilis and related gram-positive bacteria. Annual Review of Microbiology. 2007;61:215–236. pmid:18035607
- 23. Nannapaneni P, Hertwig F, Depke M, Hecker M, Mäder U, Völker U, et al. Defining the structure of the general stress regulon of Bacillus subtilis using targeted microarray analysis and random forest classification. Microbiology. 2012;158(3):696–707. pmid:22174379
- 24. Rodriguez Ayala F, Bartolini M, Grau R. The Stress-Responsive Alternative Sigma Factor SigB of Bacillus subtilis and Its Relatives: An Old Friend With New Functions. Frontiers in Microbiology. 2020;11:1761. pmid:33042030
- 25. van der Steen JB, Ávila-Pérez M, Knippert D, Vreugdenhil A, van Alphen P, Hellingwerf KJ. Differentiation of function among the RsbR paralogs in the general stress response of bacillus subtilis with regard to light perception. Journal of Bacteriology. 2012;194(7):1708–1716. pmid:22287516
- 26. Voelker U, Voelker A, Haldenwang WG. Reactivation of the Bacillus subtilis anti-σB antagonist, RsbV, by stress- or starvation-induced phosphatase activities. Journal of Bacteriology. 1996;178(18):5456–5463. pmid:8808936
- 27. Vijay K, Brody MS, Fredlund E, Price CW. A PP2C phosphatase containing a PAS domain is required to convey signals of energy stress to the sigmaB transcription factor of Bacillus subtilis. Molecular Microbiology. 2000;35(1):180–188. pmid:10632888
- 28. Benson AK, Haldenwang WG. Characterization of a regulatory network that controls sigma B expression in Bacillus subtilis. Journal of Bacteriology. 1992;174(3):749–757. pmid:1732211
- 29. Benson AK, Haldenwang WG. Bacillus subtilis sigma B is regulated by a binding protein (RsbW) that blocks its association with core RNA polymerase. Proceedings of the National Academy of Sciences of the United States of America. 1993;90:2330–2334. pmid:8460143
- 30. Dufour A, Haldenwang WG. Interactions between a Bacillus subtilis anti-sigma factor (RsbW) and its antagonist (RsbV). Journal of Bacteriology. 1994;176(7):1813–1820. pmid:8144446
- 31. Alper S, Dufour A, Garsin DA, Duncan L, Losick R. Role of adenosine nucleotides in the regulation of a stress-response transcription factor in Bacillus subtilis. Journal of Molecular Biology. 1996;260(2):165–177. pmid:8764398
- 32. Kalman S, Duncan ML, Thomas SM, Price CW. Similar organization of the sigB and spoIIA operons encoding alternate sigma factors of Bacillus subtilis RNA polymerase. Journal of Bacteriology. 1990;172(10):5575–5585. pmid:2170324
- 33. Wise AA, Price CW. Four additional genes in the sigB operon of Bacillus subtilis that control activity of the general stress factor sigma B in response to environmental signals. Journal of Bacteriology. 1995;177(1):123–133. pmid:8002610
- 34. Muzzey D, van Oudenaarden A. Qunatitative Time-Lapse Flourescence Microscopy in Single Cells. Annual Review of Cell and Developmental Biology. 2011;25:301–327.
- 35. Young JW, Locke JCW, Altinok A, Rosenfeld N, Bacarian T, Swain PS, et al. Measuring single-cell gene expression dynamics in bacteria using fluorescence time-lapse microscopy. Nature Protocols. 2011;7(1):80–88. pmid:22179594
- 36. Locke JCW, Young JW, Fontes M, Hernández MJ, Elowitz MB. Stochastic Pulse Regulation in Bacterial Stress Response. Science. 2011;334(6054):366–369. pmid:21979936
- 37. Cabeen MT, Russell JR, Paulsson J, Losick R. Use of a microfluidic platform to uncover basic features of energy and environmental stress responses in individual cells of Bacillus subtilis. PLoS Genetics. 2017;13(7):e1006901. pmid:28727759
- 38. Park J, Dies M, Lin Y, Hormoz S, Smith-Unna SE, Quinodoz S, et al. Molecular Time Sharing through Dynamic Pulsing in Single Cells. Cell Systems. 2018;6(2):216–229. pmid:29454936
- 39. Patange O, Schwall C, Jones M, Villava C, Griffith DA, Phillips A, et al. Escherichia coli can survive stress by noisy growth modulation. Nature Communications. 2018;9(1):5333. pmid:30559445
- 40. Nadezhdin E, Murphy N, Dalchau N, Phillips A, Locke JCW. Stochastic pulsing of gene expression enables the generation of spatial patterns in Bacillus subtilis biofilms. Nature Communications. 2020;11(1):950. pmid:32075967
- 41. Schwall CP, Loman TE, Martins BMC, Cortijo S, Villava C, Kusmartsev V, et al. Tunable phenotypic variability through an autoregulatory alternative sigma factor circuit. Molecular Systems Biology. 2021;17(7):1–16. pmid:34286912
- 42. Hamm CW, Butler DR, Cabeen MT. Bacillus subtilis Stressosome Sensor Protein Sequences Govern the Ability To Distinguish among Environmental Stressors and Elicit Different σB Response Profiles. mBio. 2022;13(6):e0200122. pmid:36409125
- 43. Veening JW, Smits WK, Kuipers O. Bistability, Epigenetics, and Bet-Hedging in Bacteria. Annual review of microbiology. 2008;62:193–210. pmid:18537474
- 44. Young JW, Locke JCW, Elowitz MB. Rate of environmental change determines stress response specificity. Proceedings of the National Academy of Sciences of the United States of America. 2013;110(10):4140–4145. pmid:23407164
- 45. Narula J, Tiwari A, Igoshin OA. Role of Autoregulation and Relative Synthesis of Operon Partners in Alternative Sigma Factor Networks. PLoS Computational Biology. 2016;12(12):e1005267. pmid:27977677
- 46. Gillespie DT. The chemical Langevin equation. Journal of Chemical Physics. 2000;113(1):297–306.
- 47. Mavroudis PD, Scheff JD, Calvano SE, Lowry SF, Androulakis IP. Entrainment of peripheral clock genes by cortisol. Physiol Genomics. 2012;44(11):607–621. pmid:22510707
- 48. Perez-Carrasco R, Guerrero P, Briscoe J, Page KM. Intrinsic Noise Profoundly Alters the Dynamics and Steady State of Morphogen-Controlled Bistable Genetic Switches. PLoS Computational Biology. 2016;12(10):1–23. pmid:27768683
- 49. Manning CS, Biga V, Boyd J, Kursawe J, Ymisson B, Spiller DG, et al. Quantitative single-cell live imaging links HES5 dynamics with cell-state and fate in murine neurogenesis. Nature Communications. 2019;10(1):2835. pmid:31249377
- 50. Perez-Carrasco R, Barnes CP, Schaerli Y, Isalan M, Briscoe J, Page KM. Combining a Toggle Switch and a Repressilator within the AC-DC Circuit Generates Distinct Dynamical Behaviors. Cell Systems. 2018;6(4):521–530.e3. pmid:29574056
- 51. Gillespie DT. A General Method for Numerically Simulating the Stochastic Time Evolution of Coupled Chemical Reactions. Journal of Computational Physics. 1976;22:403–434.
- 52. Gillespie DT. Exact stochastic simulation of coupled chemical reactions. Journal of Physical Chemistry. 1977;81(25):2340–2361.
- 53. Wallace EWJ, Gillespie DT, Sanft KR, Petzold LR. Linear noise approximation is valid over limited times for any chemical system that is sufficiently large. IET Systems Biology. 2012;6(4):102–115. pmid:23039691
- 54. Higham DJ. Modeling and simulating chemical reactions. SIAM Review. 2008;50(2):347–368.
- 55. Shahrezaei V, Ollivier JF, Swain PS. Colored extrinsic fluctuations and stochastic gene expression. Molecular systems biology. 2008. pmid:18463620
- 56.
Lindner B, García-Ojalvo J, Neiman A, Schimansky-Geier L. Effects of noise in excitable systems; 2004.
- 57. Rao CV, Wolf DM, Arkin AP. Control, exploitation and tolerance of intracellular noise. Nature. 2002;420(6912):231–237. pmid:12432408
- 58. Rué P, Süel GM, Garcia-Ojalvo J. Optimizing periodicity and polymodality in noise-induced genetic oscillators. Phys Rev E. 2011;83:061904. pmid:21797400
- 59. Golding I, Paulsson J, Zawilski SM, Cox EC. Real-time kinetics of gene activity in individual bacteria. Cell. 2005;123(6):1025–1036. pmid:16360033
- 60. Friedman N, Cai L, Xie XS. Linking stochastic dynamics to population distribution: An analytical framework of gene expression. Physical Review Letters. 2006;97(16):1–4. pmid:17155441
- 61. Ray C, Igoshin O. Interplay of Gene Expression Noise and Ultrasensitive Dynamics Affects Bacterial Operon Organization. PLoS computational biology. 2012;8:e1002672. pmid:22956903
- 62. Beentjes CHL, Perez-Carrasco R, Grima R. Exact solution of stochastic gene expression models with bursting, cell cycle and replication dynamics. Phys Rev E. 2020;101:032403. pmid:32290003
- 63. Süel GM, Garcia-Ojalvo J, Liberman LM, Elowitz MB. An excitable gene regulatory circuit induces transient cellular differentiation. Nature. 2006;440(7083):545–550. pmid:16554821
- 64. Gillespie DT. Stochastic simulation of chemical kinetics. Annual Review of Physical Chemistry. 2007;58:35–55. pmid:17037977
- 65. Warne DJ, Baker RE, Simpson MJ. Simulation and inference algorithms for stochastic biochemical reaction networks: From basic concepts to state-of-the-art. Journal of the Royal Society Interface. 2019;16(151):20180943. pmid:30958205
- 66. Altschuler SJ, Wu LF. Cellular Heterogeneity: Do Differences Make a Difference? Cell. 2010;141(4):559–563. pmid:20478246
- 67. Balázsi G, van Oudenaarden A, Collins J. Cellular Decision Making and Biological Noise: From Microbes to Mammals. Cell. 2011;144(6):910–925. pmid:21414483
- 68. Tsimring LS. Noise in biology. Reports on Progress in Physics. 2014;77(2):026601. pmid:24444693
- 69. Rhodius VA, Segall-Shapiro TH, Sharon BD, Ghodasara A, Orlova E, Tabakh H, et al. Design of orthogonal genetic switches based on a crosstalk map of σs, anti-σs, and promoters. Molecular Systems Biology. 2013;9:702. pmid:24169405
- 70. Bervoets I, Van Brempt M, Van Nerom K, Van Hove B, Maertens J, De Mey M, et al. A sigma factor toolbox for orthogonal gene expression in Escherichia coli. Nucleic Acids Research. 2018;46(4):2133–2144. pmid:29361130
- 71. Pinto D, Vecchione S, Wu H, Mauri M, Mascher T, Fritz G. Engineering orthogonal synthetic timer circuits based on extracytoplasmic function factors. Nucleic Acids Research. 2018;46(14):7450–7464. pmid:29986061
- 72. Loman TE, Ma Y, Ilin V, Gowda S, Korsbo N, Yewale N, et al. Catalyst: Fast Biochemical Modeling with Julia. bioRxiv. 2022.
- 73. Rackauckas C, Nie Q. DifferentialEquations.jl—A Performant and Feature-Rich Ecosystem for Solving Differential Equations in Julia. Journal of Open Research Software. 2017;5(15):15.
- 74.
Isaacson SA, Ilin V, Rackauckas CV. JumpProcesses.jl; 2022. https://github.com/SciML/JumpProcesses.jl/.
- 75. Rackauckas C, Nie Q. Adaptive methods for stochastic differential equations via natural embeddings and rejection sampling with memory. Discrete Continuous Dyn Syst Ser B. 2017;22(7):2731–2761. pmid:29527134
- 76.
Veltz R. BifurcationKit.jl; 2020. https://hal.archives-ouvertes.fr/hal-02902346.