Network design principle for robust oscillatory behaviors with respect to biological noise

Abstract
Editor's evaluation
Introduction
Results
Discussion
Methods
Data availability
References
Article and author information
Metrics

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 three-node 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 activator-inhibitor 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 activator-inhibitor oscillator, which leads to the emergence of high-robustness motif—the repressilator with positive autoregulation—improves the oscillation accuracy in comparison to the topology with only an activator-inhibitor oscillator. These design principles may be applicable to other oscillatory circuits.

Editor's evaluation

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.

https://doi.org/10.7554/eLife.76188.sa0

Introduction

Oscillatory behaviors have been observed in a broad range of biological processes, such as cell cycle (Ferrell et al., 2011; Tyson, 1991), circadian rhythms (Partch et al., 2014), and mitotic wave in Drosophila embryo (Deneke et al., 2016). Oscillatory features, including period and amplitude, can encode functional information, which plays an essential role in coordinating gene regulation (Cai et al., 2008) or transmitting distinct stimuli (Hao and O’Shea, 2012; Heltberg et al., 2019). In past decades, negative feedback, time delay, and nonlinearity have been identified as key mechanisms for biochemical oscillation (Novák and Tyson, 2008), following which researchers artificially synthesized biochemical networks capable of oscillation (Atkinson et al., 2003; Chen et al., 2015; Elowitz and Leibler, 2000; Potvin-Trottier et al., 2016; Stricker et al., 2008; Tigges et al., 2010; Zhang et al., 2017). Repressilator (Elowitz and Leibler, 2000) and activator-inhibitor oscillator (Atkinson et al., 2003) are the most famous of these synthetic oscillators.

While many synthetic biological circuits can oscillate, their dynamics are typically irregular, owing to ubiquitous biological noise such as fluctuations in the microenvironment and inherent stochasticity of chemical reactions (Elowitz et al., 2002; Li et al., 2009; Potvin-Trottier et al., 2016; Raser and O’Shea, 2004; Swain et al., 2002; Yu et al., 2018). Thus, a natural question is how the biological systems achieve accurate oscillation in the presence of noise. Previous studies revealed that many kinetic parameters can influence the robustness of the biological oscillators, such as the system size and degree of cooperativity of reactions (Gonze et al., 2002a), timescale of the promoter interaction (Forger and Peskin, 2005), repressor degradation rate (Potvin-Trottier et al., 2016), free energy cost measured by ATP/ADP ratios (Cao et al., 2015; Fei et al., 2018; Qin et al., 2021), and kinetic parameter-determined oscillation mechanisms (i.e., limit cycle or force driving) (Monti et al., 2018). Moreover, growing evidence suggests the existence of the relationship between network configurations and noise buffering capabilities for biochemical oscillators. For example, in a synthetic microbial consortium oscillator composed of two different types of bacteria, adding negative autoregulation to the negative feedback loop increases the parameter space to oscillate persistently in the face of noise (Chen et al., 2015); an additional positive feedback loop in the biochemical oscillator consisting of the negative feedback loop can decrease the coefficient of variation (CV) of period when considering the stochasticity of reactions (Mather et al., 2009) and possess nearly constant period when varying the synthesis rate (Stricker et al., 2008).

Instead of exploring mechanisms to achieve accurate oscillation case by case, we try to understand the general network design principles of accurate oscillation using the bottom-up approach (Ma et al., 2009; Qiao et al., 2019) and discover the specific network topologies that can oscillate and attenuate noise simultaneously. Here, we systematically explore the relationship between the network topology and robustness to different sources of noise in both two- and three-node networks. We first perform an exhausting search of two- and three-node network topologies to identify those capable of oscillation in the absence of noise, and then investigate the abilities of those oscillatory topologies to achieve accurate oscillation in the presence of different sources of noise. Two different sources are considered: parameters are perturbed by noise terms whose magnitudes are proportional to parameters (i.e., extrinsic noise); chemical reactions induce stochasticity due to a small copy number of proteins (i.e., intrinsic noise). We classify all oscillatory topologies according to what core motifs they include, and then compare the ability to execute accurate oscillation in the presence of noise among different categories. Two categories whose core motifs include a repressilator with a positive feedback perform better than others. Importantly, the existence of positive autoregulation always enhances the performance. While these results hold regardless of what source of noise exists, mechanisms to attenuate different sources of noise are distinct: long period buffers the extrinsic noise, and high amplitude attenuates the intrinsic noise. Moreover, we experimentally validate that adding a repressilator to the activator-inhibitor topology in synthetic NF-κB signaling circuits can improve the performance to buffer noise, indicating the important role of the repressilator with a positive autoregulation in filtering noise.

Results

Searching for topologies robustly executing accurate oscillation

Index for measuring the oscillation accuracy

To measure the accuracy of the oscillatory behavior, we use the dimensionless correlation time, which is the correlation time $τ$ divided by the period $T$ . The correlation time $τ$ describes how fast the autocorrelation function $C (t)$ exponentially decays. To be specific, for a noisy dynamic trajectory of the oscillator, $C (t)$ displays a damped oscillation (Figure 1A):

C (t) \equiv \frac{{〈(x (t + s) - 〈x〉) (x (s) - 〈x〉)〉}_{s}}{〈x^{2}〉 - {〈x〉}^{2}} = \exp (- t / τ) \times c o s (2 π t / T)

Figure 1

Download asset Open asset

Searching all possible two-node and three-node network topologies for oscillatory topologies with high accuracy in the presence of noise.

(A) The accuracy of the oscillatory behavior against the noise is measured by the ratio of the correlation time $τ$ to the period $T$ . (B) Possible links in two-node and three-node network topologies. (C) Two typical topologies of biological oscillators: the activator-inhibitor and the repressilator topologies. (D) The deterministic model. (E) Stochastic models where the extrinsic and intrinsic noise are considered separately. (F) Illustration of searching oscillatory topologies and measuring the robustness to noise for a given oscillatory topology. See ‘Methods’ and Supplementary file 1a for parameter ranges and magnitudes of extrinsic and intrinsic noise.

where ${〈\dots〉}_{s}$ is defined by ${〈f (s)〉}_{s} = \lim_{S \to \infty} \frac{1}{S} \int_{0}^{S} f (s) d s$ , and $〈\dots〉$ is the ensemble average; $T$ is the period (time needed from one peak to the next peak). If fluctuations of the noisy trajectory are small, the autocorrelation decays slowly, leading to a large value of $τ$ . The correlation time $τ$ has the same unit as that of the period, so $τ / T$ is dimensionless. Therefore, instead of using $τ$ , we utilize $τ / T$ to measure the accuracy, which is equal to the quantity that previous work has used (Cao et al., 2015) except a constant factor.

Network topology space

We limit ourselves to network topologies with two or three nodes (Figure 1B) and search for topologies capable of accurate oscillation using a bottom-up concept. While the signaling pathway of the oscillator in nature is complex, the core motif executing functions may be simple (Lim et al., 2013; Ma et al., 2009; Novák and Tyson, 2008; Qiao et al., 2019), and thus two- or three-node networks might be enough to capture key features. Besides, the number of all two- and three-node network topologies is 3⁹ because there are nine links in total and each link has three options: activation, inhibition, or does not exist; however, by excluding topologies with isolated nodes or symmetrical to existing topologies, the number of possible two- and three-node network topologies is reduced from 3⁹ to 1955. Here, two typical oscillatory topologies are shown in Figure 1C: the activator-inhibitor and repressilator topologies. For the activator-inhibitor topology, the activator (node A) has a positive autoregulation and positively regulates the inhibitor (node B), but is negatively regulated by the inhibitor; for the repressilator topology, each node acts as a repressor to inhibit its next node, thus constituting a cyclic negative-feedback loop.

Mathematical modeling

To model two- and three-node network topologies, we use transcriptional regulations to describe interactions among nodes (Figure 1D; see ‘Methods’). In a transcriptional regulatory network, nodes and links represent genes’ products and transcriptional regulations, respectively; genes’ products work as transcription factors to interact with the regulatory sequence of other genes and activate or inhibit the transcription, regulating the production rates of other genes’ products, that is, other nodes. Moreover, when multiple transcription factors regulate the same gene simultaneously, the competitive inhibition logic is adopted: those transcription factors compete for the same binding sites. Thus, the transcriptional activity of a gene depends on the relative weights of transcription factors activating this gene and those inhibiting this gene. Figure 1D illustrates the ordinary differential equation describing dynamics of node A when node A not only activates itself but also is inhibited by node B. In this equation, the variable $A$ represents the concentration of the product of gene A; $k_{b a s a l}$ is the basal production rate (much smaller than other terms); $v_{A A}$ is the maximum production rate caused by product A; $K_{A A}$ and $K_{B A}$ are binding affinities of products $A$ and $B$ to gene A, respectively; $r_{A}$ is the degradation rate; $n$ is the Hill coefficient; and the production rate is determined by relative weights of $A$ and $B$ .

Based on the above deterministic model, we develop stochastic models to describe the oscillatory behavior in the presence of noise. According to the source of noise, the biological noise can be decomposed into extrinsic and intrinsic components. On the one hand, we model the extrinsic noise as the variability of parameters including the maximum production rate $v$ and the degradation rate $r$ (Figure 1E; see ‘Methods’): each of these parameters is added by a noise term with zero mean, and the standard deviation of the noise term is proportional to the value of the kinetic parameter. On the other hand, the intrinsic noise, generated by the stochasticity of discrete chemical reactions, is modeled by directly simulating the dynamics of molecular numbers rather than concentrations. To this end, we introduce the cell volume $V$ , and naturally the molecular number of each node is the product of the cell volume $V$ and the concentration. As reactions progress, the molecular numbers would randomly increase or decrease by one at some time point (Figure 1E; see ‘Methods’), and the waiting time of the increase and decrease obeys exponential distributions with parameters determined by the production and decay rates in the deterministic model, respectively. This stochastic process can be exactly solved by the Gillespie algorithm, which has been widely used in previous studies (Liu et al., 2020; Thattai and van Oudenaarden, 2001; Veliz-Cuba et al., 2015; Zhao et al., 2021); however, the computation cost is high, and thus we use chemical Langevin equations as approximations to reduce the cost (Gillespie, 2000). Although the biological noise in nature usually has the extrinsic and intrinsic components simultaneously, we only consider the case where only one source of noise exists for simplicity, that is, only extrinsic noise exists or only intrinsic noise exists.

Procedures to search for network topologies robustly executing accurate oscillation

To search for two- and three-node network topologies that can robustly achieve accurate oscillation (i.e., high dimensionless correlation time $τ / T$ ), two steps are performed (Figure 1F): the first step is to identify topologies capable of oscillation in the whole network topology space (the upper panel in Figure 1F); the second step is to use the 90-percentile value of $τ / T$ to quantify the robustness of each oscillatory network topology to achieve accurate oscillation (the lower panel in Figure 1F). For a given topology, the 90-percentile value of $τ / T$ is defined as the value of $τ / T$ below which 90% of $τ / T$ ’s fall when 1000 parameter sets are randomly assigned. We refer the reader to ‘Methods’ for details, and here we only show major procedures. In the first step (the upper panel in Figure 1F), to obtain oscillatory network topologies in the whole network topology space, we randomly assign 10,000 parameter sets for each network topology and simulate the deterministic dynamics. The oscillatory network topology is chosen by the following two criteria: the network topology without repressilator is regarded as an oscillatory network topology if at least 80 parameter sets are capable of oscillation; the network topology with repressilator is defined as an oscillatory network topology if at least 10 parameter sets achieve oscillation. In this way, we finally obtain 474 oscillatory network topologies, and nearly 35% of them are with the repressilator. If we used the threshold of 80 for all network topologies, oscillatory network topologies with repressilator only occupy 20% of all oscillatory network topologies, which may lose the generality of conclusions about the repressilator. In the second step (the lower panel in Figure 1F), for each of these 474 oscillatory network topologies, we sample 1000 parameter sets capable of oscillation in the absence of noise and calculate the 90-percentile value of $τ / T$ in the presence of extrinsic noise or intrinsic noise. This value measures the robustness of the given topology against noise: the higher the value is, the larger probability to achieve accurate oscillation the topology has.

The robustness of accurate oscillation against extrinsic noise for different network topologies

Classification of all 474 oscillatory network topologies

We start by classifying all 474 oscillatory network topologies according to five types of core motifs. These five types of core motifs are as follows: the first core motif (shown in brown in Figure 2A) is composed of the repressilator and a positive autoregulation, but the node with the positive autoregulation is not allowed to have a positive incoming link; the second core motif (shown in orange in Figure 2A) is similar to the first core motif except that the positive incoming link to the positive autoregulated node is required; the third type of core motifs include the activator-inhibitor topology and its two variants (shown in green in Figure 2A); the fourth and fifth core motifs are the repressilator and delayed negative feedback (Figure 2A), respectively. Based on the identification of five types of core motifs, we define C1 category as the network topologies that contain only the first type of core motif, and so do the C2 category, C3 category, C4 category, and C5 category. The above five categories constitute near 59% of all 474 oscillatory network topologies, while the rest of topologies are those containing at least two of these five types of core motifs. Note that these oscillatory network topologies all have a negative feedback structure, which is consistent with previous studies (Glass and Pasternack, 1978; Novák and Tyson, 2008).

Figure 2 with 2 supplements see all

Download asset Open asset

The relationship between network topology and robustness to extrinsic noise.

(A) Venn diagram of all 474 oscillatory network topologies. The C1, C2, C3, C4, and C5 categories are nonoverlapping collections of network topologies that contain core motifs in brown, orange, green, blue, and dark blue, respectively. The number in red is the number of topologies in each region. (B) The violin plots of 90 percentiles of dimensionless correlation time ( $τ / T$ ) for C1, C2, …, and C5 categories present in (A). Each category corresponds to one violin plot, and topologies are denoted by dots. The Wilcoxon rank-sum tests (one-tailed) are applied to adjacent categories (***p<0.001; ****p<0.0001). (C) The violin plots of 90 percentiles of $τ / T$ for topologies that do not belong to any of C1, C2, …, and C5 categories. Taken C13 category as an example, the C13 category is the collection of the topology containing the first and third types of motifs simultaneously. (D) Illustration of calculations of the period average (denoted as $T_{o p t}$ ) and the amplitude average (denoted as $A_{o p t}$ ) for a given topology. For a given topology, $T_{o p t}$ (or $A_{o p t}$ ) is the expectation of the best-fit exponential distribution of 1000 periods (or 1000 amplitudes) with random assigned parameters. (E) The scatter plots of 90 percentiles of $τ / T$ versus the ranks according to $T_{o p t}$ (left) and $A_{o p t}$ (right). Each oscillatory network topology is denoted by a dot, with color determined by its category. Spearman coefficients are illustrated on the top of each plot, and the plot with a relatively high Spearman coefficient is highlighted by the black box. (F) The violin plots of $T_{o p t}$ for the five network topology categories present in (A). The Wilcoxon rank-sum tests (one-tailed) are applied to adjacent categories (ns, not significant; ***p<0.001; ****p<0.0001). (G) Hypothesis of topology–period–robustness relation. The network topology may affect robustness to extrinsic noise by acting on an intermediate quantity – period. (H) Typical dynamics in the presence of extrinsic noise and corresponding autocorrelation functions for the five topologies present in (A). From top to bottom, the period length decreases, and the autocorrelation also decreases. The amplitude for those dynamics in the absence of noise is almost the same (near 12). The 90 percentile of $τ / T$ is the average from five replicates.

The topologies containing repressilator with positive autoregulation perform better than those containing the activator-inhibitor topology when facing extrinsic noise, and the positive autoregulation enhances the robustness against extrinsic noise

Next we compare the robustness of accurate oscillation against extrinsic noise among above C1, C2, …, and C5 categories. Figure 2B shows the violin plots of 90 percentiles of $τ / T$ in the presence of extrinsic noise for the five categories, where each violin corresponds to one category. By applying one-tailed Wilcoxon rank-sum tests to adjacent two categories, we find that 90 percentiles of $τ / T$ for C1 category are significantly larger than those for C2 category, and this relation also holds between C2 and C3 categories, between C3 and C4 categories, and between C4 and C5 categories. These findings indicate that the order of these five categories according to the robustness of accurate oscillation to extrinsic noise is C1 > C2 > C3 > C4 > C5. The facts that C1 > C3 and C2 > C3 demonstrate that topologies containing the repressilator with positive autoregulation achieve higher robustness against extrinsic noise than those containing the activator-inhibitor topology. Besides, core motifs in both C1 and C2 categories have an extra positive autoregulation in comparison to the core motif (the repressilator) in C4 category, suggesting that the higher robustness of C1 and C2 categories than C4 category against extrinsic noise is due to the effect of positive autoregulation in improving the robustness to extrinsic noise. This effect is also validated by the comparison of the robustness between C3 and C5 categories.

The above analyses focused on the oscillatory network topologies in C1, C2, …, and C5 categories, which account for nearly 59% of all oscillatory network topologies. To perform a complete research, we also investigate the robustness to extrinsic noise for the remaining 41% oscillatory network topologies. These topologies all contain at least two types of core motifs and can be classified into seven categories, based on what core motifs the topology has. Then we compare each of them with its ‘component’ category (i.e., C1, C2, …, or C5 categories). For example, C13 category is composed of topologies that contain both the first and third types of core motifs, and its two ‘component’ categories are defined as C1 and C3 categories. The comparison is made in Figure 2C, where each group of violin plots separated by dashed lines represents the 90 percentiles of $τ / T$ (in the presence of extrinsic noise) for the category with combined core motifs and its ‘component’ categories. It can be seen that the category with combined core motifs usually shows intermediate robustness among its ‘component’ categories. That is to say, if a network topology has low robustness against extrinsic noise, adding a high-robustness core motif usually improves the robustness, but the combined topology cannot outperform the added high-robustness core motif.

Topologies with long period achieve high robustness against extrinsic noise

The above analyses indicate that network topologies differ widely in their robustness to achieve accurate oscillation, then we ask what mechanisms cause these differences. 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 birth-and-death Markov process, Paulsson obtained the variance of the protein in gene networks and found that it is related to the network’s elasticity, which is calculated from the deterministic model. Based on these observations, we explore two important characteristics for the oscillator: period and amplitude. Instead of focusing on a specific oscillatory network topology, we consider all 474 oscillatory network topologies and study what period and amplitude each topology prefers. To be precise, for each network topology, we calculate the distributions of period and amplitude from 1000 randomly sampled oscillation parameter sets and approximate mean values of period and amplitude by $T_{o p t}$ and $A_{o p t}$ , respectively. Here, we refer to the amplitude as the maximal peak value among nodes A–C; $T_{o p t}$ (or $A_{o p t}$ ) is defined as the expectation of the best-fit exponential distribution of 1000 periods (or 1000 amplitudes) (Figure 2D). Therefore, the topology with large $T_{o p t}$ tends to oscillate with long period, and the topology with large $A_{o p t}$ usually indicates an oscillation with high amplitude. Note that these two quantities are calculated in the noise-free system, and thus are not affected by the amplitude of the noise source or the type of noise.

To investigate the role of above two quantities $T_{o p t}$ and $A_{o p t}$ in the robustness of accurate oscillation against extrinsic noise, we calculate Spearman coefficients between these two quantities and 90 percentiles of $τ / T$ for all 474 oscillatory network topologies (Figure 2E). In Figure 2E, each dot represents an oscillatory network topology, with the x-axis representing the ranking according to $T_{o p t}$ (left panel) or $A_{o p t}$ (right panel). The Spearman coefficient between $T_{o p t}$ and 90 percentile of $τ / T$ for all 474 oscillatory network topologies is 0.94, which is larger than that between $A_{o p t}$ and 90 percentile of $τ / T$ (0.88). This result not only holds for all 474 oscillatory network topologies, but also holds within each of C1, C2, …, and C5 categories (Figure 2—figure supplement 1). These findings indicate that the robustness to extrinsic noise is more highly correlated with long period rather than high amplitude.

Since the long period benefits the robustness to extrinsic noise, then we ask how network topologies affect the period and whether those topologies with long period indeed show high robustness to extrinsic noise. To answer these questions, we analyze $T_{o p t}$ for C1, C2, …, and C5 categories (Figure 2F). The ranking of these five categories according to $T_{o p t}$ is C1 > C2 ≈ C3 > C4 > C5, which is obtained by the one-tailed Wilcoxon rank-sum tests for each adjacent two categories. This ranking is almost the same as that according to the robustness of accurate oscillation to extrinsic noise (C1 > C2 > C3 > C4 > C5) except rankings for C2 and C3 categories, suggesting that the topology with long period usually leads to high robustness of accurate oscillation to extrinsic noise (Figure 2G). The only inconsistency is that C2 and C3 categories differ in the robustness but have no significant difference in the probability to achieve long period. That is to say, C2 category might show better robustness to extrinsic noise than C3 category though they have the same period. Figure 2H shows typical dynamics for five different topologies when extrinsic noise exists. Those topologies from the top panel to the bottom panel belong to categories C1 to C5, respectively. Their dynamics have almost the same amplitude, but the period, as well as the autocorrelation, decreases when categories vary from C1 to C5. These findings suggest that topologies with prolonged period tend to have good performance to filter extrinsic noise, and this correlation is less likely due to that they have different amplitude.

The robustness of accurate oscillation against intrinsic noise for different network topologies

In the presence of only intrinsic noise, the repressilator with positive autoregulation is still better than the activator-inhibitor, and the advantage of positive autoregulation still holds

Unlike considering the robustness of accurate oscillation against parameter variability in the previous section, we next study the case where only intrinsic noise exists. With the same oscillatory network topology categories present in Figure 2A, 90 percentiles of the dimensionless correlation time ( $τ / T$ ) in the presence of only intrinsic noise also show a roughly similar trend from C1 to C5 categories (Figure 3A) except that C2 category exhibits the same robustness with C1 category to intrinsic noise while C2 category shows lower robustness than C1 in the presence of extrinsic noise. Moreover, the higher robustness of C1 and C2 categories compared with C3 category validates the better performance of the repressilator with positive autoregulation than the activator-inhibitor topology; the improvement of robustness from C4 category to C1 (or C2) category indicates the effect of positive autoregulation on the robustness to intrinsic noise; the comparison between the robustness of C5 and C3 categories also implies the advantage of the positive autoregulation. These findings are consistent with those when noise is only originated from the extrinsic noise. Another consistency is that the hybrid of core motifs imparts an intermediate robustness not only in the presence of the extrinsic noise (Figure 2C) but also in the presence of the intrinsic noise (Figure 3B).

Figure 3 with 4 supplements see all

Download asset Open asset

The relationship between network topology and robustness to intrinsic noise.

(**A–C**) The same plots as in Figure 2B, C and E, except that the noise comes from the stochasticity of biochemical reactions. The 90 percentile of dimensionless correlation time ( $τ / T$ ) shows a higher correlation with $A_{o p t}$ rather than $T_{o p t}$ . (D) The violin plots of $A_{o p t}$ for the five network topology categories present in Figure 2A. The Wilcoxon rank-sum tests (one-tailed) are applied to adjacent categories (ns, not significant; ***p<0.001; ****p<0.0001). (E) Hypothesis of topology–amplitude–robustness relation. The correlation between $A_{o p t}$ and 90 percentiles of $τ / T$ in the presence of intrinsic noise indicates that the network topology may influence the robustness against intrinsic noise through the amplitude. (F) The same plot as Figure 2H except in the presence of intrinsic noise. From top to bottom, while the period in the absence of noise is very close (near 12), the amplitude decreases, and the autocorrelation also decreases. The 90 percentile of $τ / T$ is the average from five replicates.

Topologies with the high amplitude enable high robustness against intrinsic noise

Similar to the analysis of robustness to extrinsic noise, we try to answer whether the period or amplitude is highly correlated with the robustness to intrinsic noise. In the presence of only intrinsic noise, the Spearman correlation coefficient of 90 percentiles of $τ / T$ and $T_{o p t}$ for all 474 oscillatory network topologies is 0.72, which is smaller than that of 90 percentiles of $τ / T$ and $A_{o p t}$ (0.81) (Figure 3C). These findings suggest that unlike the case of extrinsic noise where the robustness is more strongly correlated with period, the robustness of accurate oscillation against intrinsic noise is more highly correlated with amplitude. In other words, the topology with the high amplitude has a larger probability to achieve high robustness against intrinsic noise than that with long period. However, it should be noted that the correlation coefficient between 90 percentiles of $τ / T$ and $A_{o p t}$ for all oscillatory network topologies is not very close to 1 (the right panel in Figure 3C), and it is also much smaller than 1 even in each category (Figure 3—figure supplement 1), implying that the relation between the amplitude and robustness to intrinsic noise is not very strong, and that some topologies with small amplitude may perform better than those with high amplitude. Therefore, there might exist other mechanisms to attenuate intrinsic noise.

Furthermore, by applying Wilcoxon rank-sum tests to amplitude average ( $A_{o p t}$ ) for neighboring two categories, we find that the ranking of five network topology categories according to amplitude is C1 > C2≈ C3 > C4 > C5 (Figure 3D). This ranking is almost the same as that according to the robustness to intrinsic noise (C1 ≈ C2 > C3 > C4 > C5) (Figure 3E), implying that the amplitude might link the topology category and the robustness to intrinsic noise. The only exception is C2 category: because of the fact that C1 > C2 ≈ C3 according to $A_{o p t}$ and the fact that the amplitude strongly correlates with the robustness, the C2 category is supposed to show the same robustness to intrinsic noise as C3 category and exhibit lower robustness than C1 category; however, the robustness to intrinsic noise for C2 category is actually at the same level of C1 category, further demonstrating that the high amplitude is not the only mechanism to enhance the robustness to intrinsic noise (Figure 3—figure supplement 1). Figure 3F shows typical dynamics in the presence of intrinsic noise whose topologies belong to distinct categories. Those dynamics exhibit near period, but their amplitudes and autocorrelations decrease from category C1 to category C5, which supplies a possibility to enhance the robustness to intrinsic noise through varying topologies while maintaining period.

Simulations using the Gillespie algorithm lead to similar conclusions

The above analyses are based on simulations for chemical Langevin equations, which can only give approximate solutions of the dynamical behavior in the presence of intrinsic noise. To test whether this approximation is feasible, we use the Gillespie algorithm to exactly solve the stochastic dynamical behavior when facing intrinsic noise, and then conduct similar analyses (Figure 3—figure supplement 2) as the previous section has done. According to the robustness rankings for C1–C5 categories, the repressilator with positive autoregulation performs better than the activator-inhibitor, and the topologies with positive autoregulation are better than that without positive autoregulation. Besides, the robustness is more correlated with the mean amplitude rather than the mean period, and the order of the five categories sorted by the mean amplitude is almost the same as that sorted by robustness, indicating the bridge role of amplitude to link topologies and the robustness to intrinsic noise. These results are consistent with the conclusions based on chemical Langevin equations. We also find that the Gillespie algorithm leads to higher dimensionless correlation times than chemical Langevin equations since the maximal correlation in Figure 3A is near 6 and that in Figure 3—figure supplement 2A is 40. However, this difference does not indicate that chemical Langevin equations are bad approximations: when the system behaves normal noise filtering capability, these two methods give similar dimensionless correlation times (Figure 3—figure supplement 3A); when the system buffer noise perfectly, dimensionless correlation times calculated through these two methods differ a lot, but autocorrelation functions remain similar (Figure 3—figure supplement 3B), which indicates that chemical Langevin equations still capture the system’s ability to buffer noise. The reason why large and extremely large dimensionless correlation times result in almost same correlations might be that doubling long correlation time cannot increase autocorrelation efficiently due to the property of the exponential function.

Relations between period/amplitude and oscillation accuracy against noise are validated by analytical approaches

The above simulations revealed relations between two important features (i.e., period and amplitude) of the oscillator and the oscillator’s robustness to noise. However, these results just showed the correlation rather than the causal relationship. Besides, because the period and amplitude are usually positively correlated (Figure 4—figure supplement 1), it is hard to control one feature and analyze the effect of the other feature. Fortunately, these two problems can be solved by introducing the timescale or rescaling parameters. In this way, we can change one feature while maintaining the other feature, and then analytically derive causal relations between period or amplitude and the oscillation accuracy. We will illustrate these methods and corresponding results below.

To study the relation between period and oscillation accuracy against noise, we maintain the amplitude and tune the period through changing the factor $M$ on the right-hand side in ordinary differential equations (Figure 4A), and then analyze the phase noise through the analytical approach proposed by Demir et al., 2000. Varying $M$ can be regarded as the rescaling of time $t$ , so the period is changed while maintaining the amplitude, and thus we can focus on the effect of period on the oscillation accuracy. In order to analyze the system with variable $M$ , we first summarize Demir et al.’s work. They carried out nonlinear perturbation analysis for oscillators and obtained an exact equation for phase deviation. We only summarize the main results below. The dynamics of a perturbed oscillator can be described as a set of differential equations:

\dot{x} = f (x) + B (x) ξ (t)

Figure 4 with 1 supplement see all

Download asset Open asset

Effects of tuning period or amplitude on oscillation accuracies.

(**A, B**) Relations between period (A) or amplitude (B) and oscillation accuracies. For the activator-inhibitor topology, we arbitrarily choose an oscillation parameter set (see Supplementary file 1b for parameters), and then we add $M$ or $N$ in the way shown in the upper panel to tune the period or amplitude, respectively. The lower panels are simulation results of the period, amplitude, and dimensionless correlation time when $M$ or $N$ increases.

where $x \in R^{3}$ , $f (\cdot) : R^{3} \to R^{3}$ , $B (\cdot) : R^{3} \to R^{3 \times 3}$ and $ξ (t) \in R^{3}$ is random noise. The unperturbed system $\dot{x} = f (x)$ has a periodic solution $x_{s} (t)$ (with period $T$ ). It can be proved that the variance of the phase deviation $σ^{2} (t)$ satisfies $σ^{2} (t) = c t$ , and $c$ is as follows:

\begin{matrix} c = \frac{1}{T} \int_{0}^{T} v_{1}^{T} (t^{`}) B (x_{s} (t^{`})) B^{T} (x_{s} (t^{`})) v_{1} (t^{`}) d t^{`}, \end{matrix}

where $v_{1}^{T} (t)$ is the first row of the matrix $V (t)$ . Here, the first column of $V^{- 1} (t)$ is ${\dot{x}}_{s} (t)$ , and $V^{- 1} (t) d i a g [μ_{1}, μ_{2}, μ_{3}] V (t)$ is the state transition matrix for $\dot{w} = A (t) w (t)$ where $μ_{i}$ ’s are Floquet exponents and $A (t) = {\frac{\partial f (x)}{\partial x} |}_{x = x_{s} (t)}$ (see ‘Methods’ for details). The Equation 01 gives an analytic expression describing the phase noise, so we use the dimensionless $c$ , that is, $c / T$ , to measure the oscillation accuracy instead of $τ / T$ . To calculate $c / T$ for the systems with $M$ , we use $T$ to denote the period for the system without $M$ , and then the period for this new system is $M T$ . Besides, $v_{1} (t)$ becomes $M v_{1} (\frac{t}{M})$ (see ‘Methods’ for details). For the noise term, we merge $\frac{1}{M}$ with kinetic parameters $v_{i j}, δ_{i}, r_{i}$ , that is, these parameters become $\frac{1}{M}$ of original values, and then we model the extrinsic and intrinsic noise as that in Figure 1E. Therefore, $B (x)$ becomes $B (x) \frac{1}{M}$ when facing extrinsic noise as the magnitude of noise source is proportional to the kinetic parameters. However, in the presence of only intrinsic noise, $B (x)$ becomes $B (x) \sqrt{\frac{1}{M}}$ because the noise term in the chemical Langevin equation is usually the square root of reaction rates. Then we can calculate the ratio of the slope of the variance of the phase noise to the period ( $c / T$ ) using the Equation 01 (see ‘Methods’ for details). We find that the $c / T$ in the presence of only extrinsic noise is proportional to the $1 / M$ , and that in the presence of only intrinsic noise is not affected by $M$ . Note that the smaller the $c / T$ is, the more accurate the oscillation is. Thus, large $M$ enhances the oscillation accuracy against extrinsic noise, which is also numerically validated by the trend of dimensionless correlation times for the system with different $M$ (right lower panel in Figure 4A). Since large $M$ leads to long period but has no influence on amplitude, the prolonged period might be the reason for high oscillation accuracy in the presence of extrinsic noise.

For the relation between amplitude and oscillation accuracy against noise, we keep the period and tune the amplitude through the rescaling parameter $N$ and then analyze the rescaled system. For a fixed topology with a set of oscillation kinetic parameters, we replace the variables $A$ , $B$ , and $C$ with $\tilde{A} / N$ , $\tilde{B} / N$ , and $\tilde{C} / N$ , respectively (Figure 4B). This rescaling makes amplitudes of $\tilde{A}$ $N$ times as high as that of $A$ , and so do $\tilde{B}$ and $\tilde{C}$ . However, this rescaling has no influence on the period, so we can focus on the role of amplitude in the oscillation accuracy. The system with rescaled variables $\tilde{A}$ , $\tilde{B}$ , and $\tilde{C}$ shows unchanged oscillation accuracy against extrinsic noise with varied $N$ , but the oscillation accuracy against intrinsic noise increases with increased $N$ (see ‘Methods’). Taken together, large $N$ not only increases the amplitude but also improves the oscillation accuracy to intrinsic noise while maintaining the period. These results are consistent with numerical simulations for tendencies of period, amplitude, and dimensionless correlation times (lower panel in Figure 4B). These results indicate that the improvement of the oscillation accuracy to intrinsic noise may due to the high amplitude rather than period.

Analyses of synthetic NF-κB signaling circuits demonstrate the improvement of the oscillation accuracy when adding a repressilator topology to the activator-inhibitor

In previous sections, we have used two- and three-node networks to approximate biological systems and focused on the noise coming from the variability in kinetic parameters or chemical reactions. Though biological systems are more complex than two- or three-node networks and face noise from various sources besides the above noises, the investigation of a specific biological system— a synthetic NF-κB signaling circuit is consistent with the theoretical results in previous sections. As described in our previous work, we implement the design of negative feedback-only circuit 1 (Figure 5A) by integrating the synthetic RelA-IκBα signaling circuit into the yeast MAPK pathway. The nuclear-to-cytoplasmic RelA oscillations can be triggered by inducing the degradation of IκBα through the activation of yeast MAP kinase Fus3, and we can monitor these single-cell oscillations for up to 10 hr. Base on this simple circuit, we then modify its structure by adding extra regulations. One modification is adding constantly expressed IκBα. This copy of IκBα also inhibits RelA and is inhibited by Fus3, so it provides another pathway from Fus3 to activate RelA (the orange link in circuit 2 in Figure 5A). Another modification is adding a yeast MAPK phosphatase Msg5 (the orange link in circuit 3 in Figure 5A), which is activated by RelA and can dephosphorylate Fus3. In circuit 3, Msg5, RelA, and IκBα form a repressilator topology. To obtain the single-cell time trajectories for these three circuits, we employ time-lapse microscopy to track the RelA nuclear localization dynamics for over 10 hr. The period lengths are determined as the time intervals between the successive peaks of these trajectories, and then we calculate the CV of those period lengths for over 50 cells. We find that the circuit 2 shows similar CV of period as that for the circuit 1, but that circuit 3 exhibits lower CV of period than circuit 1 (Figure 5B). These results suggest that the additional repressilator topology indeed facilitates the noise buffering capability for the activator-inhibitor topology.

Figure 5

Download asset Open asset

Experimental evidence of the improvement of oscillation accuracy when adding a repressilator topology to the activator-inhibitor.

(A) Schematic showing the synthetic NF-κB signaling circuits (circuit 1) and other two circuits (circuits 2 and 3) with more complicated structures. In circuit 1, RelA activates IκBα, but the latter inhibits RelA, constituting a negative feedback; we regard circuit 1 as the activator-inhibitor where the positive autoregulation is replaced by multiple reactions. Based on circuit 1, we construct circuit 2 by adding a copy of IκBα to form a new link from Fus3 to NF-κB (RelA), where RelA, IκBα, and Fus3 cannot constitute a repressilator. Circuit 3 is circuit 1 with additional negative regulation from RelA to Fus3, which is constructed by adding a yeast MAPK phosphatase Msg5 that is activated by RelA and dephosphorylate Fus3. Circuit 3 includes a repressilator consisting of RelA, IκBα, and Fus3. (B) Bar plots of coefficients of variation (CVs) of period for three synthetic circuits present in (A). For each circuit, time required from the peak of RelA to next peak from 0 to 10 hr for over 50 cells is recorded to calculate the CV. Two replicates are performed for each circuit. (C) Procedures to study how adding repressilator topology to activator-inhibitor affect the oscillation accuracy by using mathematical models present in Figure 1. $K_{C A}$ represents the binding affinity of protein C to gene A, and large $1 / K_{C A}$ indicates strong inhibition from protein C to protein A. (**D–G**) Simulation results for period (D), amplitude (E), and dimensionless correlation time in the presence of extrinsic (F) or intrinsic (G) noise when $1 / K_{C A}$ increases. See Supplementary file 1c for parameters. The marker with x coordinate 0 denotes the value for activator-inhibitor. The error bar represents the standard deviation for 100 repeated simulations. Increasing $1 / K_{C A}$ (i.e., strengthening the inhibition from protein C to protein A) prolongs period, increases the amplitude, and improves the oscillation accuracy.

Such improvement of the oscillation accuracy when adding a repressilator topology to the activator-inhibitor in the synthetic NF-κB circuit can also be validated using the mathematical models as described in Figure 1. We use nodes A, B, and C to denote IκBα, RelA, and Fus3, respectively, and thus circuits 1 and 3 in Figure 5A are networks shown in upper and lower panels in Figure 5C, respectively. These two networks are interconvertible by tuning $K_{C A}$ (the binding affinity of protein C to gene A): (near) zero $1 / K_{C A}$ indicates the little effect of protein C to protein A, and thus protein A is not affected by protein C, leading to the activator-inhibitor; non-zero $1 / K_{C A}$ implies the existence of the inhibition from protein C to protein A, resulting in the network with an activator-inhibitor and a repressilator. For a given oscillation parameter set for the activator-inhibitor, we first set $1 / K_{C A}$ to be (near) zero and calculate period, amplitude, and dimensionless correlation time in the presence of extrinsic or intrinsic noise (the first points in Figure 5D–G); then we increase $1 / K_{C A}$ , and calculate the corresponding quantities (i.e., period, amplitude, and dimensionless correlation time, as shown from the second points in Figure 5D–G). We can find that the period, amplitude, and dimensionless correlation time for the activator-inhibitor are usually smaller than those with an additional repressilator, and this gap enlarges with increased $1 / K_{C A}$ , that is, the increased strength of the negative regulation from the additional node C to the inhibitor node A (Figure 5C). Therefore, we demonstrate that adding the repressilator to the activator-inhibitor enhances the oscillation accuracy. This is consistent with that C1 and C2 categories exhibit higher robustness than C3 categories. Moreover, the prolonged period and increased amplitude, which are also observed in C1 and C2 categories, may be the reason for such enhancement (Figure 5D and E).

Discussion

It remains the major challenge in biology to understand how living systems perform complex behaviors accurately in the presence of inevitable noise. Instead of studying biological networks case by case, we try to answer whether there exist general network design principles for living systems to execute biological functions by using a bottom-up strategy (Lim et al., 2013; Zhang and Tang, 2019).

Here, we systematically explore the network design principles for accurate oscillation in both two- and three-node networks. We identify several key motifs that have distinct robustness to noise. The motif —— a repressilator with positive autoregulation behaves better than other motifs present in Figure 2A in most cases, especially the activator-inhibitor oscillator; the additional positive autoregulation can improve the robustness. These results are consistent in spite of sources of noise. However, different sources of noise utilize distinct mechanisms to filter noise: the variability of parameters, a type of extrinsic noise, is largely filtered through long period, and the intrinsic noise is buffered by high amplitude.

Interestingly, investigations of three engineered NF-κB signaling circuits partly validate our simulation results. For the negative-feedback loop circuit, if the additional new regulations form a repressilator, low variance of period will occur, but if no repressilator emerges, the variance of period shows no significant change. These findings show the advantage of the repressilator against noise.

While modifying network topology and changing regulation strength for a fixed topology are both options to improve the robustness of accurate oscillation, each network’s robustness is an indicator of the probability of this network topology achieving accurate oscillation with varied regulatory strengths (Figure 2—figure supplement 2, Figure 3—figure supplement 4): the network topology with high robustness tends to show high dimensionless autocorrelation time when varying regulatory strengths, that is, accurate oscillation (first 10 bars in Figure 2—figure supplement 2 and Figure 3—figure supplement 4); the network topology with low robustness displays a bad performance of oscillation accuracy in the whole parameter space (last 10 bars in Figure 2—figure supplement 2 and Figure 3—figure supplement 4). Besides, our work also suggests that tuning network topology is more efficient than changing regulatory strength. This is based on the observations that network topologies with low robustness (last 10 bars in Figure 2—figure supplement 2 and Figure 3—figure supplement 4) cannot have a high oscillation accuracy even when searching all kinetic parameter space, but changing topologies may increase the probability of high oscillation accuracy. So we suggest that a feasible way to improve the oscillation accuracy in synthetic networks is to first modify the topology to avoid low-robustness ones and then tune the regulation strength, as illustrated in Figure 5C.

Mechanisms to buffer different sources of noise in the oscillator can be dramatically different. On the one hand, long period is able to attenuate extrinsic noise, which is also called the time-averaging strategy. This strategy has been widely studied in nonoscillatory networks, such as circuits that are sensitive to the stimulus (Hornung and Barkai, 2008), circuits with ‘‘switch-like’’ behaviors (Wang et al., 2010), and adaptive circuits (Nie et al., 2020; Sartori and Tu, 2011). For these nonoscillatory circuits, fluctuations in output have been proven to be related to some key timescales, and long timescales often result in the output with small variance. On the other hand, the intrinsic noise is hard to be attenuated through time averaging, such as the adaptive incoherent feed-forward loop (Sartori and Tu, 2011). Actually, the right way to buffer intrinsic noise in biological oscillators was found to depend on levels of molecules. For example, the importance of protein numbers has also been demonstrated in the work of Potvin-Trottier et al., 2016. They found that increasing the peak and bottom values decreases the CV in the decay phase of the oscillator. Based on these results, it is suggestive that the network topology with long period and high amplitude may enable good robustness to both extrinsic and intrinsic noise. Interestingly, it is usually not hard to obtain long period and high amplitude simultaneously since the long period tends to allow the protein number to climb to a high level (Figure 4—figure supplement 1).

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 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 toggle-switch system and excitable system (Lindner et al., 2004; Terebus et al., 2019; Zaks et al., 2005). The possible reason might be the noise-induced 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.

In this work, the extrinsic noise is assumed only from fluctuations in kinetic parameters, and its magnitude linearly depends on the level of the parameter. Except this type of extrinsic noise, cells also face the random partitioning that occurs during cell division, noisy stimulus, and so on (Monti et al., 2018; Veliz-Cuba et al., 2015). Since two different types of noise studied in this work require different mechanisms to buffer, other sources of noise may also need new mechanisms to filter. Thus, some unknown principles need to be further revealed and incorporated into network design as the increasingly improved complexity and multiple sources of noise.

Another limitation of our work is that we did not decompose the reactions in the deterministic model into detailed elementary reaction steps when using the Gillespie algorithm. The advantage of simulating nonelementary reactions with Hill-type 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., 2002b; 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 Hill-type reaction rate is based on the quasi-steady-state approximation, which does not hold when binding/unbinding of TF to the promoter is slow or comparable to the timescales of protein production or decay (Choi et al., 2008). Furthermore, this method neglects detailed reaction in gene regulatory networks, and thus fails to study the 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, 2008; Terebus et al., 2019). We anticipate the need for a more detailed model where every reaction of Hill-type 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 three-node 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 the 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.

Methods

Mathematical modeling

Deterministic model

To model two-node and three-node network topologies, we use transcriptional regulatory networks and assume competitive inhibition among multiple transcription factors. The competitive inhibition means that multiple transcription factors compete for the same binding sites if they regulate one gene simultaneously (Shi et al., 2017). So the gene expression depends on the relative weight of transcription factors inhibiting this gene and that activating this gene. The following set of ordinary equations is used to describe the deterministic dynamics of a three-node transcriptional regulatory network:

{\begin{cases} \frac{d A}{d t} = k_{b a s a l} + \frac{\sum_{i} v_{i A} {(\frac{x_{i}}{K_{i A}})}^{3} + δ_{A}}{1 + \sum_{i} {(\frac{x_{i}}{K_{i A}})}^{3} + \sum_{i} {(\frac{y_{i}}{K_{i A}})}^{3}} - r_{A} A \\ \frac{d B}{d t} = k_{b a s a l} + \frac{\sum_{i} v_{i B} {(\frac{x_{i}}{K_{i B}})}^{3} + δ_{B}}{1 + \sum_{i} {(\frac{x_{i}}{K_{i B}})}^{3} + \sum_{i} {(\frac{y_{i}}{K_{i B}})}^{3}} - r_{B} B \\ \frac{d C}{d t} = k_{b a s a l} + \frac{\sum_{i} v_{i C} {(\frac{x_{i}}{K_{i C}})}^{3} + δ_{C}}{1 + \sum_{i} {(\frac{x_{i}}{K_{i C}})}^{3} + \sum_{i} {(\frac{y_{i}}{K_{i C}})}^{3}} - r_{C} C \end{cases}

where $A$ , $B$ , and $C$ are the concentrations of proteins A, B, and C. $x_{i} = A, B$ or $C$ in each equation denotes the concentration of protein activating the gene, and $y_{i} = A, B$ or $C$ the concentration of protein inhibiting the gene. The production rate constant $v_{i j}$ represents the maximal production rate of protein $i$ regulated by protein $j$ , with $K_{i j}$ binding affinity. If there exist proteins activating gene i, $δ_{i}$ is zero; if no protein activates gene i, $δ_{i}$ is non-zero and represents the production rate caused by other proteins. $k_{b a s a l}$ is the basal production rate. The equations for the two-node transcriptional regulatory network can be obtained in a similar way.

To provide a better explanation about the nonlinear reaction term in above equations, we took the following case as an example: protein A (i.e., TF) binds to gene B to inhibit the gene expression, and protein B binds to the same site in gene B to activate the gene expression. We assumed that (1) there are three binding sites in gene B, which once protein A (or B) binds to, then B (or A) cannot. The elementary reactions are described as follows:

A + G_{B} ⇌ G_{B} ∙ A, A + G_{B} ∙ A ⇌ G_{B} ∙ A_{2}, A + G_{B} ∙ A_{2} ⇌ G_{B} ∙ A_{3},

B + G_{B} ⇌ G_{B} ∙ B, B + G_{B} ∙ B ⇌ G_{B} ∙ B_{2}, B + G_{B} ∙ B_{2} ⇌ G_{B} ∙ B_{3},

where $G_{B}, A, B$ denote gene B, protein A, and protein B, respectively. The dissociation rates ( $k_{r e v e r s e} / k_{f o r w a r d}$ ) for these six reactions are $K_{1}, K_{2}, \dots, K_{6}$ . Therefore, the fraction of the gene B at the state $G_{B} ∙ B_{3}$ in equilibrium is given by

\frac{\frac{1}{K_{4} K_{5} K_{6}} B^{3}}{1 + \frac{1}{K_{1}} A + \frac{1}{K_{1} K_{2}} A^{2} + \frac{1}{K_{1} K_{2} K_{3}} A^{3} + \frac{1}{K_{4}} B + \frac{1}{K_{4} K_{5}} B^{2} + \frac{1}{K_{4} K_{5} K_{6}} B^{3}} .

Furthermore, we assumed that (2) $K_{6} ≪ K_{5}, K_{4}$ and $K_{3} ≪ K_{2}, K_{1}$ , so that this fraction can be rewritten as

\frac{\frac{1}{K_{4} K_{5} K_{6}} B^{3}}{1 + \frac{1}{K_{1} K_{2} K_{3}} A^{3} + \frac{1}{K_{4} K_{5} K_{6}} B^{3}} .

We further assumed (3) that only gene B staying at the state $G_{B} ∙ B_{3}$ can lead to transcription and subsequent translation for protein B, and (4) that the binding/unbinding of TFs to a gene can achieve a rapid equilibrium as TF levels change, and thus the production rate of protein B is modeled as

v_{B B} \frac{\frac{1}{K_{4} K_{5} K_{6}} B^{3}}{1 + \frac{1}{K_{1} K_{2} K_{3}} A^{3} + \frac{1}{K_{4} K_{5} K_{6}} B^{3}},

where $v_{B B}$ is the maximal production rate when gene B is bound with three protein B. This form is the same as those in (1) if $K_{1} K_{2} K_{3}$ and $K_{4} K_{5} K_{6}$ are substituted by $K_{A B}^{3}$ and $K_{B B}^{3}$ , respectively.

Stochastic model in the presence of extrinsic noise

To generate extrinsic noise, we perturb each kinetic parameter $v_{i j}, δ_{i}, r_{i}$ by multiplying the sum of 1 and an independent temporal noise term, and obtain a new system described by the following stochastic differential equations:

{\begin{cases} \frac{d A}{d t} = k_{b a s a l} + \frac{\sum_{i} v_{i A} (1 + ε η_{i A}) {(\frac{x_{i}}{K_{i A}})}^{3} + δ_{A} (1 + ε ξ_{A})}{1 + \sum_{i} {(\frac{x_{i}}{K_{i A}})}^{3} + \sum_{i} {(\frac{y_{i}}{K_{i A}})}^{3}} - r_{A} A (1 + ε ζ_{A}) \\ \frac{d B}{d t} = k_{b a s a l} + \frac{\sum_{i} v_{i B} (1 + ε η_{i B}) {(\frac{x_{i}}{K_{i B}})}^{3} + δ_{B} (1 + ε ξ_{B})}{1 + \sum_{i} {(\frac{x_{i}}{K_{i B}})}^{3} + \sum_{i} {(\frac{y_{i}}{K_{i B}})}^{3}} - r_{B} B (1 + ε ζ_{B}) \\ \frac{d C}{d t} = k_{b a s a l} + \frac{\sum_{i} v_{i C} (1 + ε η_{i C}) {(\frac{x_{i}}{K_{i C}})}^{3} + δ_{C} (1 + ε ξ_{C})}{1 + \sum_{i} {(\frac{x_{i}}{K_{i C}})}^{3} + \sum_{i} {(\frac{y_{i}}{K_{i C}})}^{3}} - r_{C} C (1 + ε ζ_{C}) \end{cases}

Here, the control parameter $ε$ indicates magnitude of perturbation of kinetic parameters, and large $ε$ represents big perturbation of kinetic parameters. $η_{i j}$ , $ξ_{i}$ , $ζ_{i}$ are independent noise terms and all modeled by the Ornstein–Uhlenbeck process:

τ^{n o i s e} d z = - z d t + σ d W_{t}

where $W_{t}$ is standard Wiener processes. This equation implies that $z (t)$ has zero mean and variance $\frac{σ^{2}}{2 τ^{n o i s e}}$ .

Stochastic model in the presence of intrinsic noise

To induce intrinsic noise, we replace the concentration of protein with the number of protein by introducing the cell volume $V$ and assume that production events and degradation events occur independently and randomly. To be precise, we use $X_{A}, X_{B}, X_{C}$ to denote numbers of proteins A, B, C, respectively, and then we replace $A, B, C$ in Equation 2 by $X_{A} / V, X_{B} / V, X_{C} / V$ , respectively. Therefore, the dynamics of protein numbers $X_{A}, X_{B}, X_{C}$ are described by the following reactions:

X_{A} \to_{X_{A} \overset{r_{A} X_{A}}{\to} X_{A} - 1}^{V k_{b a s a l} + V (\sum_{i} v_{i A} {(\frac{X_{i}}{K_{i A}})}^{3} + δ_{A} V^{3}) / (V^{3} + \sum_{i} {(\frac{X_{i}}{K_{i A}})}^{3} + \sum_{i} {(\frac{Y_{i}}{K_{i A}})}^{3})} X_{A} + 1

X_{B} \to_{X_{B} \overset{r_{B} X_{B}}{\to} X_{B} - 1}^{V k_{b a s a l} + V (\sum_{i} v_{i B} {(\frac{X_{i}}{K_{i B}})}^{3} + δ_{B} V^{3}) / (V^{3} + \sum_{i} {(\frac{X_{i}}{K_{i B}})}^{3} + \sum_{i} {(\frac{Y_{i}}{K_{i B}})}^{3})} X_{B} + 1

X_{C} \to_{X_{C} \overset{r_{C} X_{C}}{\to} X_{C} - 1}^{V k_{b a s a l} + V (\sum_{i} v_{i C} {(\frac{X_{i}}{K_{i C}})}^{3} + δ_{C} V^{3}) / (V^{3} + \sum_{i} {(\frac{X_{i}}{K_{i C}})}^{3} + \sum_{i} {(\frac{Y_{i}}{K_{i C}})}^{3})} X_{C} + 1

We used the following two algorithms to simulate the above system.

Gillespie algorithm

We used the standard Gillespie algorithm to simulate the system. There are six reactions in total (as shown above), and the propensity functions are the reaction rates listed above the arrow. Note that we did not decompose the reactions with the Hill function rate into the elementary reactions; the reaction rate with the Hill function type has also been applied to other discrete stochastic models (Gonze and Goldbeter, 2006; Gonze et al., 2002b; Veliz-Cuba et al., 2015; Wang et al., 2019; Zhao et al., 2021) and proven to be an accurate approximation for the model composed of all elementary reactions under certain circumstances (Kim et al., 2014; Sanft et al., 2011). In our simulations, each of the six reactions occurs with the random waiting time, which obeys an exponential distribution with mean of the inverse of the propensity function. For example, the reaction of decreasing protein A by 1 has the propensity function r_AX_A, and increasing protein A by 1 corresponds to $k_{b a s e l} + V (Σ_{i} v_{i A} {(\frac{X_{i}}{K_{i A}})}^{3} + δ_{A} V^{3}) / (V^{3} + Σ_{i} {(\frac{X_{i}}{K_{i A}})}^{3} + Σ_{i} {(\frac{Y_{i}}{K_{i A}})}^{3})$ . Once we get one trajectory, we can calculate the autocorrelation time. See Figure 3—figure supplement 2 and Figure 3—figure supplement 3 for simulation results.

Langevin equations

We also used the Langevin equation, a good approximation of this system under certain conditions (Gillespie, 2000), to model the system. The corresponding chemical Langevin equations are as follows:

{\begin{cases} d X_{A} = (V k_{b a s a l} + V \frac{\sum_{i} v_{i A} {(\frac{X_{i}}{K_{i A}})}^{3} + δ_{A} V^{3}}{V^{3} + \sum_{i} {(\frac{X_{i}}{K_{i A}})}^{3} + \sum_{i} {(\frac{Y_{i}}{K_{i A}})}^{3}} - r_{A} X_{A}) d t + \sqrt{V k_{b a s a l} + V \frac{\sum_{i} v_{i A} {(\frac{X_{i}}{K_{i A}})}^{3} + δ_{A} V^{3}}{V^{3} + \sum_{i} {(\frac{X_{i}}{K_{i A}})}^{3} + \sum_{i} {(\frac{Y_{i}}{K_{i A}})}^{3}} + r_{A} X_{A}} d W_{t}^{A} \\ d X_{B} = (V k_{b a s a l} + V \frac{\sum_{i} v_{i B} {(\frac{X_{i}}{K_{i B}})}^{3} + δ_{B} V^{3}}{V^{3} + \sum_{i} {(\frac{X_{i}}{K_{i B}})}^{3} + \sum_{i} {(\frac{Y_{i}}{K_{i B}})}^{3}} - r_{B} X_{B}) d t + \sqrt{V k_{b a s a l} + V \frac{\sum_{i} v_{i B} {(\frac{X_{i}}{K_{i B}})}^{3} + δ_{B} V^{3}}{V^{3} + \sum_{i} {(\frac{X_{i}}{K_{i B}})}^{3} + \sum_{i} {(\frac{Y_{i}}{K_{i B}})}^{3}} + r_{B} X_{B}} d W_{t}^{B} \\ d X_{C} = (V k_{b a s a l} + V \frac{\sum_{i} v_{i C} {(\frac{X_{i}}{K_{i C}})}^{3} + δ_{C} V^{3}}{V^{3} + \sum_{i} {(\frac{X_{i}}{K_{i C}})}^{3} + \sum_{i} {(\frac{Y_{i}}{K_{i C}})}^{3}} - r_{C} X_{C}) d t + \sqrt{V k_{b a s a l} + V \frac{\sum_{i} v_{i C} {(\frac{X_{i}}{K_{i C}})}^{3} + δ_{C} V^{3}}{V^{3} + \sum_{i} {(\frac{X_{i}}{K_{i C}})}^{3} + \sum_{i} {(\frac{Y_{i}}{K_{i C}})}^{3}} + r_{C} X_{C}} d W_{t}^{C} \end{cases}

where $X_{i}$ is the number of protein $i$ , and $W_{t}^{i}$ is the standard Wiener process. The control parameter $V$ reflects the magnitude of stochasticity of biological reactions. The big $V$ indicates small degree of stochasticity of biological reactions. See the next section for settings in the numerical simulation.

Numerical simulations

Numerical simulations for deterministic models were carried out in MATLAB (see https://github.com/LingxiaQiao/oscillation, (copy archived at swh:1:rev:72a2d3d1146b14e7988c1cc06208fe1252e9a6f5; Qiao, 2022) for MATLAB scripts). We use the solver ode15s to simulate the dynamics. Simulations for stochastic models were also implemented in MATLAB. In the presence of extrinsic noise, we used the Milstein scheme (Kloeden and Platen, 2013) to numerically solve the noise term $η_{i j}, ξ_{j}, ζ_{j}$ and used the Euler scheme to solve the dynamics of proteins’ concentrations. To be specific, the noise term $z$ ( $z = η_{i j}, ξ_{j}, ζ_{j}$ ) at $n + 1$ time step is determined by the following manner ( $τ^{n o i s e} = 1$ ):

z^{(n + 1)} = z^{(n)} - z^{(n)} ∆ t + σ {δ W}_{n} + \frac{1}{2} σ^{2} [{({δ W}_{n})}^{2} - ∆ t]

where $∆ t$ is the time step, and ${δ W}_{n}$ obeys the normal distribution with mean zero and variance $∆ t$ . Then, the protein’s concentration is solved by the Euler scheme (taking A as an example):

{[A]}^{n + 1} = {[A]}^{n} + d t (k_{b a s a l} + \frac{\sum_{i} v_{i A} (1 + ε η_{i A}^{(n)}) {(\frac{x_{i}^{(n)}}{K_{i A}})}^{3} + δ_{A} (1 + ε ξ_{A}^{(n)})}{1 + \sum_{i} {(\frac{x_{i}^{(n)}}{K_{i A}})}^{3} + \sum_{i} {(\frac{y_{i}^{(n)}}{K_{i A}})}^{3}} - r_{A} (1 + ε ζ_{A}^{(n)}) A^{(n)}) .

In the presence of intrinsic noise, we also used the Milstein scheme to numerically solve the dynamics of proteins’ copy numbers. Taking $X_{A}$ as an example, its value at $n + 1$ time step is as follows:

X_{A}^{(n + 1)} = X_{A}^{(n)} + (V k_{b a s a l} + V \frac{\sum_{i} v_{i A} {(\frac{X_{i}^{(n)}}{K_{i A}})}^{3} + δ_{A} V^{3}}{V^{3} + \sum_{i} {(\frac{X_{i}^{(n)}}{K_{i A}})}^{3} + \sum_{i} {(\frac{Y_{i}^{(n)}}{K_{i A}})}^{3}} - r_{A} X_{A}^{(n)}) ∆ t + σ_{A} {δ W}_{A}^{(n)} + \frac{1}{2} σ_{A}^{2} [{({δ W}_{A}^{(n)})}^{2} - ∆ t],

σ_{A} = \sqrt{V k_{b a s a l} + V \frac{\sum_{i} v_{i A} {(\frac{X_{i}^{(n)}}{K_{i A}})}^{3} + δ_{A} V^{3}}{V^{3} + \sum_{i} {(\frac{X_{i}^{(n)}}{K_{i A}})}^{3} + \sum_{i} {(\frac{Y_{i}^{(n)}}{K_{i A}})}^{3}} + r_{A} X_{A}^{(n)}}

where

Searching for topologies capable of accurate oscillation

There are two steps for searching for topologies robustly capable of accurate oscillation. The first step is to select network topologies that are able to robustly oscillate among all two- and three-node network topologies. For each topology, 10,000 sets of kinetic parameters are assigned randomly, with ranges shown in Supplementary file 1a; for each parameter set, the initial state of the protein concentration is set to be 0, and we use ode15s in MATLAB to simulate the deterministic dynamics in the time interval [0, 1000]. The dynamics is regarded as oscillation if the following two requirements are satisfied: every protein concentration cannot maintain unchanged in the time interval [700, 1000]; peaks in the time interval [700, 1000] cannot differ a lot. The first requirement excludes the dynamics reaching the steady state, and the second the damping oscillator. We record the number of oscillatory dynamics for each topology, and then regard the topology with this number larger than 80 as the oscillatory topology. But for the topology with the repressilator, if the number of oscillatory dynamics exceeds 10, we still regard this topology as the oscillatory topology. This loose threshold ensures enough oscillatory topologies with the repressilator. In this way, we get 474 oscillatory topologies.

The second step is to explore the robustness of accurate oscillation for the above 474 oscillatory topologies. For each oscillatory topology, we sample enough parameter sets until there are 1000 parameter sets capable of oscillation. For each of these 1000 parameter sets, we record the period T and the amplitude (the maximal peak value among all protein concentrations) from deterministic behavior; next we simulate the stochastic behavior in the time interval [0, 100T]. In the presence of only extrinsic noise, the initial concentration is set to be the state when the concentration B reaches the highest value in a period, but in the presence of intrinsic noise, the initial concentration is converted to the copy number by multiplying the concentration with the cell volume $V$ . We use schemes mentioned in the previous sections to numerically solve the stochastic dynamics, with the time step in Supplementary file 1a. For a given noisy trajectory, the dimensionless autocorrelation time $τ / T$ is $- 1 / \log (c)$ , where $c$ is the autocorrelation coefficient at T. Since there are two or three trajectories each of which corresponds to a type of protein, so there are two or three dimensionless autocorrelation times, and we use the largest one as the final dimensionless autocorrelation time. Finally, we use the 90 percentile of dimensionless autocorrelation time to measure the robustness of this topology against noise. The 90 percentile is averaged over five repeated simulations.

Analytical results for the relation between robustness and period when tuning the timescale M

Phase noise in Demir et al.’s study

In this section, we briefly summarize Demir et al., 2000 study about the phase noise. We consider the dynamics described by the following equations:

\dot{x} = f (x) + B (x) ξ (t)

where $x \in R^{3}$ , $f (\cdot) : R^{3} \to R^{3}$ , $B (\cdot) : R^{3} \to R^{3 \times 3}$ and $ξ (t) \in R^{3}$ is random noise. Note that the noise amplitude $B$ is only related to $x$ , which is not affected by the time $t$ . The unperturbed system $\dot{x} = f (x)$ has a periodic solution $x_{s} (t)$ (with period $T$ ). Linearizing the noise-perturbed system around $x_{s} (t)$ gives the following system:

\dot{w} = A (t) w (t) + B (x_{s} (t)) ξ (t)

where $w (t) = x (t) - x_{s} (t)$ , $A (t) = {\frac{\partial f (x)}{\partial x} |}_{x = x_{s} (t)}$ is $T$ –periodic. From Floquet theory, the state transition matrix $Φ (t, s)$ for $\dot{w} = A (t) w (t)$ is given by

Φ (t, s) = U (t) \exp (D (t - s)) V (s) = \int_{i = 1}^{3} u_{i} (t) e x p (μ_{i} (t - s)) v_{i}^{T} (s)

where $U (t)$ is $T$ -periodic nonsingular matrix with columns denoted by $u_{i} (t)$ , and $V (t)$ with rows denoted by $v_{i}^{T} (t)$ is $U^{- 1} (t)$ . The $μ_{i}$ ’s are Floquet exponents. We can further choose $μ_{1} = 0$ and $u_{1} (t) = {\dot{x}}_{s} (t)$ . Then corresponding $v_{1} (t)$ will play an important role in calculating the phase noise.

From the nonlinear perturbation analysis, $z (t) = x_{s} (t + α (t)) + y (t)$ solves Equation 5 for a small $y (t)$ . The $α (t)$ and $y (t)$ are called as phase noise and deviation noise, respectively. It can be proved that the variance of the phase noise $α (t)$ increases linearly with time $t$ , that is,

V a r (α (t)) = c t

c = \frac{1}{T} \int_{0}^{T} v_{1}^{T} (t^{`}) B (x_{s} (t^{`})) B^{T} (x_{s} (t^{`})) v_{1} (t^{`}) d t^{`} .

where $c$ is given by

Since $c$ has the same unit as $T$ , we divide $c$ by $T$ to ensure a dimensionless index when measuring the phase noise.

The vector $v_{1} (t)$ when tuning $M$

We first consider the oscillator governed by the following equation:

\dot{x} = f (x)

We still use notations in the previous section to denote the quantities for this system. For example, the solution $x_{s} (t)$ , the period $T$ , the Jacobi $A (t)$ of $f (x)$ at the solution $x_{s} (t)$ . We assume the state transition matrix $Φ (t, s) = U (t) \exp (D (t - s)) V (s)$ , which satisfies the first Floquet exponent is 0 and the first column of $U (t)$ is the time derivative of $x_{s} (t)$ . Then the first row of $V (t)$ is denoted as $v_{1}^{T} (t)$ , which can be used to calculate the variance of phase noise.

Next, we explore how $v_{1} (t)$ changes when the right-hand term is divided by the timescale $M$ . By this way, we obtain

\begin{matrix} \dot{x} = f (x) \frac{1}{M} \end{matrix}

It is easy to verify that the system governed by Equation 10 has a periodic solution $x_{s} (t / M)$ with period $M T$ . The linearization of this system gives

\dot{w} = \frac{1}{M} A (\frac{t}{M}) w (t)

where $w (t) = x (t) - x_{s} (t / M)$ . According to the definition of $Φ (t, s)$ , $Φ (t, s)$ satisfies

\frac{d Φ (t, s)}{d t} = A (t) Φ (t, s), Φ (s, s) = I .

So $Φ (\frac{t}{M}, \frac{s}{M})$ satisfies

\frac{d Φ (\frac{t}{M}, \frac{s}{M})}{d t} = \frac{1}{M} A (\frac{t}{M}) Φ (\frac{t}{M}, \frac{s}{M}), Φ (\frac{s}{M}, \frac{s}{M}) = I .

Therefore, $Φ (\frac{t}{M}, \frac{s}{M})$ is the state transition matrix for $\dot{w} = \frac{1}{M} A (\frac{t}{M}) w (t)$ . Since $Φ (t, s) = U (t) \exp (D (t - s)) V (s)$ , we can obtain

Φ (\frac{t}{M}, \frac{s}{M}) = \frac{1}{M} U (\frac{t}{M}) \exp (\frac{D (t - s)}{M}) M V (\frac{s}{M})

where the first term $\frac{1}{M}$ in the right-hand term is to ensure the first column of $\frac{1}{M} U (\frac{t}{M})$ is the time derivative of $x_{s} (\frac{t}{M})$ , that is, $\frac{1}{M} {\dot{x}}_{s} (\frac{t}{M})$ . Thus, the first row of $M V (\frac{t}{M})$ is $M v_{1}^{T} (\frac{t}{M})$ , which can be used to calculate the variance of phase noise.

Oscillation accuracy against extrinsic noise when tuning the time scale

For the system governed by Equation 2, we add $\frac{1}{M}$ to the right-hand side and perturb the $v_{i j}, δ_{i}, r_{i}$ to introduce the extrinsic noise, thus leading to the following equations:

{\begin{array}{ll} \frac{d A}{d t} = (k_{b a s a l} + \frac{\sum_{i} v_{i A} {(\frac{x_{i}}{K_{i A}})}^{3} + δ_{A}}{1 + \sum_{i} {(\frac{x_{i}}{K_{i A}})}^{3} + \sum_{i} {(\frac{y_{i}}{K_{i A}})}^{3}} - r_{A} A) \frac{1}{M} + ε \frac{1}{M} (\frac{\sum_{i} v_{i A} η_{i A} {(\frac{x_{i}}{K_{i A}})}^{3} + δ_{A} ξ_{A}}{1 + \sum_{i} {(\frac{x_{i}}{K_{i A}})}^{3} + \sum_{i} {(\frac{y_{i}}{K_{i A}})}^{3}} - r_{A} A ζ_{A}) \\ \frac{d B}{d t} = (k_{b a s a l} + \frac{\sum_{i} v_{i B} {(\frac{x_{i}}{K_{i B}})}^{3} + δ_{B}}{1 + \sum_{i} {(\frac{x_{i}}{K_{i B}})}^{3} + \sum_{i} {(\frac{y_{i}}{K_{i B}})}^{3}} - r_{B} B) \frac{1}{M} + ε \frac{1}{M} (\frac{\sum_{i} v_{i B} η_{i B} {(\frac{x_{i}}{K_{i B}})}^{3} + δ_{B} ξ_{B}}{1 + \sum_{i} {(\frac{x_{i}}{K_{i B}})}^{3} + \sum_{i} {(\frac{y_{i}}{K_{i B}})}^{3}} - r_{B} B ζ_{B}) \\ \frac{d C}{d t} = (k_{b a s a l} + \frac{\sum_{i} v_{i C} {(\frac{x_{i}}{K_{i C}})}^{3} + δ_{C}}{1 + \sum_{i} {(\frac{x_{i}}{K_{i C}})}^{3} + \sum_{i} {(\frac{y_{i}}{K_{i C}})}^{3}} - r_{C} C) \frac{1}{M} + ε \frac{1}{M} (\frac{\sum_{i} v_{i C} η_{i C} {(\frac{x_{i}}{K_{i C}})}^{3} + δ_{C} ξ_{C}}{1 + \sum_{i} {(\frac{x_{i}}{K_{i C}})}^{3} + \sum_{i} {(\frac{y_{i}}{K_{i C}})}^{3}} - r_{C} C ζ_{C}) \end{array},

For simplicity, we still use $x$ and $f (x)$ to denote $(A, B, C)$ and the terms in the first brackets in the right-hand terms, respectively. Thus, the above equations can be rewritten as

\dot{x} = f (x) \frac{1}{M} + \frac{1}{M} B_{e x} (x) Π (t)

where $Π (t)$ = $(η_{A A}, η_{B A}, η_{C A,} ξ_{A}, ζ_{A}, η_{A B}, η_{B B}, η_{C B,} ξ_{B}, ζ_{B}, η_{A C}, η_{B C}, η_{C C,} ξ_{C}, ζ_{C})$ , and $B_{e x} (x)$ is the matrix whose elements are coefficients of the random noise when $M = 1$ . Recall that the ‘ $v_{1} (t)$ ’ in Equation 8 is $M v_{1} (\frac{t}{M})$ for the system $\dot{x} = f (x) \frac{1}{M}$ , so the slope of variance of phase noise over period for the system with timescale $M$ is

\frac{1}{M T} \frac{1}{M T} \int_{0}^{M T} M v_{1}^{T} (\frac{t^{`}}{M}) \frac{1}{M} B_{e x} (x_{s} (t^{`} / M)) \frac{1}{M} B_{e x}^{T} (x_{s} (t^{`} / M)) M v_{1} (\frac{t^{`}}{M}) d t^{`}

By replacing $t^{`}$ with $t^{` `} M$ and using $t^{`}$ to denote $t^{` `}$ , we obtain

\frac{1}{M T^{2}} \int_{0}^{T} v_{1}^{T} (t^{`}) B (x_{s} (t^{`})) B^{T} (x_{s} (t^{`})) v_{1} (t^{`}) d t^{`}

So it can be concluded that large $M$ causes small normalized phase noise in the presence of extrinsic noise. As large $M$ also leads to long period but has no effect on the amplitude, the long period might be the reason for high oscillation accuracy in the presence of extrinsic noise.

Oscillation accuracy against intrinsic noise when tuning the timescale

Similarly, for the system governed by Equation 2, we add $\frac{1}{M}$ to the right-hand side and introduce the cell volume $V$ to incorporate the intrinsic noise, thus leading to the following chemical Langevin equations:

{\begin{cases} d X_{A} = (V k_{b a s a l} + V \frac{\sum_{i} v_{i A} {(\frac{X_{i}}{K_{i A}})}^{3} + δ_{A} V^{3}}{V^{3} + \sum_{i} {(\frac{X_{i}}{K_{i A}})}^{3} + \sum_{i} {(\frac{Y_{i}}{K_{i A}})}^{3}} - r_{A} X_{A}) \frac{1}{M} d t + \sqrt{\frac{1}{M}} \sqrt{V k_{b a s a l} + V \frac{\sum_{i} v_{i A} {(\frac{X_{i}}{K_{i A}})}^{3} + δ_{A} V^{3}}{V^{3} + \sum_{i} {(\frac{X_{i}}{K_{i A}})}^{3} + \sum_{i} {(\frac{Y_{i}}{K_{i A}})}^{3}} + r_{A} X_{A}} d W_{t}^{A} \\ d X_{B} = (V k_{b a s a l} + V \frac{\sum_{i} v_{i B} {(\frac{X_{i}}{K_{i B}})}^{3} + δ_{B} V^{3}}{V^{3} + \sum_{i} {(\frac{X_{i}}{K_{i B}})}^{3} + \sum_{i} {(\frac{Y_{i}}{K_{i B}})}^{3}} - r_{B} X_{B}) \frac{1}{M} d t + \sqrt{\frac{1}{M}} \sqrt{V k_{b a s a l} + V \frac{\sum_{i} v_{i B} {(\frac{X_{i}}{K_{i B}})}^{3} + δ_{B} V^{3}}{V^{3} + \sum_{i} {(\frac{X_{i}}{K_{i B}})}^{3} + \sum_{i} {(\frac{Y_{i}}{K_{i B}})}^{3}} + r_{B} X_{B}} d W_{t}^{B} \\ d X_{C} = (V k_{b a s a l} + V \frac{\sum_{i} v_{i C} {(\frac{X_{i}}{K_{i C}})}^{3} + δ_{C} V^{3}}{V^{3} + \sum_{i} {(\frac{X_{i}}{K_{i C}})}^{3} + \sum_{i} {(\frac{Y_{i}}{K_{i C}})}^{3}} - r_{C} X_{C}) \frac{1}{M} d t + \sqrt{\frac{1}{M}} \sqrt{V k_{b a s a l} + V \frac{\sum_{i} v_{i C} {(\frac{X_{i}}{K_{i C}})}^{3} + δ_{C} V^{3}}{V^{3} + \sum_{i} {(\frac{X_{i}}{K_{i C}})}^{3} + \sum_{i} {(\frac{Y_{i}}{K_{i C}})}^{3}} + r_{C} X_{C}} d W_{t}^{C} \end{cases}

We use $X$ and $f (X)$ to denote $(X_{A}, X_{B}, X_{C})$ and the terms in the first brackets in the right-hand terms, respectively. Thus, the above equations can be rewritten as

\dot{X} = f (X) \frac{1}{M} + \sqrt{\frac{1}{M}} B_{i n} (X) Λ (t)

where $Λ (t)$ is $(d W_{A} / d t, d W_{B} / d t, d W_{C} / d t)$ , and $B_{i n} (X)$ is the matrix whose elements are coefficients of the random noise when $M = 1$ . If we use $v_{1}^{i n} (t)$ to represent the ‘ $v_{1} (t)$ ’ for the system $\dot{X} = f (X)$ , the ‘ $v_{1} (t)$ ’ for the system $\dot{X} = f (X) \frac{1}{M}$ is $M v_{1}^{i n} (\frac{t}{M})$ . So the slope of variance of phase noise over period for the system with cell volume $V$ is

\frac{1}{M T} \frac{1}{M T} \int_{0}^{M T} M {(v_{1}^{i n} (\frac{t^{`}}{M}))}^{T} \sqrt{\frac{1}{M}} B_{i n} (x_{s} (t^{`} / M)) \sqrt{\frac{1}{M}} B_{i n}^{T} (x_{s} (t^{`} / M)) M v_{1}^{i n} (\frac{t^{`}}{M}) d t^{`}

By replacing $t^{`}$ with $t^{` `} M$ and using $t^{`}$ to denote $t^{` `}$ , we obtain

\frac{1}{T^{2}} \int_{0}^{T} {(v_{1}^{i n} (t^{`}))}^{T} B (x_{s} (t^{`})) B^{T} (x_{s} (t^{`})) v_{1}^{i n} (t^{`}) d t^{`}

It can be seen that $M$ has no effect on the normalized slope of variance of phase noise, so the long period might not influence the noise in proteins when facing intrinsic noise.

Analytical results for the relation between robustness and amplitude when tuning the rescaling parameter $N$

Deterministic model with rescaled variables

To analyze the relation between amplitude and oscillation accuracy against noise, we replace $(A, B, C)$ in Equation 2 with $(\tilde{A} / N, \tilde{B} / N, \tilde{C} / N)$ , which allows us to tune the amplitude by varying $N$ . After this rescaling, we obtain the following equations for $\tilde{A}$ , $\tilde{B}$ and $\tilde{C}$ :

{\begin{cases} \frac{d \tilde{A}}{d t} = N k_{b a s a l} + \frac{\sum_{i} N v_{i A} {(\frac{{\tilde{x}}_{i}}{N K_{i A}})}^{3} + N δ_{A}}{1 + \sum_{i} {(\frac{{\tilde{x}}_{i}}{N K_{i A}})}^{3} + \sum_{i} {(\frac{{\tilde{y}}_{i}}{N K_{i A}})}^{3}} - r_{A} \tilde{A} \\ \frac{d \tilde{B}}{d t} = N k_{b a s a l} + \frac{\sum_{i} N v_{i B} {(\frac{{\tilde{x}}_{i}}{N K_{i B}})}^{3} + N δ_{B}}{1 + \sum_{i} {(\frac{{\tilde{x}}_{i}}{N K_{i B}})}^{3} + \sum_{i} {(\frac{{\tilde{y}}_{i}}{N K_{i B}})}^{3}} - r_{B} \tilde{B} \\ \frac{d \tilde{C}}{d t} = N k_{b a s a l} + \frac{\sum_{i} N v_{i C} {(\frac{{\tilde{x}}_{i}}{N K_{i C}})}^{3} + N δ_{C}}{1 + \sum_{i} {(\frac{{\tilde{x}}_{i}}{N K_{i C}})}^{3} + \sum_{i} {(\frac{{\tilde{y}}_{i}}{N K_{i C}})}^{3}} - r_{C} \tilde{C} \end{cases}

where ${\tilde{x}}_{i} = \frac{x_{i}}{N}, {\tilde{y}}_{i} = \frac{y_{i}}{N}$ . If $N = 1$ , this equation is the same as Equation 2, so $\tilde{A}$ , $\tilde{B}$ , and $\tilde{C}$ show same amplitudes with $A, B$ , and $C$ , respectively. Nevertheless, if $N \neq 1$ , the amplitude of $\tilde{A}$ , $\tilde{B}$ , or $\tilde{C}$ is $N$ times as high as that of $A, B$ , or $C$ . Note that $N$ has no effect on the period.

Oscillation accuracy against extrinsic noise when tuning the rescaling parameter

In the system governed by Equation 11, we assume that $N$ causes binding affinities, $v_{i j}$ ’s, and $δ_{i}$ ’s to $N$ times their original values, but the $r_{i}$ ’s remain unchanged. Next, we consider the system described by Equation 11 in the presence of only extrinsic noise. We perturb each kinetic parameter $v_{i j}, δ_{i}, r_{i}$ by the same method as mentioned in the section ‘Mathematical modeling’ and obtain a new system described by the following equations:

{\begin{cases} \frac{d \tilde{A}}{d t} = N k_{b a s a l} + \frac{\sum_{i} N v_{i A} {(\frac{{\tilde{x}}_{i}}{N K_{i A}})}^{3} + N δ_{A}}{1 + \sum_{i} {(\frac{{\tilde{x}}_{i}}{N K_{i A}})}^{3} + \sum_{i} {(\frac{{\tilde{y}}_{i}}{N K_{i A}})}^{3}} - r_{A} \tilde{A} + ε (\frac{\sum_{i} N v_{i A} η_{i A} {(\frac{{\tilde{x}}_{i}}{N K_{i A}})}^{3} + N δ_{A} ξ_{A}}{1 + \sum_{i} {(\frac{{\tilde{x}}_{i}}{N K_{i A}})}^{3} + \sum_{i} {(\frac{{\tilde{y}}_{i}}{N K_{i A}})}^{3}} - r_{A} \tilde{A} ζ_{A}) \\ \frac{d \tilde{B}}{d t} = N k_{b a s a l} + \frac{\sum_{i} N v_{i B} {(\frac{{\tilde{x}}_{i}}{N K_{i B}})}^{3} + N δ_{B}}{1 + \sum_{i} {(\frac{{\tilde{x}}_{i}}{N K_{i B}})}^{3} + \sum_{i} {(\frac{{\tilde{y}}_{i}}{N K_{i B}})}^{3}} - r_{B} \tilde{B} + ε (\frac{\sum_{i} N v_{i B} η_{i A} {(\frac{{\tilde{x}}_{i}}{N K_{i B}})}^{3} + N δ_{B} ξ_{A}}{1 + \sum_{i} {(\frac{{\tilde{x}}_{i}}{N K_{i B}})}^{3} + \sum_{i} {(\frac{{\tilde{y}}_{i}}{N K_{i B}})}^{3}} - r_{B} \tilde{B} ζ_{A}) \\ \frac{d \tilde{C}}{d t} = N k_{b a s a l} + \frac{\sum_{i} N v_{i C} {(\frac{{\tilde{x}}_{i}}{N K_{i C}})}^{3} + N δ_{C}}{1 + \sum_{i} {(\frac{{\tilde{x}}_{i}}{N K_{i C}})}^{3} + \sum_{i} {(\frac{{\tilde{y}}_{i}}{N K_{i C}})}^{3}} - r_{C} \tilde{C} + ε (\frac{\sum_{i} N v_{i C} η_{i A} {(\frac{{\tilde{x}}_{i}}{N K_{i C}})}^{3} + N δ_{C} ξ_{A}}{1 + \sum_{i} {(\frac{{\tilde{x}}_{i}}{N K_{i C}})}^{3} + \sum_{i} {(\frac{{\tilde{y}}_{i}}{N K_{i C}})}^{3}} - r_{C} \tilde{C} ζ_{A}) \end{cases}

where $η_{i j}$ , $ξ_{i}$ , $ζ_{i}$ are independent noise terms and all modeled by Equation 4 By multiplying $1 / N$ to both sides of above equations, we get

{\begin{cases} \frac{1}{N} \frac{d \tilde{A}}{d t} = k_{b a s a l} + \frac{\sum_{i} v_{i A} {(\frac{{\tilde{x}}_{i}}{N K_{i A}})}^{3} + δ_{A}}{1 + \sum_{i} {(\frac{{\tilde{x}}_{i}}{N K_{i A}})}^{3} + \sum_{i} {(\frac{{\tilde{y}}_{i}}{N K_{i A}})}^{3}} - r_{A} \frac{1}{N} \tilde{A} + ε (\frac{\sum_{i} v_{i A} η_{i A} {(\frac{{\tilde{x}}_{i}}{N K_{i A}})}^{3} + δ_{A} ξ_{A}}{1 + \sum_{i} {(\frac{{\tilde{x}}_{i}}{N K_{i A}})}^{3} + \sum_{i} {(\frac{{\tilde{y}}_{i}}{N K_{i A}})}^{3}} - r_{A} \frac{1}{N} \tilde{A} ζ_{A}) \\ \frac{1}{N} \frac{d \tilde{B}}{d t} = k_{b a s a l} + \frac{\sum_{i} v_{i B} {(\frac{{\tilde{x}}_{i}}{N K_{i B}})}^{3} + δ_{B}}{1 + \sum_{i} {(\frac{{\tilde{x}}_{i}}{N K_{i B}})}^{3} + \sum_{i} {(\frac{{\tilde{y}}_{i}}{N K_{i B}})}^{3}} - r_{B} \frac{1}{N} \tilde{B} + ε (\frac{\sum_{i} v_{i B} η_{i A} {(\frac{{\tilde{x}}_{i}}{N K_{i B}})}^{3} + δ_{B} ξ_{A}}{1 + \sum_{i} {(\frac{{\tilde{x}}_{i}}{N K_{i B}})}^{3} + \sum_{i} {(\frac{{\tilde{y}}_{i}}{N K_{i B}})}^{3}} - r_{B} \frac{1}{N} \tilde{B} ζ_{A}) \\ \frac{1}{N} \frac{d \tilde{C}}{d t} = k_{b a s a l} + \frac{\sum_{i} v_{i C} {(\frac{{\tilde{x}}_{i}}{N K_{i C}})}^{3} + δ_{C}}{1 + \sum_{i} {(\frac{{\tilde{x}}_{i}}{N K_{i C}})}^{3} + \sum_{i} {(\frac{{\tilde{y}}_{i}}{N K_{i C}})}^{3}} - r_{C} \frac{1}{N} \tilde{C} + ε (\frac{\sum_{i} v_{i C} η_{i A} {(\frac{{\tilde{x}}_{i}}{N K_{i C}})}^{3} + δ_{C} ξ_{A}}{1 + \sum_{i} {(\frac{{\tilde{x}}_{i}}{N K_{i C}})}^{3} + \sum_{i} {(\frac{{\tilde{y}}_{i}}{N K_{i C}})}^{3}} - r_{C} \frac{1}{N} \tilde{C} ζ_{A}) \end{cases}

Let $\tilde{\tilde{A}} = \frac{\tilde{A}}{N},$ $\tilde{\tilde{B}} = \frac{\tilde{B}}{N}$ , and $\tilde{\tilde{B}} = \frac{\tilde{B}}{N}$ , the set of equations for $\tilde{\tilde{A}}$ , $\tilde{\tilde{B}}$ , and $\tilde{\tilde{C}}$ is the same as that in Equation 3, in which $N$ do not appear. So the dynamics of $\tilde{\tilde{A}}$ , $\tilde{\tilde{B}}$ , or $\tilde{\tilde{C}}$ will not change with varied $N$ , thus leading to the same accuracy of oscillation when varying $N$ . Based on $\tilde{\tilde{A}} N = \tilde{A},$ $\tilde{\tilde{B}} N = \tilde{B}$ , and $\tilde{\tilde{C}} N = \tilde{C}$ and the fact that the rescaling has no effect on the correlation function, $\tilde{A}$ shows the same accuracy of oscillation as $\tilde{\tilde{A}}$ , and so does $\tilde{B}$ or $\tilde{C}$ . Therefore, in the system for $\tilde{A}$ , $\tilde{B}$ , and $\tilde{C}$ , its oscillation accuracy remains the same with varied $N$ . Since $N$ influences the amplitude while maintaining the period, the change in the amplitude will not affect the oscillation noise against extrinsic noise.

Oscillation accuracy against intrinsic noise when tuning the rescaling parameter

The dynamics of the system described by Equation 11. in the presence of only intrinsic noise is governed by

{\begin{cases} d X_{\tilde{A}} = (V N k_{b a s a l} + V \frac{\sum_{i} N v_{i A} {(\frac{X_{\tilde{i}}}{N K_{i A}})}^{3} + N δ_{A} V^{3}}{V^{3} + \sum_{i} {(\frac{X_{\tilde{i}}}{N K_{i A}})}^{3} + \sum_{i} {(\frac{Y_{\tilde{i}}}{N K_{i A}})}^{3}} - r_{A} X_{\tilde{A}}) d t \\ + \sqrt{V N k_{b a s a l} + V \frac{\sum_{i} N v_{i A} {(\frac{X_{\tilde{i}}}{N K_{i A}})}^{3} + N δ_{A} V^{3}}{V^{3} + \sum_{i} {(\frac{X_{\tilde{i}}}{N K_{i A}})}^{3} + \sum_{i} {(\frac{Y_{\tilde{i}}}{N K_{i A}})}^{3}} + r_{A} X_{\tilde{A}}} d W_{t}^{A} \\ d X_{\tilde{B}} = (V N k_{b a s a l} + V \frac{\sum_{i} N v_{i B} {(\frac{X_{\tilde{i}}}{N K_{i B}})}^{3} + N δ_{B} V^{3}}{V^{3} + \sum_{i} {(\frac{X_{\tilde{i}}}{N K_{i B}})}^{3} + \sum_{i} {(\frac{Y_{\tilde{i}}}{N K_{i B}})}^{3}} - r_{B} X_{\tilde{B}}) d t \\ + \sqrt{V N k_{b a s a l} + V \frac{\sum_{i} N v_{i B} {(\frac{X_{\tilde{i}}}{N K_{i B}})}^{3} + N δ_{B} V^{3}}{V^{3} + \sum_{i} {(\frac{X_{\tilde{i}}}{N K_{i B}})}^{3} + \sum_{i} {(\frac{Y_{\tilde{i}}}{N K_{i B}})}^{3}} + r_{B} X_{\tilde{B}}} d W_{t}^{B} \\ d X_{\tilde{C}} = (V N k_{b a s a l} + V \frac{\sum_{i} N v_{i C} {(\frac{X_{\tilde{i}}}{N K_{i C}})}^{3} + N δ_{C} V^{3}}{V^{3} + \sum_{i} {(\frac{X_{\tilde{i}}}{N K_{i C}})}^{3} + \sum_{i} {(\frac{Y_{\tilde{i}}}{N K_{i C}})}^{3}} - r_{C} X_{\tilde{C}}) d t \\ + \sqrt{V N k_{b a s a l} + V \frac{\sum_{i} N v_{i C} {(\frac{X_{\tilde{i}}}{N K_{i C}})}^{3} + N δ_{C} V^{3}}{V^{3} + \sum_{i} {(\frac{X_{\tilde{i}}}{N K_{i C}})}^{3} + \sum_{i} {(\frac{Y_{\tilde{i}}}{N K_{i C}})}^{3}} + r_{C} X_{\tilde{C}}} d W_{t}^{C} \end{cases}

where $X_{\tilde{A}} = V \tilde{A}$ , $X_{\tilde{B}} = V \tilde{B}$ , and $X_{\tilde{C}} = V \tilde{C}$ , and $V$ is the cell volume. By multiplying $1 / N$ to both sides of above equations, we get

{\begin{cases} \frac{1}{N} d X_{\tilde{A}} = (V k_{b a s a l} + V \frac{\sum_{i} v_{i A} {(\frac{X_{\tilde{i}}}{N K_{i A}})}^{3} + δ_{A} V^{3}}{V^{3} + \sum_{i} {(\frac{X_{\tilde{i}}}{N K_{i A}})}^{3} + \sum_{i} {(\frac{Y_{\tilde{i}}}{N K_{i A}})}^{3}} - r_{A} \frac{1}{N} X_{\tilde{A}}) d t \\ + \sqrt{\frac{1}{N}} \sqrt{V k_{b a s a l} + V \frac{\sum_{i} v_{i A} {(\frac{X_{\tilde{i}}}{N K_{i A}})}^{3} + δ_{A} V^{3}}{V^{3} + \sum_{i} {(\frac{X_{\tilde{i}}}{N K_{i A}})}^{3} + \sum_{i} {(\frac{Y_{\tilde{i}}}{N K_{i A}})}^{3}} + r_{A} \frac{1}{N} X_{\tilde{A}}} d W_{t}^{A} \\ \frac{1}{N} d X_{\tilde{B}} = (V k_{b a s a l} + V \frac{\sum_{i} v_{i B} {(\frac{X_{\tilde{i}}}{N K_{i B}})}^{3} + δ_{B} V^{3}}{V^{3} + \sum_{i} {(\frac{X_{\tilde{i}}}{N K_{i B}})}^{3} + \sum_{i} {(\frac{Y_{\tilde{i}}}{N K_{i B}})}^{3}} - r_{B} \frac{1}{N} X_{\tilde{B}}) d t \\ + \sqrt{\frac{1}{N}} \sqrt{V k_{b a s a l} + V \frac{\sum_{i} v_{i B} {(\frac{X_{\tilde{i}}}{N K_{i B}})}^{3} + δ_{B} V^{3}}{V^{3} + \sum_{i} {(\frac{X_{\tilde{i}}}{N K_{i B}})}^{3} + \sum_{i} {(\frac{Y_{\tilde{i}}}{N K_{i B}})}^{3}} + r_{B} \frac{1}{N} X_{\tilde{B}}} d W_{t}^{B} \\ \frac{1}{N} d X_{\tilde{C}} = (V k_{b a s a l} + V \frac{\sum_{i} v_{i C} {(\frac{X_{\tilde{i}}}{N K_{i C}})}^{3} + δ_{C} V^{3}}{V^{3} + \sum_{i} {(\frac{X_{\tilde{i}}}{N K_{i C}})}^{3} + \sum_{i} {(\frac{Y_{\tilde{i}}}{N K_{i C}})}^{3}} - r_{C} \frac{1}{N} X_{\tilde{C}}) d t \\ + \sqrt{\frac{1}{N}} \sqrt{V k_{b a s a l} + V \frac{\sum_{i} v_{i C} {(\frac{X_{\tilde{i}}}{N K_{i C}})}^{3} + δ_{C} V^{3}}{V^{3} + \sum_{i} {(\frac{X_{\tilde{i}}}{N K_{i C}})}^{3} + \sum_{i} {(\frac{Y_{\tilde{i}}}{N K_{i C}})}^{3}} + r_{C} \frac{1}{N} X_{\tilde{C}}} d W_{t}^{C} \end{cases}

Let $X_{\tilde{\tilde{A}}} = \frac{X_{\tilde{A}}}{N},$ $X_{\tilde{\tilde{B}}} = \frac{X_{\tilde{B}}}{N},$ and $X_{\tilde{\tilde{C}}} = \frac{X_{\tilde{C}}}{N}$ , equations for $X_{\tilde{\tilde{A}}}$ , $X_{\tilde{\tilde{B}}}$ , and $X_{\tilde{\tilde{C}}}$ are

{\begin{cases} d X_{\tilde{\tilde{A}}} = (V k_{b a s a l} + V \frac{\sum_{i} v_{i A} {(\frac{X_{\tilde{\tilde{i}}}}{K_{i A}})}^{3} + δ_{A} V^{3}}{V^{3} + \sum_{i} {(\frac{X_{\tilde{\tilde{i}}}}{K_{i A}})}^{3} + \sum_{i} {(\frac{Y_{\tilde{\tilde{i}}}}{K_{i A}})}^{3}} - r_{A} X_{\tilde{\tilde{A}}}) d t \\ + \sqrt{\frac{1}{N}} \sqrt{V k_{b a s a l} + V \frac{\sum_{i} v_{i A} {(\frac{X_{\tilde{\tilde{i}}}}{K_{i A}})}^{3} + δ_{A} V^{3}}{V^{3} + \sum_{i} {(\frac{X_{\tilde{\tilde{i}}}}{K_{i A}})}^{3} + \sum_{i} {(\frac{Y_{\tilde{\tilde{i}}}}{K_{i A}})}^{3}} + r_{A} X_{\tilde{\tilde{A}}}} d W_{t}^{A} \\ d X_{\tilde{\tilde{B}}} = (V k_{b a s a l} + V \frac{\sum_{i} v_{i B} {(\frac{X_{\tilde{\tilde{i}}}}{K_{i B}})}^{3} + δ_{B} V^{3}}{V^{3} + \sum_{i} {(\frac{X_{\tilde{\tilde{i}}}}{K_{i B}})}^{3} + \sum_{i} {(\frac{Y_{\tilde{\tilde{i}}}}{K_{i B}})}^{3}} - r_{B} X_{\tilde{\tilde{B}}}) d t \\ + \sqrt{\frac{1}{N}} \sqrt{V k_{b a s a l} + V \frac{\sum_{i} v_{i B} {(\frac{X_{\tilde{\tilde{i}}}}{K_{i B}})}^{3} + δ_{B} V^{3}}{V^{3} + \sum_{i} {(\frac{X_{\tilde{\tilde{i}}}}{K_{i B}})}^{3} + \sum_{i} {(\frac{Y_{\tilde{\tilde{i}}}}{K_{i B}})}^{3}} + r_{B} X_{\tilde{\tilde{B}}}} d W_{t}^{B} \\ d X_{\tilde{\tilde{C}}} = (V k_{b a s a l} + V \frac{\sum_{i} v_{i C} {(\frac{X_{\tilde{\tilde{i}}}}{K_{i C}})}^{3} + δ_{C} V^{3}}{V^{3} + \sum_{i} {(\frac{X_{\tilde{\tilde{i}}}}{K_{i C}})}^{3} + \sum_{i} {(\frac{Y_{\tilde{\tilde{i}}}}{K_{i C}})}^{3}} - r_{C} X_{\tilde{\tilde{C}}}) d t \\ + \sqrt{\frac{1}{N}} \sqrt{V k_{b a s a l} + V \frac{\sum_{i} v_{i C} {(\frac{X_{\tilde{\tilde{i}}}}{K_{i C}})}^{3} + δ_{C} V^{3}}{V^{3} + \sum_{i} {(\frac{X_{\tilde{\tilde{i}}}}{K_{i C}})}^{3} + \sum_{i} {(\frac{Y_{\tilde{\tilde{i}}}}{K_{i C}})}^{3}} + r_{C} X_{\tilde{\tilde{C}}}} d W_{t}^{C} \end{cases}

In above equations, $N$ only negatively affects the magnitude of noise term, so the oscillation accuracies of $X_{\tilde{\tilde{A}}}$ , $X_{\tilde{\tilde{B}}}$ , and $X_{\tilde{\tilde{C}}}$ increased with increased $N$ . Thus, the oscillation accuracies of $X_{\tilde{A}}$ , $X_{\tilde{B}}$ , and $X_{\tilde{C}}$ also increased with increased $N$ because the correlation function is not affected by the rescaling operation. Besides, large $N$ increase the amplitude while maintaining the period. Taken together, the high amplitude may enhance the oscillation noise against intrinsic noise.

Data availability

The current manuscript is a computational study. Modelling code and NFKB data for plotting are uploaded to GitHub at https://github.com/LingxiaQiao/oscillation, (copy archived at swh:1:rev:72a2d3d1146b14e7988c1cc06208fe1252e9a6f5).

References

(2003) Development of genetic circuitry exhibiting toggle switch or oscillatory behavior in Escherichia coli
Cell 113:597–607.

https://doi.org/10.1016/s0092-8674(03)00346-5
- PubMed
- Google Scholar
(2008) Frequency-modulated nuclear localization bursts coordinate gene regulation
Nature 455:485–490.

https://doi.org/10.1038/nature07292
- PubMed
- Google Scholar
1. Cao Y
2. Lu HM
3. Liang J
(2010) Probability landscape of heritable and robust epigenetic state of lysogeny in phage lambda
PNAS 107:18445–18450.

https://doi.org/10.1073/pnas.1001455107
- PubMed
- Google Scholar
1. Cao Y
2. Wang H
3. Ouyang Q
4. Tu Y
(2015) The free energy cost of accurate biochemical oscillations
Nature Physics 11:772–778.

https://doi.org/10.1038/nphys3412
- PubMed
- Google Scholar
1. Cao Y
2. Terebus A
3. Liang J
(2016) Accurate chemical master equation solution using multi-finite buffers
Multiscale Modeling & Simulation 14:923–963.

https://doi.org/10.1137/15M1034180
- Google Scholar
1. Cao Y
2. Lei X
3. Ribeiro RM
4. Perelson AS
5. Liang J
(2018) Probabilistic control of HIV latency and transactivation by the tat gene circuit
PNAS 115:12453–12458.

https://doi.org/10.1073/pnas.1811195115
- Google Scholar
1. Chen Y
2. Kim JK
3. Hirning AJ
4. Josić K
5. Bennett MR
(2015) SYNTHETIC biology: emergent genetic oscillations in a SYNTHETIC microbial consortium
Science 349:986–989.

https://doi.org/10.1126/science.aaa3794
- PubMed
- Google Scholar
1. Choi PJ
2. Cai L
3. Frieda K
4. Xie XS
(2008) A stochastic single-molecule event triggers phenotype switching of A bacterial cell
Science 322:442–446.

https://doi.org/10.1126/science.1161427
- PubMed
- Google Scholar
(2000) Phase noise in oscillators: a unifying theory and numerical methods for characterization
IEEE Transactions on Circuits and Systems I 47:655–674.

https://doi.org/10.1109/81.847872
- Google Scholar
(2016) Waves of cdk1 activity in S phase synchronize the cell cycle in Drosophila embryos
Developmental Cell 38:399–412.

https://doi.org/10.1016/j.devcel.2016.07.023
- PubMed
- Google Scholar
1. Duncan A
2. Liao S
3. Vejchodský T
4. Erban R
5. Grima R
(2015) Noise-induced multistability in chemical systems: discrete versus continuum modeling
Physical Review. E, Statistical, Nonlinear, and Soft Matter Physics 91:042111.

https://doi.org/10.1103/PhysRevE.91.042111
- PubMed
- Google Scholar
1. Elowitz MB
2. Leibler S
(2000) A synthetic oscillatory network of transcriptional regulators
Nature 403:335–338.

https://doi.org/10.1038/35002125
- PubMed
- Google Scholar
(2002) Stochastic gene expression in a single cell
Science 297:1183–1186.

https://doi.org/10.1126/science.1070919
- PubMed
- Google Scholar
1. Fei C
2. Cao Y
3. Ouyang Q
4. Tu Y
(2018) Design principles for enhancing phase sensitivity and suppressing phase fluctuations simultaneously in biochemical oscillatory systems
Nature Communications 9:1434.

https://doi.org/10.1038/s41467-018-03826-4
- PubMed
- Google Scholar
(2011) Modeling the cell cycle: why do certain circuits oscillate?
Cell 144:874–885.

https://doi.org/10.1016/j.cell.2011.03.006
- PubMed
- Google Scholar
1. Forger DB
2. Peskin CS
(2005) Stochastic simulation of the mammalian circadian clock
PNAS 102:321–324.

https://doi.org/10.1073/pnas.0408465102
- PubMed
- Google Scholar
1. Gillespie DT
(2000) The chemical langevin equation
J Chem Phys 113:297–306.

https://doi.org/10.1063/1.481811
- Google Scholar
1. Gillespie DT
2. Petzold LR
(2003) Improved leap-size selection for accelerated stochastic simulation
J Chem Phys 119:8229–8234.

https://doi.org/10.1063/1.1613254
- Google Scholar
1. Glass L
2. Pasternack JS
(1978) Stable oscillations in mathematical models of biological control systems
Journal of Mathematical Biology 6:207–223.

https://doi.org/10.1007/BF02547797
- Google Scholar
(2002a) Biochemical clocks and molecular noise: theoretical study of robustness factors
J Chem Phys 116:10997–11010.

https://doi.org/10.1063/1.1475765
- Google Scholar
(2002b) Deterministic versus stochastic models for circadian rhythms
Journal of Biological Physics 28:637–653.

https://doi.org/10.1023/A:1021286607354
- PubMed
- Google Scholar
1. Gonze D
2. Goldbeter A
(2006) Circadian rhythms and molecular noise
Chaos 16:026110.

https://doi.org/10.1063/1.2211767
- PubMed
- Google Scholar
1. Hao N
2. O’Shea EK
(2012) Signal-dependent dynamics of transcription factor translocation controls gene expression
Nature Structural & Molecular Biology 19:31–39.

https://doi.org/10.1038/nsmb.2192
- PubMed
- Google Scholar
(2019) On chaotic dynamics in transcription factors and the associated effects in differential gene regulation
Nature Communications 10:71.

https://doi.org/10.1038/s41467-018-07932-1
- PubMed
- Google Scholar
1. Hornung G
2. Barkai N
(2008) Noise propagation and signaling sensitivity in biological networks: a role for positive feedback
PLOS Computational Biology 4:e8.

https://doi.org/10.1371/journal.pcbi.0040008
- PubMed
- Google Scholar
(2014) The validity of quasi-steady-state approximations in discrete stochastic simulations
Biophysical Journal 107:783–793.

https://doi.org/10.1016/j.bpj.2014.06.012
- PubMed
- Google Scholar
Book
1. Kloeden PE
2. Platen E
(2013)
Numerical Solution of Stochastic Differential Equations

Berlin: Springer-Verlag.
- Google Scholar
1. Li X
2. Cassidy JJ
3. Reinke CA
4. Fischboeck S
5. Carthew RW
(2009) A microrna imparts robustness against environmental fluctuation during development
Cell 137:273–282.

https://doi.org/10.1016/j.cell.2009.01.058
- PubMed
- Google Scholar
1. Lim WA
2. Lee CM
3. Tang C
(2013) Design principles of regulatory networks: searching for the molecular algorithms of the cell
Molecular Cell 49:202–212.

https://doi.org/10.1016/j.molcel.2012.12.020
- PubMed
- Google Scholar
(2004) Effects of noise in excitable systems
Physics Reports 392:321–424.

https://doi.org/10.1016/j.physrep.2003.10.015
- Google Scholar
Preprint
1. Liu L
2. Wang P
3. Nair MS
4. Yu J
5. Rapp M
6. Wang Q
7. Luo Y
8. Chan JFW
9. Sahi V
10. Figueroa A
11. Guo XV
12. Cerutti G
13. Bimela J
14. Gorman J
15. Zhou T
16. Chen Z
17. Yuen KY
18. Kwong PD
19. Sodroski JG
20. Yin MT
21. Sheng Z
22. Huang Y
23. Shapiro L
24. Ho DD
(2020) Potent Neutralizing Antibodies Directed to Multiple Epitopes on Sars-Cov-2 Spike
bioRxiv.

https://doi.org/10.1101/2020.06.17.153486
- Google Scholar
1. Ma W
2. Trusina A
3. El-Samad H
4. Lim WA
5. Tang C
(2009) Defining network topologies that can achieve biochemical adaptation
Cell 138:760–773.

https://doi.org/10.1016/j.cell.2009.06.013
- PubMed
- Google Scholar
(2009) Delay-induced degrade-and-fire oscillations in small genetic circuits
Physical Review Letters 102:068105.

https://doi.org/10.1103/PhysRevLett.102.068105
- PubMed
- Google Scholar
(2018) Robustness of clocks to input noise
Physical Review Letters 121:078101.

https://doi.org/10.1103/PhysRevLett.121.078101
- PubMed
- Google Scholar
1. Munsky B
2. Khammash M
(2006) The finite state projection algorithm for the solution of the chemical master equation
J Chem Phys 124:044104.

https://doi.org/10.1063/1.2145882
- PubMed
- Google Scholar
1. Nie Q
2. Qiao L
3. Qiu Y
4. Zhang L
5. Zhao W
(2020) Noise control and utility: from regulatory network to spatial patterning
Science China Mathematics 63:425–440.

https://doi.org/10.1007/s11425-019-1633-1
- Google Scholar
1. Novák B
2. Tyson JJ
(2008) Design principles of biochemical oscillators
Nature Reviews. Molecular Cell Biology 9:981–991.

https://doi.org/10.1038/nrm2530
- PubMed
- Google Scholar
(2014) Molecular architecture of the mammalian circadian clock
Trends in Cell Biology 24:90–99.

https://doi.org/10.1016/j.tcb.2013.07.002
- PubMed
- Google Scholar
1. Paulsson J
(2004) Summing up the noise in gene networks
Nature 427:415–418.

https://doi.org/10.1038/nature02257
- PubMed
- Google Scholar
(2016) Synchronous long-term oscillations in a synthetic gene circuit
Nature 538:514–517.

https://doi.org/10.1038/nature19841
- PubMed
- Google Scholar
1. Qiao L
2. Zhao W
3. Tang C
4. Nie Q
5. Zhang L
(2019) Network topologies that can achieve dual function of adaptation and noise attenuation
Cell Systems 9:271–285.

https://doi.org/10.1016/j.cels.2019.08.006
- PubMed
- Google Scholar
Software
1. Qiao L
(2022) Oscillation, version swh:1:rev:72a2d3d1146b14e7988c1cc06208fe1252e9a6f5
Software Heritage.

https://archive.softwareheritage.org/swh:1:dir:58ab93976660d8a97ad3327aee3df93c774b12d6;origin=https://github.com/LingxiaQiao/oscillation;visit=swh:1:snp:0ebee09d7315cf9147b78e307079d8c7cc676299;anchor=swh:1:rev:72a2d3d1146b14e7988c1cc06208fe1252e9a6f5
1. Qin BW
2. Zhao L
3. Lin W
(2021) A frequency-amplitude coordinator and its optimal energy consumption for biological oscillators
Nature Communications 12:5894.

https://doi.org/10.1038/s41467-021-26182-2
- PubMed
- Google Scholar
1. Raser JM
2. O’Shea EK
(2004) Control of stochasticity in eukaryotic gene expression
Science 304:1811–1814.

https://doi.org/10.1126/science.1098641
- PubMed
- Google Scholar
(2011) Legitimacy of the stochastic michaelis-menten approximation
IET Systems Biology 5:58.

https://doi.org/10.1049/iet-syb.2009.0057
- PubMed
- Google Scholar
1. Sartori P
2. Tu Y
(2011) Noise filtering strategies in adaptive biochemical signaling networks: application to E. coli chemotaxis
Journal of Statistical Physics 142:1206–1217.

https://doi.org/10.1007/s10955-011-0169-z
- PubMed
- Google Scholar
1. Shahrezaei V
2. Swain PS
(2008) Analytical distributions for stochastic gene expression
PNAS 105:17256–17261.

https://doi.org/10.1073/pnas.0803850105
- PubMed
- Google Scholar
1. Shi W
2. Ma W
3. Xiong L
4. Zhang M
5. Tang C
(2017) Adaptation with transcriptional regulation
Scientific Reports 7:42648.

https://doi.org/10.1038/srep42648
- PubMed
- Google Scholar
1. Stricker J
2. Cookson S
3. Bennett MR
4. Mather WH
5. Tsimring LS
6. Hasty J
(2008) A fast, robust and tunable synthetic gene oscillator
Nature 456:516–519.

https://doi.org/10.1038/nature07389
- PubMed
- Google Scholar
(2002) Intrinsic and extrinsic contributions to stochasticity in gene expression
PNAS 99:12795–12800.

https://doi.org/10.1073/pnas.162041399
- PubMed
- Google Scholar
1. Terebus A
2. Liu C
3. Liang J
(2019) Discrete and continuous models of probability flux of switching dynamics: uncovering stochastic oscillations in a toggle-switch system
J Chem Phys 151:185104.

https://doi.org/10.1063/1.5124823
- PubMed
- Google Scholar
(2021) Exact probability landscapes of stochastic phenotype switching in feed-forward loops: phase diagrams of multimodality
Frontiers in Genetics 12:645640.

https://doi.org/10.3389/fgene.2021.645640
- PubMed
- Google Scholar
1. Thattai M
2. van Oudenaarden A
(2001) Intrinsic noise in gene regulatory networks
PNAS 98:8614–8619.

https://doi.org/10.1073/pnas.151588598
- PubMed
- Google Scholar
(2010) A synthetic low-frequency mammalian oscillator
Nucleic Acids Research 38:2702–2711.

https://doi.org/10.1093/nar/gkq121
- PubMed
- Google Scholar
1. Tyson JJ
(1991) Modeling the cell division cycle: cdc2 and cyclin interactions
PNAS 88:7328–7332.

https://doi.org/10.1073/pnas.88.16.7328
- PubMed
- Google Scholar
1. Veliz-Cuba A
2. Hirning AJ
3. Atanas AA
4. Hussain F
5. Vancia F
6. Josić K
7. Bennett MR
(2015) Sources of variability in a synthetic gene oscillator
PLOS Computational Biology 11:e1004674.

https://doi.org/10.1371/journal.pcbi.1004674
- PubMed
- Google Scholar
1. Vellela M
2. Qian H
(2009) Stochastic dynamics and non-equilibrium thermodynamics of a bistable chemical system: the schlögl model revisited
Journal of the Royal Society, Interface 6:925–940.

https://doi.org/10.1098/rsif.2008.0476
- PubMed
- Google Scholar
1. Wang L
2. Xin J
3. Nie Q
(2010) A critical quantity for noise attenuation in feedback systems
PLOS Computational Biology 6:e1000764.

https://doi.org/10.1371/journal.pcbi.1000764
- PubMed
- Google Scholar
1. Wang DG
2. Wang S
3. Huang B
4. Liu F
(2019) Roles of cellular heterogeneity, intrinsic and extrinsic noise in variability of p53 oscillation
Scientific Reports 9:5883.

https://doi.org/10.1038/s41598-019-41904-9
- PubMed
- Google Scholar
1. Yu P
2. Nie Q
3. Tang C
4. Zhang L
(2018) Nanog induced intermediate state in regulating stem cell differentiation and reprogramming
BMC Systems Biology 12:22.

https://doi.org/10.1186/s12918-018-0552-3
- PubMed
- Google Scholar
(2005) Noise induced complexity: from subthreshold oscillations to spiking in coupled excitable systems
Chaos 15:26117.

https://doi.org/10.1063/1.1886386
- PubMed
- Google Scholar
1. Zhang ZB
2. Wang QY
3. Ke YX
4. Liu SY
5. Ju JQ
6. Lim WA
7. Tang C
8. Wei P
(2017) Design of tunable oscillatory dynamics in a synthetic NF-κB signaling circuit
Cell Systems 5:460–470.

https://doi.org/10.1016/j.cels.2017.09.016
- PubMed
- Google Scholar
1. Zhang MY
2. Tang C
(2019) Bi-functional biochemical networks
Physical Biology 16:016001.

https://doi.org/10.1088/1478-3975/aae74c
- Google Scholar
1. Zhao W
2. Qiao L
3. Yan S
4. Nie Q
5. Zhang L
(2021) Mathematical modeling of histone modifications reveals the formation mechanism and function of bivalent chromatin
IScience 24:102732.

https://doi.org/10.1016/j.isci.2021.102732
- PubMed
- Google Scholar

Article and author information

Author details

Lingxia Qiao

Beijing International Center for Mathematical Research, Peking University, Beijing, China

Contribution
Formal analysis, Validation, Investigation, Methodology, Writing – original draft, Writing – review and editing

Contributed equally with
Zhi-Bo Zhang and Wei Zhao

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0003-2639-6851
Zhi-Bo Zhang
1. Center for Quantitative Biology, Peking University, Beijing, China
2. Peking-Tsinghua Joint Center for Life Sciences, Academy for Advanced Interdisciplinary Studies, Peking University, Beijing, China
Contribution
Formal analysis, Validation, Investigation, Methodology, Writing – original draft, Writing – review and editing

Contributed equally with
Lingxia Qiao and Wei Zhao

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-3316-2634
Wei Zhao

Center for Quantitative Biology, Peking University, Beijing, China

Contribution
Formal analysis, Validation, Investigation, Methodology, Writing – original draft, Writing – review and editing

Contributed equally with
Lingxia Qiao and Zhi-Bo Zhang

Competing interests
No competing interests declared
Ping Wei
1. Center for Quantitative Biology, Peking University, Beijing, China
2. Center for Cell and Gene Circuit Design, CAS Key Laboratory of Quantitative Engineering Biology, Shenzhen Institute of Synthetic Biology, Shenzhen Institutes of Advanced Technology, Chinese Academy of Sciences, Shenzhen, China
Contribution
Supervision, Funding acquisition, Project administration, Writing – review and editing

For correspondence
ping.wei@siat.ac.cn

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-1275-7145
Lei Zhang
1. Beijing International Center for Mathematical Research, Peking University, Beijing, China
2. Center for Quantitative Biology, Peking University, Beijing, China
Contribution
Conceptualization, Supervision, Funding acquisition, Methodology, Project administration, Writing – review and editing

For correspondence
zhangl@math.pku.edu.cn

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0001-9972-2051

Funding

National Key Research and Development Program of China (2021YFF1200500)

Lei Zhang

National Natural Science Foundation of China (12050002)

Lei Zhang

National Key Basic Research Program of China (2018YFA0902800)

Ping Wei

National Natural Science Foundation of China (31622022)

Ping Wei

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

LZ was partly supported by the National Key Research and Development Program of China 2021YFF1200500 and National Natural Science Foundation of China No. 12050002. PW was partly supported by the National Key Basic Research Program of China 2018YFA0902800 and the National Natural Science Foundation of China 31622022.

Copyright

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.