Time-scale separation and stochasticity conspire to impact phenotypic dynamics in the canonical and inverted Bacillus subtilis core genetic regulation circuits

Abstract

Background: In this work, we study two seemingly unrelated aspects of core genetic nonlinear dynamical control of the competence phenotype in Bacillus subtilis, a common Gram-positive bacterium living in the soil.

Methods: We focus on hitherto unchartered aspects of the dynamics by exploring the effect of time-scale separation between transcription and translation and, as well, the effect of intrinsic molecular stochasticity. We consider these aspects of regulatory control as two possible evolutionary handles.

Results: Hence, using theory and computations, we study how the onset of oscillations breaks the excitability-based competence phenotype in two topologically close evolutionary-competing circuits: the canonical “wild-type” regulation circuit selected by Evolution and the corresponding indirect-feedback inverted circuit that failed to be selected by Evolution, as was shown elsewhere, due to dynamical reasons.

Conclusions: Relying on in-silico perturbation of the living state, we show that the canonical core genetic regulation of excitability-based competence is more robust against switching to phenotype-breaking oscillations than the inverted feedback organism. We show how this is due to time-scale separation and stochasticity.

Graphical abstract

Keywords

Bacillus subtilis / competence / gene regulation / deterministic dynamics / stochastic dynamics

INTRODUCTION

Living systems are at times full of randomness thus unpredictable and, at other times, orderly and predictable. These competing viewpoints are far from being incompatible with each other and seem entrenched in the subtlety of not only how, but why life organizes in the way we observe or increasingly infer by relying on theory. Which of these aspects is most important depends entirely on context that relates to the size of fluctuations. At the core of genetic control, there exists a fine interplay of intrinsic randomness with deterministic dynamics that deeply impacts phenotype. It is fair to say that, at the core genetic control level, living systems are neither fully predictable nor fully random, but in fact reflect the impact of both dynamical aspects. The functional role of biochemical fluctuations – noise – in genetic circuits is reviewed in Ref. [¹]. Here, to explore these aspects further, but in an expanded scope that includes dynamical time-scales, we focus on the canonical Evolution selected Bacillus subtilis bacterium control circuit and on an inverted mutant circuit that Evolution did not favor [²]. The core regulatory gene circuits that underlie the behavior of these organisms are shown in Figure 1. While these systems have been theoretically [³] and experimentally studied at length already [⁴–⁶], in the present work, we expand the scope of study to focus on subtle and as yet unexplored aspects of the interplay of biochemical noise – stochasticity – with time-scale separation between DNA transcription (the synthesis of mRNA) and the subsequent translation of mRNA.

RESULTS

We conduct our investigations as much shepherded as bolstered by the concept demonstrated in Ref. [²] that dynamic oscillations are essentially competence phenotype breaking, and that the Evolution proofed core genetic regulation circuit of B. subtilis must therefore operate in the excitatory regime. We thus consider these two dynamically close yet not biologically equivalent regimes (excitability and oscillation), in the light of accessible evolutionary handles offered by stochasticity manifested as biochemical noise sourced in the paucity of some key molecular components in the system, and by time-scale separation between transcription and translation as it relates to effective time delays in the system. Time delays between transcription and translation occur because DNA transcription and mRNA translation do not take place at the same location in the cell. The mRNA must relocate from the DNA transcription site where it is synthesized to a site where it is translated to protein. This delay is particularly acute for eukaryotes where transcription occurs inside the nucleus while translation occurs outside the nucleus, in the cytoplasm. For prokaryotes such as bacteria, both transcription and translation necessarily occur in the cytoplasm, so delays can be expected to be smaller, but not necessarily negligible.

Canonical B. subtilis regulation

Figure 2A and 2B show the phase portraits of the “2D” infinite time-scale separation dynamics of the canonical B. subtilis organism.ComK and ComS are the core regulators of competence, the ability of the organism to accept exogenous DNA to evolve when ComK exceeds a certain threshold [⁷,⁸]. The depicted dynamics is referred to as “2D” dynamics because, by stipulating the 4D sub-manifold of the parent 6D system (ComK, ComS, mRNA_ComK, mRNA_ComS, MecA_ComK, MecA_ComS) to be at rest (i.e., time derivatives of mRNA_ComK, mRNA_ComS, MecA_ComK, MecA_ComS are set to zero) the parent 6D system collapses into the two remaining dimensions for the concentrations of protein ComK and ComS. Here, MecA is a protease that is integral to the core regulatory system because it enables a competitive process for the degradation of ComK and ComS; MecA_ComK and MecA_ComS are bound states of ComK with MecA and ComS with MecA, respectively. The reactions forming this 4D sub-manifold are faster than all others in the circuit, hence upon adiabatic approximation, only protein ComK and ComS appear in the 2D system [²,⁶]. We refer to the time-scale separation (between transcription and translation) in the 2D system as “infinite”.

In Figure 2A–2C each diagram plots the genetic regulators ComK on the x-axis and ComS on the y-axis, and each plot features the two nullclines of the 2D dynamics (the ComK nullcline in red and the ComS nullcline in green). Nullclines are the loci of dynamical points on the phase diagram where the sub-manifold of ComK dynamics and the sub-manifold of ComS dynamics are individually, that is to say, separately, at rest. Therefore, at the intersections of nullclines, it is clear that the dynamics of the entire system is at rest. Here, the nullclines intersect at three fixed points: #1 a stable spiral, #2 a saddle node, and #3 an unstable spiral. Several tracks emanating at different initial conditions of ComK and ComS (black dots) illustrate the basic dynamical behaviors. We study the dynamics as a function of the key parameter of the feedback regulation, k_1. The parameter k₁ parameterizes the negative feedback of ComK onto ComS (see Methods). As can be verified directly in Figure 2, comparing phase diagrams of Figure 2A to that of Figure 2B, changes in the feedback regulation via k₁ slightly affect the shape of the ComS nullcline visibly impacting mostly the location of fixed point #3. More specifically, in Figure 2A at k₁=0.0333, magenta tracks fall into fixed point #1. Brown tracks started on the other side of the separatrix drawn in cyan, undergo a classic excitatory excursion around unstable fixed point #3 and eventually also reach fixed point #1. But in Figure 2B at k₁=0.0405, the excitatory regime is replaced by a oscillatory regime; the brown track started on the right side of the separatrix now falls into a large amplitude limit cycle. In Figure 2C, a blowup region of Figure 2B, the details of the vicinity of the saddle node are shown in the case of the oscillatory track with k₁=.0405. Whereas the excitatory biological regime of the canonical B. subtilis is as in Figure 2A, we notice that, actually, the excitatory regime lies not far at all from the oscillatory regime shown in Figure 2A and 2C. Indeed, regular large-amplitude oscillations suddenly replace excitatory trajectories as k₁ increases from 0.0333 to 0.0405. It is clear that the implied phenotypes are drastically different between excitatory and oscillatory core genetic control. From the point of view of biology, the phenotype of interest is the appearance of core genetic regulation induced sustained yet transitory intervals of high levels of ComK. During high ComK levels, B. subtilis is competent for accepting exogenous DNA [⁷,⁹,¹⁰]. In the excitatory case, as depicted in Figure 2A, random biochemical fluctuations must first occur in order to trigger the phenotype’s randomly occurring and random-duration dynamical excursions into high ComK (brown tracks in Figure 2A). Indeed, it is easy to see that intrinsic randomness is required because if the system is started anywhere on the left side of the separatrix (cyan line), dynamics will only bring it to the lower fixed point #1, without excursions into the high ComK regime. Only upon a fluctuation that brings the starting point across to the right side of the separatrix, is the system then permitted to complete a large ComK excursion around fixed point #3. In the oscillatory case however, as shown in Figure 2B and 2C, high ComK intervals will also occur but at completely predictable and regularly spaced intervals of time. Whereas in the excitatory case, stochasticity is central to the dynamics, in the oscillatory case however, appears to be relegated to a blurring function.

No other tool than a fixed-point continuation analysis of the dynamics as shown in Figure 2D can summarize the behavior better. Here, k₁ is the continuation parameter on the x-axis, and in black and orange are shown the locations and stability (solid: stable, dash: unstable) of the ComK and ComS upper fixed point #3. The values of ComK and ComS are jointly shown on the y-axis. There is a Hopf bifurcation denoted by “H” at k₁~0.06. To the left of the Hopf bifurcation, oscillations are seen to develop. The upper/lower extrema of these oscillations are shown by the purple/blue and green/red lines for ComS and ComK, respectively. The amplitudes of these oscillations grow steadily from the Hopf as k₁ diminishes until oscillations disappear abruptly towards the left of the figure. We refer to this point as the appearance/disappearance of oscillations because as k₁ increases across this point, suddenly maximum amplitude oscillations appear, whereas as k₁ is decreased across this point, the oscillations suddenly disappear. Oscillations occur due to the presence of the saddle, as detailed in Figure 2C. To the left of the oscillation appearance/disappearance point (as in Figure 2A, the system is in excitatory mode. To the right (as in Figure 2B and 2C), the system is in oscillatory mode. As stated above, to the left of this point, stochasticity is paramount to the canonical phenotype. To the right, is secondary.

Looking for oscillations: the effects of time-scale separation and stochasticity

In this work, we are interested to study out how time-scale separation and stochasticity conspire to affect phenotype. We focus on the oscillatory regime because the phenotype is simple to identify, and appearance of oscillations in an otherwise excitatory regime represents competence-phenotype-breaking dynamical behavior and thus has a distinctly biological meaning. Figure 3 show side by side a purely deterministic 2D (infinite time-scale separation) phase portrait and a discrete event stochastic simulation of imperfect (i.e., finite) time-scale separation. In both cases, k₁=0.05, a value that is expected to produce healthy oscillations of large amplitude. See Figure 2C for details. The limit-cycle in Figure 3A is unmistakable. In Figure 3B, no limit cycle behavior can be seen. The stochastic dynamics is devoid of rotation; the system is still excitatory. We reasoned that the primary difference between these two regimes must be the specific amount of time-scale separation in the simulations shown in Figure 3B; specifically that there might be insufficient time-scale separation in the stochastic system.

In order to study the effect of time-scale separation unfettered by the possible influence of fluctuations, we devised a fluctuation-free version of the stochastic system that we refer to as the “6D system”. The 6D system is, in fact, the parent system to the 2D system that was mentioned earlier. This continuous dynamical system is derived from the set of discrete-event reaction equations underlying the stochastic system, but completely removing the fluctuations. In effect, the 6D system is the infinite number of molecule limit (i.e., so-called thermodynamic limit) of the stochastic system; hence 6D is a deterministic system. More details are given in the Methods section and references [²,⁶]. Since our goal is to use adjustable time-scale separation to study its impact on the competence phenotype, we introduce the concept of “adjustable” time-scale separation to the 6D system. This is in the sense that increasingly slowing down the entire protein sub-manifold (ComK, ComS, MecA_ComK, MecA_ComS) in the 6D system will effectively increasingly separate the timescales of transcription and translation in an adjustable manner, thus mimicking a tunable time delay between the two processes. Note that in practice, in order to fairly install adjustable time-scale separation into the 6D model of Bacillus subtilis, in addition to the actual mRNA translation rates into protein, other rates related the MecA protease dynamics as well as protein degradation rates must also be equally slowed down, so that the entire protein sub-manifold is uniformly slowed down with respect to transcription. Hence, with the entire protein sub-manifold slowed down, while mRNA quickly reaches a steady state, its concomitant protein takes an longer time to settle into its own steady-state. The effective resulting tunable delay is the time difference between mRNA reaching its steady state and the later time when concomitant protein reaches its own steady state. Adjustable time-scale separation is thus implemented using slow down factors acting on translation and related protein manifold rates. Although we refer to such factors as translation slowdown factors, they act on the entire protein dynamics.

We aimed to study the impact of time-scale separation on the dynamics and to do so we conducted bifurcation studies versus two translation slowdown factors (SDT_ComK and SDT_ComS) introduced in the dynamics. See Methods for implementation details. Figure 4 show the 6D system’s fixed point locations as we vary the translation slowdown factors. The presence of the Hopf bifurcation testifies that oscillations are clearly induced by slowing down mRNA translation and associated protein manifold dynamics. We picked SDT_ComS = 0.9 and SDT_ComK = 0.14 from these studies because these factors correspond to steady mid-range healthy oscillations well away from critical points.

In Figure 5A we show, in magenta, a 2D track (with infinite time-scale separation) falling into the expected limit cycle, and we also show (in cyan) a 6D track with SDT_ComS =0.9 and SDT_ComK = 0.14 resulting in similar oscillations. Thus the 2D system and the 6D system are in agreement. Note that a 6D track without translation slowdown (i.e., with translation slow down factors SDT_ComS and SDT_ComK set to 1) fails to fall into oscillations, as is demonstrated in Supplementary Figure S1. Figure 5B shows the equivalent stochastic computation with the slow down factors applied. Despite the presence of translation slowdowns in the stochastic simulation, it is abundantly clear that no oscillations are induced yet; time-scale separation that brings the dynamics to oscillations in the deterministic infinite-number of molecule limit represented by the 6D system is clearly insufficient to bring the corresponding stochastic system with finite stochasticity into similar oscillations. Clearly stochasticity level, being the other parameter in the system, must be impacting the phenotype. This is shown by the simulation presented in Figure 5C. Here, as in Figure 5B slowdown of translation is in effect, but in addition, the system’s intrinsic stochasticity level is lowered by roughly a factor of 3 compared to Figure 5B. We observe the clear unmistakable occurrence of oscillations.

Thus both time-scale separation and stochasticity conspire to mask the oscillation phenotype; or, from a biological viewpoint, time-scale separation and stochasticity work together to maintain the system in an excitatory regime exhibiting the competence phenotype.

Inverted B. subtilis regulation

Hereon, we shift focus to the inverted B. subtilis regulation. This circuit was not selected by Evolution for reasons explained in Ref. [²]. However, because it was genetically re-engineered for investigation, the experimental behavior of the inverted B. subtilis circuit is well understood. In Figure 6A and 6B present the phase portrait of the 2D regulation in which the two regulators are ComK and MecA. Nullclines for these two regulators are shown in red and green, respectively. The three fixed points of the dynamics, located at their mutual intersections, are shown by the labels “1”, “2” and “3”. They are stable spiral, saddle point, and unstable spiral, respectively. The separatrix is shown in cyan. In Figure 6A, k_m=0.1474 (excitatory regime, no limit cycle) and in Figure 6B, k_m=0.1475 (oscillatory regime, large limit cycle). The parameter k_m is key to the feedback regulation as it parametrizes the positive feedback of ComK onto MecA (see Methods). Changes in k_m affect the shape of the MecA nullcline resulting in a slight shift in the location of fixed point #3. Computed trajectories are shown in magenta. They originate at the location of the black dot. Note that the direction of motion (here, counterclockwise) is opposite that of the canonical circuit (clockwise, see Figure 2, for example). Figure 6C shows the bifurcation analysis of the 2D system vs. k_m. It depicts the overall dynamics and, and in particular explains the sudden appearance of the large limit cycle observed between k_m=0.1474 and k_m=0.1475. There is a Hopf bifurcation denoted by “H” at k_m ~ 0.1559 where oscillations are born. The ComK and MecA locations of fixed point #3 are indicated by black and orange lines, respectively. Solid lines are used to indicate stability while dotted lines are for unstable. As k_m decreases, the amplitude of the limit cycle born at the Hopf increases as is denoted by showing the limit cycle’s upper and lower limits: purple and blue lines for MecA; green and orange lines for ComK. On the left side of the diagram, these oscillations terminate abruptly. Again, as in the case of the canonical B. subtilis regulation circuit, we refer to this point at the oscillation appearance/disappearance point depending on whether crossing is from the left to right (increasing k_m) or opposite direction. In fact these oscillations appear in a similar dynamical manner as for the canonical circuit, that is to say, they appear at the saddle node, hence similar to Figure 2C). However, the detailed view for the inverted circuit is omitted for brevity.

Figure 7A–7C show corresponding stochastic simulations at k_m= 0.1, k_m = 1.5 and k_m = 0.2 respectively. Note that since units on the stochastic simulation plots are dimensioned, a scaling factor “sf” appears such that the effective k_m = k_m× sf (see Methods). Supplementary Figure S2 details how these probability density plots are obtained from combining a large number of statistically independent histories (see Methods). In accordance to the 2D bifurcation analysis shown in Figure 6C, we observe the stochastic simulations shown in Figure 7A–7C to be in excitatory, oscillatory and bi-stable regimes, respectively. Thus, there is general dynamical agreement between the two viewpoints: deterministic and stochastic. However, the deterministic dynamics was computed, by design, in the limit of the mRNA sub-manifold being at rest (infinite time-scale separation) whereas, the stochastic regimes shown in Figure 7 do not strictly enforce this limit. In Figure 7, the stochastic simulations have finite time-scale separation, not infinite. Below we focus on some detailed consequences of this difference: the two viewpoints (deterministic and stochastic) do, in fact, lack specific dynamical agreement.

Persisting oscillations in finite time-scale separation

We may ask how the finite time-scale separation affects the stochastic dynamics by focusing on the oscillation appearance/disappearance point (left side of Figure 6C). For reference, in Figure 8A we show the finite time-scale separation stochastic simulation at k_m = 0.125, that is to say, in robust excitatory regime (refer to Figure 6C. There is no limit cycle in the dynamics. In Figure 8B and 8C, k_m = 0.1474 and k_m = 0.1475, respectively. From the 2D bifurcation analysis, at k_m = 0.1474 there should be no oscillations, and at k_m =0.1475 large amplitude oscillations are expected. However, Figure 8B and 8C of the stochastic dynamics actually depict similar oscillatory behavior at these two values of k_m: the clear oscillations around the upper fixed point seen in Figure 8C also persist to the lower k_m regime of Figure 8B. We seek the source of the discrepancy.

Focusing on the time-scale separation difference between the 2D computation and the stochastic simulation, we show in Figure 8D a stochastic simulation that now includes a translation slowdown factor SDT= 0.5 (common to both ComK and MecA) (see Methods). Comparing Figure 8B and 8D, it is emphatically clear that increasing time-scale separation makes the limit cycle disappear from the stochastic simulation: there is no cycling about the upper fixed point anymore. This behavior is as expected as the stochastic simulation’s time-scale separation is now closer to that of the 2D ideal. Could stochasticity alone account for the difference? Supplementary Figure S3 shows a stochastic simulation without translation slowdown, but for which the number of molecules and volume are simultaneously increased by a factor of 10 thus leaving concentrations unchanged but reducing intrinsic fluctuations by about a factor of 3. A factor of 10 in number of molecules is arguably biologically maximal. Comparing Supplementary Figure S3A–S3C to the corresponding ones of Figure 8, we conclude that reduced biochemical noise alone cannot be the source of disagreement for the location of the oscillation appearance/disappearance point in the stochastic simulations. In fact, as we have shown above, it is time-scale separation that constitutes the source of disagreement (again, compare Figure 8B and 8D). Hence, in the inverted B. subtilis system, finite time scale separation extends the oscillatory regime range, or, from a biological point of view, finite time scale separation reduces the range of excitability.

DISCUSSION

In this work, we tackled the difficult problem of pinning down the respective roles of time-scale separation and of biochemical fluctuations in situations where these seemingly unrelated aspects of the dynamics conspire to affect the competence phenotype in Bacillus subtilis, and might be individually or collectively harnessed by Evolution. We took the approach that Evolution would harness the power of any “handle” and thus we investigated seemingly unrelated phenotype-breaking aspects of the core regulation of competence. We focused on parameterization departures from the living state, rather than global parameter exploration, thus offering a unique perspective. Comparing Figure 2 and Figure 6, we noticed that both the canonical and the inverted circuits share a hitherto unchartered phenotype flip ability from the excitability-based competence phenotype to an oscillation-based competence-breaking phenotype, and that this ability is within the reach of only a very slight feedback parameter adjustment. Hence we focused our investigations on this subtle but biologically impactful dynamical aspect to compare circuits. We devised an effective method to remove fluctuations from the stochastic simulations of the canonical circuit without having to conduct a prohibitive number of discrete-event simulations to estimate the mean of their probability distributions. This approach yields a system of coupled ordinary differential equations from the set of noisy discrete-event reactions. This born stochastic but now deterministic system conveniently continues to be endowed in the time scale separation underlying the stochastic dynamics. This approach simplifies studying the two conspiring aspects of the stochastic dynamics by separating their effects.

In the case of the canonical B. subtilis regulation (Figure 1A), we could then compare a 6D deterministic computation (with finite time-scale separation) to the 2D simplified view (with built-in infinite time-scale separation, the so-called adiabatic approximation). This 6D system is the parent dynamical system to the 2D system. We stress that we based our investigations on departures from a biologically proven parameterization of our models [²,⁵,⁶], so we did not rely on parameter sampling. By performing a bifurcation analysis of the 6D system with respect a “translation slowdown” parameter (keeping in mind that, in fact, this parameter parameterizes the slowdown of the entire protein sub-manifold with respect to transcription thus enables tunable time-scale separation), we looked for and found effective (unequal) translation slowdown factors in the ComK and ComS protein dynamics that could then be applied to the stochastic dynamics. The anticipated result however – the appearance of oscillations – did not have to immediately occur in a version of the stochastic dynamics now endowed with this slowed down translation. This is because stochasticity is, by design, not part of the 6D model and, yet, it could still be the source of the dynamical discrepancy. The anticipated oscillations did not, in fact, immediately occur in the translation slowed down version of the stochastic dynamics. But by concurrently reducing the amount of intrinsic biochemical noise in the system, we then discovered that this behavior is due to the conspiring effect of stochasticity. Hence, only by further reducing the size of biochemical fluctuations in the system, did the expected oscillations finally appear (Figure 5C). Note that reducing stochasticity alone in the system (leaving time-scale separation intact) was verified to be insufficient at resolving the dynamical discrepancy (data omitted for brevity).

In the case of the inverted B. subtilis regulation (Figure 1B), we showed, based again on departures from an experimentally verified parameterization [²], that oscillations that should have already vanished from an infinite time-scale separation perspective persisted only due to remaining finite time-scale separation in the stochastic system. Thus by slowing down translation on both the ComK and MecA protein manifold (here, by the same factor) in the discrete event simulation, oscillations expected to have vanished did, in fact, then vanish (compare Figure 8B and 8D). Additionally, we showed that stochasticity alone could not resolve the dynamical mismatch (Supplementary Figure S3).

In this work, we chose to use bifurcation analysis to aid in the exploration of the robustness of phenotype to parameter changes. We specifically focused on the onset of oscillations as phenotype breaking; excitability being competence phenotype-enabling. We explored core regulatory topology choices in the form of the canonical core regulation circuit and the inverted regulation circuit. We studied the impact of the operational level of biochemical noise dictated by evolutionary choices in the number (i.e., expression level) of key molecular regulators of competence present in the cell. Our modus operandi was thus to study the phenotypic impact of departures from biologically established model parameterization of the living state. In the seminal work of Ma et al.[¹¹], both network topology and parameter space were systematically numerically sampled for their respective and combined impact on the phenotype of adaptation. Applying the approach of automated massive parameter sampling, Zhang et al. [¹²] focused their own inquiry on the B. subtilis canonical and inverted core regulation circuits we studied here. Their key finding is that the canonical competence regulation circuit is more robust to parameter changes than the inverted circuit. The findings reported here are consistent with theirs since the appearance of phenotype-breaking oscillations required both the increase of time-scale separation in the canonical circuit as well as the decrease of intrinsic noise. In the case of the inverted circuit, only time-scale separation adjustment was needed for phenotype reconciliation, making it more brittle.

Suel et al. [⁶] already computationally explored the various dynamical regimes the canonical competence circuit presents under variation of some of the key regulation parameters accessible experimentally, but not the positive and negative feedback coupling parameters considered here. Robustness to all parameter variations and to noise in the competence durations was also studied in Ref. [²]. In contrast, in the present study, we focused on how both the canonical and evolutionary plausible inverted circuit can break the competence phenotype by switching from competence-enabling excitability into dynamically nearby competence-breaking oscillations upon the feedback coupling parameter perturbation, time-scale separation perturbation and intrinsic biochemical noise alteration. Our investigations can be considered local perturbations upon the living state and are complementary to the study of Zhang et al. [¹²]who used a global parameter variation approach to show that while the two networks are topologically equivalent, the corresponding systems can present different dynamical behaviors including the coexistence of stable steady state and oscillation, or two kinds of oscillations.

Is stochasticity truly an orthogonal evolutionary handle to time-scale separation?

It is conceivable that Evolution can adjust relevant constituent molecular numbers to affect stochasticity. It is also similarly conceivable that Evolution can independently affect time-scale separation by adjusting relevant rate parameters i.e., fine tuning binding constants [¹³–¹⁷] However, are these two evolutionary handles necessarily orthogonal? Evolution has of course no effect on basic physics. Here we note that cell size couples these two aspects. Physically larger cells have more molecules in a concomitant proportionately bigger volume so they are expected to be less noisy. This expectation has been confirmed experimentally by the fusion of two bacterial cells together [⁶]. However, because transcription and translation are physically separate non-co-located processes, physically larger cells would also be expected to have, on average, more time-scale separation. Therefore, for larger cells, physics suggests that noise would be reduced and time-scale separation would be increased. These are exactly the conditions we have observed in this work that led to phenotype breaking in the canonical circuit: increased time-scale separation and reduced noise caused a dynamical regime change from biologically desirable excitability to biologically undesirable oscillations. It would be interesting to experimentally probe the suggested relationship between the physical size of B. subtilis and the competence phenotype.

CONCLUSIONS

Time-scale separation is controlled by relationships between kinetic parameters in the mRNA and protein system whereas intrinsic stochasticity is controlled by the number of molecules of key regulators in the system. Therefore, these two dynamical aspects offer accessible handles (“knobs”) hat Evolution may be able to adjust independently to affect phenotype: e.g. as studied here, cause or avoid oscillations. In this work, for one evolutionary realized canonical bacterium and one hypothetical but evolutionary plausible re-engineered mutant, we specifically point out how these handles can be used to affect the feedback coupling parameter range over which oscillations occur. Elsewhere [²] it has already been shown that oscillations are detrimental to the competence phenotype. Here, consistent with Zhang et al. [¹²], we show that biologically undesirable oscillatory regimes are, in fact, present in both versions of the control circuits (canonical and inverted). Moreover, for the canonical competence circuit, onset of possibly phenotypically catastrophic large oscillations can be guarded against using both time-scale separation and fluctuation levels originating in the paucity of some of the molecular components in the system. On the other hand, for the inverted competence circuit, time-scale separation alone is sufficient to guard against unwanted oscillations. Our results reinforce the notion that the canonical regulation architecture as implemented by the living system is more robust to phenotype-breaking evolutionary choices yielding oscillations.

METHODS

Canonical and inverted B. subtilis circuits and the feedback continuation parameters

The deterministic and stochastic models of the canonical and inverted B. subtilis circuits are thoroughly developed by Cagatayet al. [²]. The canonical (i.e., native) B. subtilis core circuit was devised to explain the “Meks” bacterium behavior whereas the inverted B. subtilis circuit was devised to shed light on the “Synex” mutant behavior. In Ref. [⁶], the deterministic and stochastic models of the canonical B. subtilis circuits are also presented in details. In this work we reuse exactly the same models, equations and parameters of the living state. We see the present work as extending the original effort.

In the specific cases of the continuation parameter k₁ of the deterministic canonical model and the corresponding parameter k_m of the deterministic inverted model, we varied these parameters over a range to produce the relevant bifurcation analyses. More specifically, the role of k₁ is given by Equation S3 in the supplement of Ref. [⁶] where it parameterizes the negative feedback onto the protein ComS by the protein ComK. The role of k_m is similarly described by Equation S² in the supplement of Ref. [²]. Here the regulation is that of positive feedback onto the protein MecA by the protein ComK. These models are referred to in this work as 2D models because the dynamics is collapsed onto the protein dimensions via adiabatic approximations of the concomitant mRNA etc. sub-manifolds.

6D model: the infinite number of molecule limit of the underlying stochastic model of the canonical B. subtilis circuit

The equations for the 6D model are straightforwardly assembled, term by term (production terms are positively signed, degradation terms are negatively signed) from the set of discrete-event reactions of Ref. [²] and are given below.

dmRNA_ComK /dt= Pconstitutive_ComK × k₁ + P_ComK × f(ComK,k₂,k_k,n) − mRNA_ComK × k_7,

dmRNA_ComS /dt= Pconstitutive_ComS × k₄ + P_ComS × g(ComK,k₅,k_s,p)− mRNA_ComS × k_9,

dComK /dt= mRNA_ComK × k₃ + MecA_ComK × km₁₁ − ComK × k₈ − MecA_free × ComK × k_11,

dComS /dt= mRNA_ComS × k₆ + MecA_ComS × km₁₃ − ComS × k₁₀ − MecA_free × ComS × k_13,

dMecA_ComK /dt= MecA_free × ComK x k₁₁ − MecA_ComK × km₁₁ − MecA_ComK × k_12,

dMecA_ComS /dt= MecA_free × ComS x k₁₃ − MecA_ComS × km₁₃ − MecA_ComS × k_14,in which

f (ComK,k₂,k_k,n) = (k₂ × ComKⁿ) / (k_kⁿ + ComKⁿ),

g (ComK,k₅,k_s,p) = k₅ / (1+ (ComK/k_s)^p),

MecA_free = (M_T − MecA_ComK − MecA_ComS).

Because these equations have no innate fluctuations, they represent the infinite number of molecules version of the underlying stochastic model. They are expected to reproduce the time evolution of the first moment of the stochastic model state variable distributions that would be obtained with an infinite number of statistically independent histories. The parameters of the 6D model are consistent with the stochastic model in Ref. [²]. For convenience, we reproduce them below; in dimensioned units:

k₁: 0.0002275 s⁻¹

k₂: 0.195 s⁻¹

k₃: 0.2 s⁻¹

k₄: 0

k₅: 0.00156 s⁻¹

k₆: 0.2 s⁻¹

k₇: 0.005 s⁻¹

k₈: 0.0001 s⁻¹

k₉: 0.005 s⁻¹

k₁₀: 0.0001 s⁻¹

k₁₁: 2.02e-06 molec⁻¹ s⁻¹

km₁₁: 0.00252 s⁻¹

k₁₂: 0.05 s⁻¹

k₁₃: 4.5e-06 molec⁻¹ s⁻¹

km₁₃: 5.36e-05 s⁻¹

k₁₄: 4e-05 s⁻¹

k_k: 5200 molec

k_s: 1300 molec

n: 2

p: 3

Pconstitutive_ComK : 1

Pconstitutive_ComS: 1

P_ComK: 1

P_ComS: 1

M_T: 520 molec

The simulation volume is W = 1 molec / nMolar.

The above equations produce molecular numbers. Concentrations are obtained by dividing the number of molecules by the simulation volume e.g., [ComK] = ComK / W .The units of concentration are therefore nanomolar. Dimensionless concentrations discussed in this work are obtained from the dimensioned concentrations using scaling factors:

[ComK] _{dimensionless} = [ComK] / G_ComKG_ComK = 2.6×10⁴,

[ComS] _{dimensionless} = [ComS] / G_ComSG_ComS = 20.8.

Slowdown of translation in the canonical B. subtilis circuit

The stochastic model is described at length in Ref. [²]. The details of the model will not be reproduced here, for brevity. The slowdown of translation in the stochastic model as well as in its deterministic 6D infinite number of molecule limit model described in the preceding paragraph is implemented as follows:

k₃ ← k₃ × sdt_ComK Translation of mRNA_ComK,

k₈ ← k₈ × sdt_ComK Degradation of ComK,

k₁₁ ← k₁₁ × sdt_ComKComK binding free MecA,

km₁₁ ← km₁₁ × sdt_ComKMecA-ComK dissociation into MecA and ComK,

k₁₂ ← k₁₂ × sdt_ComKProteolytic degradation of MecA-bound ComK,

k₆ ← k₆ × sdt_ComS Translation of mRNA_ComS,

k₁₀ ← k₁₀ × sdt_ComS Degradation of ComS,

k₁₃ ← k₁₃ × sdt_ComSComS binding to free MecA,

km₁₃ ← km₁₃ × sdt_ComSMecA-ComS dissociation into MecA and ComS,

k₁₄ ← k₁₄ × sdt_ComSProteolytic degradation of MecA-bound ComS.

Slowdown of translation in the inverted B. subtilis circuit

The stochastic model is described at length in Ref. [²]. The slowdown of translation in this model is implemented as follows:

k₅ ← k₅ × sdt Translation of mRNA ComK,

k₆ ← k₆ × sdt Translation of mRNA MecA,

k₁₀ ← k₁₀ × sdt Degradation of MecA,

k₁₁ ← k₁₁ × sdt Degradation of ComK.

Scaling in the inverted B. subtilis circuit

The scaling factor “sf” relating dimensioned to dimensionless concentrations ([ ] _dimensioned = [ ]_{dimensionless} × sf) in the inverted B. subtilis circuit is sf= 25,000 [⁴].

Stochastic simulations

Stochastic phase portraits in this work are occupancy diagrams (i.e., 2D histograms) generated by combining multiple (typically 200) stochastically independent histories binned over the range of the state variables. Each time history consists of an extended trajectory simulated by use of the Gillespie algorithm [¹⁸–²⁴].

Bifurcation analysis

Bifurcation analyses were performed using Auto [²⁵], graphical interfaces XPP [²⁶] and Oscill8 [²⁷].

Simulation framework and computing platforms

All simulations were developed and performed using Matlab (The MathWorks, Natick, Massachusetts) on a Dell Intel i7-2640M dual-core Windows 7 laptop, and a Lenovo Intel i7-7700HQ quad-core laptop running Windows 10. Simulations required approximately ~3000 h (~four full-time months) of integrated time over the course of a year.

References

[1] Eldar, A. and Elowitz, M. B. (2010) Functional roles for noise in genetic circuits. Nature, 467, 167–173

[2] Cağatay, T., Turcotte, M., Elowitz, M. B., Garcia-Ojalvo, J. and Süel, G. M. (2009) Architecture-dependent noise discriminates functionally analogous differentiation circuits. Cell, 139, 512–522

[3] SchultzD., Ben Jacob E., OnuchicJ. N. and WolynesP. G. (2007) Molecular level stochastic model for competence cycles in Bacillus subtilis. Proc. Natl. Acad. Sci. USA., 104, 17582–17587

[4] TurcotteM., Garcia-Ojalvo J. and SüelG. M. (2008) A genetic timer through noise-induced stabilization of an unstable state. Proc. Natl. Acad. Sci. USA, 105, 15732–15737.

[5] SüelG. M., Garcia-Ojalvo J., LibermanL. M. and ElowitzM. B. (2006) An excitable gene regulatory circuit induces transient cellular differentiation. Nature, 440, 545–550.

[6] SüelG. M., Kulkarni R. P., DworkinJ., Garcia-OjalvoJ. and ElowitzM. B. (2007) Tunability and noise dependence in differentiation dynamics. Science, 315, 1716–1719.

[7] DubnauD. (1999) DNA uptake in bacteria. Annu. Rev. Microbiol., 53, 217–244.

[8] GrossmanA. D. (1995) Genetic networks controlling the initiation of sporulation and the development of genetic competence in Bacillus subtilis. Annu. Rev. Genet., 29, 477–508.

[9] DubnauD. (1991) Genetic competence in Bacillus subtilis. Microbiol. Rev ., 55,395

[10] DubnauD. (1991) The regulation of genetic competence in Bacillus subtilis. Mol. Microbiol., 5, 11–18.

[11] MaW., Trusina A., El-SamadH., LimW. A. and TangC. (2009) Defining network topologies that can achieve biochemical adaptation. Cell, 138, 760–773.

[12] ZhangJ., Yuan Z., LiH. X. and ZhouT. (2010) Architecture-dependent robustness and bistability in a class of genetic circuits. Biophys. J., 99, 1034–1042.

[13] LevyE. D. (2010) A simple definition of structural regions in proteins and its use in analyzing interface evolution. J. Mol. Biol., 403, 660–670.

[14] KastritisP. L. and BonvinA. M. (2013) Molecular origins of binding affinity: seeking the Archimedean point. Curr. Opin. Struct. Biol., 23, 868–877.

[15] KastritisP. L., Rodrigues J. P. G. L., Folkers G. E., BoelensR. and BonvinA. M. J. J. (2014) Proteins feel more than they see: fine-tuning of binding affinity by properties of the non-interacting surface. J. Mol. Biol., 426, 2632–2652.

[16] RosanovaA., Colliva A., OsellaM. and CaselleM. (2017) Modelling the evolution of transcription factor binding preferences in complex eukaryotes. Sci. Rep., 7, 7596

[17] EchaveJ. and Wilke C. O. (2017) Biophysical models of protein evolution: understanding the patterns of evolutionary sequence divergence. Annu. Rev. Biophys., 46, 85–103.

[18] GillespieD. T. (1976) A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. J. Comput. Phys., 22, 403–434.

[19] GillespieD. T. (1977) Exact stochastic simulation of coupled chemical reactions. J. Phys. Chem., 81, 2340–2361.

[20] GillespieD. T. (2001) Approximate accelerated stochastic simulation of chemically reacting systems. J. Chem. Phys., 115, 1716–1733.

[21] GillespieD. T. (1976) A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. J. Comput. Phys., 22, 403–434.

[22] GillespieD. T. (1977) Exact stochastic simulation of coupled chemical reactions. J. Phys. Chem., 81, 2340–2361.

[23] GillespieD. T. (1991) Markov Processes: An Introduction for Physical Scientists. Manhattan: Academic Press

[24] GillespieD. T. (2007) Stochastic simulation of chemical kinetics. Annu. Rev. Phys. Chem., 58, 35–55.

[25] DoedelE. J. (1986) AUTO: software for continuation and bifurcation problems in ordinary differential equations. California Institute of Technology, 12, 791–802

[26] ErmentroutB. (2011) XPP-Aut

[27] ConradE. D. (2011) Oscill8

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag GmbH Germany, part of Springer Nature