Network design principle for robust oscillatory behaviors with respect to biological noise
Abstract
Oscillatory behaviors, which are ubiquitous in transcriptional regulatory networks, are often subject to inevitable biological noise. Thus a natural question is how transcriptional regulatory networks can robustly achieve accurate oscillation in the presence of biological noise. Here, we search all two and threenode transcriptional regulatory network topologies for those robustly capable of accurate oscillation against the parameter variability (extrinsic noise) or stochasticity of chemical reactions (intrinsic noise). We find that, no matter what source of the noise is applied, the topologies containing the repressilator with positive autoregulation show higher robustness of accurate oscillation than those containing the activatorinhibitor oscillator, and additional positive autoregulation enhances the robustness against noise. Nevertheless, the attenuation of different sources of noise is governed by distinct mechanisms: the parameter variability is buffered by the long period, while the stochasticity of chemical reactions is filtered by the high amplitude. Furthermore, we analyze the noise of a synthetic human nuclear factor κB (NFκB) signaling network by varying three different topologies, and verify that the addition of a repressilator to the activatorinhibitor oscillator, which leads to the emergence of highrobustness motif—the repressilator with positive autoregulation, improves the oscillation accuracy in comparison to the topology with only an activatorinhibitor oscillator. These design principles may be applicable to other oscillatory circuits.
Article activity feed

Author Response
Reviewer #2 (Public Review):
The manuscript by Qiao et al studies the important problem of how to achieve accurate oscillation robustly in biological networks where noise level may be high. The authors adopted a comprehensive approach and studied how different network configurations affect osciillation. This is based on enumeration of all major architectures of 2 and 3nodes networks.
This work makes important contributions to the field, as it offers the first comprehensive survey of networks motifs capable of oscillation, with further characterization of their robustness. The authors identified core motifs of repressilator with a positive autoregulation, and activatorinhibitor oscillator. In addition, the authors have identified different mechanisms of attenuation of different sources of noises. Overall, this is …
Author Response
Reviewer #2 (Public Review):
The manuscript by Qiao et al studies the important problem of how to achieve accurate oscillation robustly in biological networks where noise level may be high. The authors adopted a comprehensive approach and studied how different network configurations affect osciillation. This is based on enumeration of all major architectures of 2 and 3nodes networks.
This work makes important contributions to the field, as it offers the first comprehensive survey of networks motifs capable of oscillation, with further characterization of their robustness. The authors identified core motifs of repressilator with a positive autoregulation, and activatorinhibitor oscillator. In addition, the authors have identified different mechanisms of attenuation of different sources of noises. Overall, this is an important study reporting many new results.
 The current stochastic model is based on the deterministic model shown in Fig 1D. However, there are a lot of assumptions in this ODE model. These include the assumption that the substrate is in instantaneous chemical equilibrium with the proteinligand/DNA complex, and the reaction involving one receptor and n identical simultaneously binding ligands, with the Hill coefficient phenomenologically characterizing cooperativity. However, it is not clear what the specifics are when applied to the networks studied, i.e., what reactions are assumed instantaneous equilibrium, or why the important phenomenon of slow TF and promoter binding can be ignored and why that is reasonable? Also, what are the receptors and what are the binding reactions of multiple ligands that gives to the Hill coefficient?
Thank you for the good suggestion. To clarify the model assumptions, we updated the Methods section Deterministic model.
As the reviewer pointed out, these assumptions may not hold in some biological systems, e.g., the case that the binding/unbinding of TF to the promotor is slow. In the revision, we added discussions on the limitations of our methods on Page 23.
 The described kinetic models based on ODE approximations may not be applicable to study strong intrinsic stochasticity arising at low copy number of molecules, where the MichelisMenton/Hilltype of ODE models are not valid. A more straightforward model would be based on mass action but more detailed reactions, without additional assumptions, from which one can write down the corresponding master equation.
In that light, it will be helpful to write out the set of equations for the stochastic networks, which is equivalent to the set of chemical master equations, from which Gillespie simulation samples.
Thank you for the suggestion. In the manuscript, we simulated the system by using Gillespie algorithm with Hilltype reaction rates (see Section “Simulations using the Gillespie algorithm lead to similar conclusions” and Figure S2). The simulation results are similar to those using Langevin equations (Figure 3 for results using Langevin equations). We updated the section "Stochastic model in the presence of intrinsic noise" to describe how we use Gillespie algorithm.
However, using Hilltype reaction rates in stochastic simulations may not be always accurate and could fail to capture the role of detailed reactions in stochasticity. Thus, we anticipate the need for a more detailed model where every reaction of Hilltype form is decomposed into the elementary reactions. The related discussions were added in the revision (see also the response to the 3rd question).
 There may be another issue with the stochastic model for the intrinsic noise, as the reaction system do not model some of the important stochasticity occurring in the system. To study transcriptional regulation under different molecules, the binding of transcription factor to the DNA/promoter is an important source of stochasticity, as copy numbers of that specific regulating protein TF may vary, at the same time it also regulates the production rate of another protein. This process cannot be adequately modeled only by copy number change of the regulated protein. That is, the binding and unbinding of the transcription factor to the promoter plays important roles in stochasticity when modeling the intrinsic noise of the biochemical reactions. It would be desirable to account for the process explicitly. If the authors decide not to model such stochasticities, potential caveats should be pointed out.
In addition, proteingene interactions often involve dimerization, that is, gene product forms a dimer, which then interact with the other gene. This dimerization may also qualitatively affect stochastic behavior of a GRN. The authors may wish to discuss all these issues.
Thank you for the good suggestion. In the revision, we added the following paragraph in the Discussion: “Another limitation of our work is that we didn’t decompose the reactions in the deterministic model into detailed elementary reaction steps when using Gillespie algorithm. The advantage of simulating nonelementary reactions with Hilltype rate functions is the low computation cost, and in some biological networks, it leads to consistent results with the model composed of all elementary reactions (Gonze et al., 2002; Kim et al., 2014; Sanft et al., 2011). However, this approach may not be always accurate, depending on the timescale separation of reactions (Kim et al., 2014; Sanft et al., 2011); for example, the Hilltype reaction rate is based on the quasisteadystate approximation, which does not hold when binding/unbinding of TF to the promotor is slow or comparable to the timescales of protein production or decay (Choi Paul et al., 2008). Furthermore, this method neglects detailed reactions in gene regulatory networks, and thus fails to study roles of these reactions in stochasticity. These detailed reactions include the binding and unbinding of the transcription factor to the promoter, dimerization of transcription factors, transcription and translation (Cao et al., 2018; Shahrezaei and Swain Peter, 2008; Terebus et al., 2019).”
 A potential drawback of the study is that the oscillation behavior hinges upon the behavior of the ODE deterministic model. It is well known even for simple networks such as transcription regulation without feedback, or when protein binding is involved there can be significant divergency between ODE model and stochastic model, where the latter exhibit multistabilities and the former none (e.g., doi.org/10.1063/1.3625958, doi.org/10.1103/PhysRevE.91.042111, doi.org/10.1063/1.5124823). There is now an increasing body of literature documenting this. This issue and potential ramifications should be discussed.
As an example, there is a new mechanism of stochastic oscillation found in toggle switch under weak promoter binding condition. This is not obvious from the corresponding ODE model and requires computation of the global map of discrete flux (doi.org/10.1063/1.5124823). It will be missed if proteinDNA binding is not modeled explicitly. It will be interesting if the authors can discuss the relationship of this type of oscillation with those based on repressilator/autoregulation and activatorinhibitor. Do they belong to perhaps different class of stochastic oscillations and if so, what are the differences?
Thank you for the good suggestion. In the revision, we added one paragraph on Page 22 to discuss this potential drawback: “Our work only focused on the effects of biological noise on oscillation accuracy, neglecting other dynamic changes caused by noise. These dynamics may include the loss of multistability and the occurrence of oscillation. Specifically, the way to model the noise may cause the loss of multistability (Duncan et al., 2015; Vellela and Qian, 2009); the presence of noise can produce oscillation even when the corresponding deterministic model cannot oscillate, which has been validated in the toggleswitch system and excitable system (Lindner et al., 2004; Terebus et al., 2019; Zaks et al., 2005). The possible reason might be the noiseinduced transition between different states. Since our work only studied network topologies whose deterministic model can generate oscillation, we did not count the topologies that cannot oscillate in the deterministic model but begin to oscillate in the stochastic model. Due to the popularity of such topologies, how these topologies buffer noise will be of interest and may lead to the discoveries of new principles.”
 The authors decided to use Chemical Langevin Equation to model the stochastic process due to computational cost. However, recent development shows that computing cost may no longer be an issue, as the finite buffer ACME algorithm can generate full probability surfaces without running costly trajectories (doi: 10.1073/pnas.1001455107, doi: 10.1137/15M1034180). In fact, this has been done for 3node networks (feedforward loops), where extensive parameter sweeps enables construction of 10^4 probability landscape, from which phase diagrams of multimodality can be constructed (doi.org/10.3389/fgene.2021.645640). I understand that the authors choose to use the Langevin model, but the probability surface governed by the chemical master equation can now be computed rather rapidly, without resorting timeconsuming Gillespie simulations. Therefore, the rationale of high computing cost may not be justifiable. This advancement should be pointed out.
Thank you for the good suggestion. In the revision, we added the following text to introduce the advancement of algorithms in the Discussion: “We anticipate the need for a more detailed model where every reaction of Hilltype form is decomposed into the elementary reactions. The recent development about stochastic algorithms with fast computation makes it feasible to simulate such detailed model for all two and threenode network topologies, for example, algorithms focusing on solving the chemical master equations (Cao et al., 2010; Cao et al., 2016; Munsky and Khammash, 2006; Terebus et al., 2021) and variants of Gillespie algorithms that directly simulate the temporal dynamics (Gillespie and Petzold, 2003). Besides, the construction of probability surfaces through these algorithms may shed light on new principles for accurate oscillation.”
 For the reason that there are many differences between master equation model, Langevin model, and ODE model, the statement on p.8 “the system responds to the noise is usually linked to the deterministic features” should be modified/qualified.
Sorry for the confusion. To clarify this, in the revision, we expanded this sentence to the following texts on Page 9: “Note that how the system responds to the noise is often linked to the deterministic features (Monti et al., 2018; Paulsson, 2004; Wang et al., 2010). For example, Monti et al. found that the circuit’s ability to sense time under input noise becomes worse when this circuit’s deterministic behavior cannot generate the limit cycle; Wang et al. adopted a similar form of noise and demonstrated the importance of signed activation time, a quantity calculated based on deterministic behavior, on the noise attenuation; by using an Ωexpansion to approximate the birthanddeath Markov process, Paulsson obtained the variance of the protein in gene networks and found it is related to network’s elasticity which is calculated from the deterministic model.”

Evaluation Summary:
The authors study the important problem of how to achieve accurate oscillation robustly in biological networks where noise level may be high. The authors adopted a comprehensive approach and study how different network configurations affect oscillation. This work makes an important contribution to the field, as it offers the first comprehensive survey of networks motifs capable of oscillation, with further characterization of their robustness.
(This preprint has been reviewed by eLife. We include the public reviews from the reviewers here; the authors also receive private feedback with suggested changes to the manuscript. The reviewers remained anonymous to the authors.)

Reviewer #1 (Public Review):
Oscillations are common in gene regulatory networks, giving cell cycles, circadian clocks etc. They are often under the influence of biological noise. By varying different topologies, the authors found that the topologies containing the repressilator with positive autoregulation show higher robustness of accurate oscillation than those containing the activatorinhibitor oscillator, and additional positive autoregulation enhances the robustness against noise. They further found that the parameter variability from extrinsic noise is buffered by the long period, while the stochasticity of chemical reactions from intrinsic noises is filtered by the high amplitude. This demonstrates the crucial role of the topology of the network for influencing the quality of the oscillations. These may serve as design …
Reviewer #1 (Public Review):
Oscillations are common in gene regulatory networks, giving cell cycles, circadian clocks etc. They are often under the influence of biological noise. By varying different topologies, the authors found that the topologies containing the repressilator with positive autoregulation show higher robustness of accurate oscillation than those containing the activatorinhibitor oscillator, and additional positive autoregulation enhances the robustness against noise. They further found that the parameter variability from extrinsic noise is buffered by the long period, while the stochasticity of chemical reactions from intrinsic noises is filtered by the high amplitude. This demonstrates the crucial role of the topology of the network for influencing the quality of the oscillations. These may serve as design principles applicable to many oscillatory circuits.

Reviewer #2 (Public Review):
The manuscript by Qiao et al studies the important problem of how to achieve accurate oscillation robustly in biological networks where noise level may be high. The authors adopted a comprehensive approach and studied how different network configurations affect osciillation. This is based on enumeration of all major architectures of 2 and 3nodes
networks.This work makes important contributions to the field, as it offers the first comprehensive survey of networks motifs capable of oscillation, with further characterization of their robustness. The authors identified core motifs of repressilator with a positive autoregulation, and activatorinhibitor oscillator. In addition, the authors have identified different mechanisms of attenuation of different sources of noises. Overall, this is an important study …
Reviewer #2 (Public Review):
The manuscript by Qiao et al studies the important problem of how to achieve accurate oscillation robustly in biological networks where noise level may be high. The authors adopted a comprehensive approach and studied how different network configurations affect osciillation. This is based on enumeration of all major architectures of 2 and 3nodes
networks.This work makes important contributions to the field, as it offers the first comprehensive survey of networks motifs capable of oscillation, with further characterization of their robustness. The authors identified core motifs of repressilator with a positive autoregulation, and activatorinhibitor oscillator. In addition, the authors have identified different mechanisms of attenuation of different sources of noises. Overall, this is an important study reporting many new results.
The current stochastic model is based on the deterministic model shown in Fig 1D. However, there are a lot of assumptions in this ODE model. These include the assumption that the substrate is in instantaneous chemical equilibrium with the proteinligand/DNA complex, and the reaction involving one receptor and n identical simultaneously binding ligands, with the Hill coefficient phenomenologically characterizing cooperativity. However, it is not clear what the specifics are when applied to the networks studied, i.e., what reactions are assumed instantaneous equilibrium, or why the important phenomenon of slow TF and promoter binding can be ignored and why that is reasonable? Also, what are the receptors and what are the binding reactions of multiple ligands that gives to the Hill coefficient?
The described kinetic models based on ODE approximations may not be applicable to study strong intrinsic stochasticity arising at low copy number of molecules, where the MichelisMenton/Hilltype of ODE models are not valid. A more straightforward model would be based on mass action but more detailed reactions, without additional assumptions, from which one can write down the corresponding master equation.
In that light, it will be helpful to write out the set of equations for the stochastic networks, which is equivalent to the set of chemical master equations, from which Gillespie simulation samples.
There may be another issue with the stochastic model for the intrinsic noise, as the reaction system do not model some of the important stochasticity occurring in the system. To study transcriptional regulation under different molecules, the binding of transcription factor to the DNA/promoter is an important source of stochasticity, as copy numbers of that specific regulating protein TF may vary, at the same time it also regulates the production rate of another protein. This process cannot be adequately modeled only by copy number change of the regulated protein. That is, the binding and unbinding of the transcription factor to the promoter plays important roles in stochasticity when modeling the intrinsic noise of the biochemical reactions. It would be desirable to account for the process explicitly. If the authors decide not to model such stochasticities, potential caveats should be pointed out.
In addition, proteingene interactions often involve dimerization, that is, gene product forms a dimer, which then interact with the other gene. This dimerization may also qualitatively affect stochastic behavior of a GRN. The authors may wish to discuss all these issues.
A potential drawback of the study is that the oscillation behavior hinges upon the behavior of the ODE deterministic model. It is well known even for simple networks such as transcription regulation without feedback, or when protein binding is involved there can be significant divergency between ODE model and stochastic model, where the latter exhibit multistabilities and the former none (e.g., doi.org/10.1063/1.3625958,doi.org/10.1103/PhysRevE.91.042111,doi.org/10.1063/1.5124823). There is now an increasing body of literature documenting this. This issue and potential ramifications should be discussed.
As an example, there is a new mechanism of stochastic oscillation found in toggle switch under weak promoter binding condition. This is not obvious from the corresponding ODE model and requires computation of the global map of discrete flux (doi.org/10.1063/1.5124823). It will be missed if proteinDNA binding is not modeled explicitly. It will be interesting if the authors can discuss the relationship of this type of oscillation with those based on repressilator/autoregulation and activatorinhibitor. Do they belong to perhaps different class of stochastic oscillations and if so, what are the differences?
The authors decided to use Chemical Langevin Equation to model the stochastic process due to computational cost. However, recent development shows that computing cost may no longer be an issue, as the finite buffer ACME algorithm can generate full probability surfaces without running costly trajectories (doi: 10.1073/pnas.1001455107, doi: 10.1137/15M1034180 ). In fact, this has been done for 3node networks (feedforward loops), where extensive parameter sweeps enables construction of 10^4 probability landscape, from which phase diagrams of multimodality can be constructed (doi.org/10.3389/fgene.2021.645640). I understand that the authors choose to use the Langevin model, but the probability surface governed by the chemical master equation can now be computed rather rapidly, without resorting timeconsuming Gillespie simulations. Therefore, the rationale of high computing cost may not be justifiable. This advancement should be pointed out.
For the reason that there are many differences between master equation model, Langevin model, and ODE model, the statement on p.8 "the system responds to the noise is usually linked to the deterministic features" should be modified/qualified.
