Network design principle for robust oscillatory behaviors with respect to biological noise
Abstract
Oscillatory behaviors, which are ubiquitous in transcriptional regulatory networks, are often subject to inevitable biological noise. Thus, a natural question is how transcriptional regulatory networks can robustly achieve accurate oscillation in the presence of biological noise. Here, we search all two and threenode transcriptional regulatory network topologies for those robustly capable of accurate oscillation against the parameter variability (extrinsic noise) or stochasticity of chemical reactions (intrinsic noise). We find that, no matter what source of the noise is applied, the topologies containing the repressilator with positive autoregulation show higher robustness of accurate oscillation than those containing the activatorinhibitor oscillator, and additional positive autoregulation enhances the robustness against noise. Nevertheless, the attenuation of different sources of noise is governed by distinct mechanisms: the parameter variability is buffered by the long period, while the stochasticity of chemical reactions is filtered by the high amplitude. Furthermore, we analyze the noise of a synthetic human nuclear factor κB (NFκB) signaling network by varying three different topologies and verify that the addition of a repressilator to the activatorinhibitor oscillator, which leads to the emergence of highrobustness motif—the repressilator with positive autoregulation—improves the oscillation accuracy in comparison to the topology with only an activatorinhibitor oscillator. These design principles may be applicable to other oscillatory circuits.
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.sa0Introduction
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; PotvinTrottier et al., 2016; Stricker et al., 2008; Tigges et al., 2010; Zhang et al., 2017). Repressilator (Elowitz and Leibler, 2000) and activatorinhibitor 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; PotvinTrottier 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 (PotvinTrottier 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 parameterdetermined 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 bottomup 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 threenode networks. We first perform an exhausting search of two and threenode 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 activatorinhibitor 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 $\tau $ divided by the period $T$. The correlation time $\mathrm{\tau}$ describes how fast the autocorrelation function $C\left(t\right)$ exponentially decays. To be specific, for a noisy dynamic trajectory of the oscillator, $C\left(t\right)$ displays a damped oscillation (Figure 1A):
where ${\u2329\cdots \u232a}_{s}$ is defined by ${\u2329f\left(s\right)\u232a}_{s}=\underset{S\to \infty}{\mathrm{lim}}\frac{1}{S}{\int}_{0}^{S}f\left(s\right)ds$, and $\u2329\cdots \u232a$ 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 $\tau $. The correlation time $\tau $ has the same unit as that of the period, so $\tau /T$ is dimensionless. Therefore, instead of using $\tau $, we utilize $\tau /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 bottomup 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 threenode networks might be enough to capture key features. Besides, the number of all two and threenode network topologies is 3^{9} 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 threenode network topologies is reduced from 3^{9} to 1955. Here, two typical oscillatory topologies are shown in Figure 1C: the activatorinhibitor and repressilator topologies. For the activatorinhibitor 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 negativefeedback loop.
Mathematical modeling
To model two and threenode 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}_{basal}$ is the basal production rate (much smaller than other terms); ${v}_{AA}$ is the maximum production rate caused by product A; ${K}_{AA}$ and ${K}_{BA}$ 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; VelizCuba 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 threenode network topologies that can robustly achieve accurate oscillation (i.e., high dimensionless correlation time $\tau /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 90percentile value of $\tau /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 90percentile value of $\tau /T$ is defined as the value of $\tau /T$ below which 90% of $\tau /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 90percentile value of $\tau /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 activatorinhibitor 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).
The topologies containing repressilator with positive autoregulation perform better than those containing the activatorinhibitor 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 $\tau /T$ in the presence of extrinsic noise for the five categories, where each violin corresponds to one category. By applying onetailed Wilcoxon ranksum tests to adjacent two categories, we find that 90 percentiles of $\tau /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 activatorinhibitor 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 $\tau /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 highrobustness core motif usually improves the robustness, but the combined topology cannot outperform the added highrobustness 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 $\mathrm{\Omega}$expansion to approximate the birthanddeath 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}_{opt}$ and ${A}_{opt}$ , respectively. Here, we refer to the amplitude as the maximal peak value among nodes A–C; ${T}_{opt}$ (or ${A}_{opt}$) is defined as the expectation of the bestfit exponential distribution of 1000 periods (or 1000 amplitudes) (Figure 2D). Therefore, the topology with large ${T}_{opt}$ tends to oscillate with long period, and the topology with large ${A}_{opt}$ usually indicates an oscillation with high amplitude. Note that these two quantities are calculated in the noisefree 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}_{opt}$ and ${A}_{opt}$ in the robustness of accurate oscillation against extrinsic noise, we calculate Spearman coefficients between these two quantities and 90 percentiles of $\tau /T$ for all 474 oscillatory network topologies (Figure 2E). In Figure 2E, each dot represents an oscillatory network topology, with the xaxis representing the ranking according to ${T}_{opt}$ (left panel) or ${A}_{opt}$ (right panel). The Spearman coefficient between ${T}_{opt}$ and 90 percentile of $\tau /T$ for all 474 oscillatory network topologies is 0.94, which is larger than that between ${A}_{opt}$ and 90 percentile of $\tau /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}_{opt}$ for C1, C2, …, and C5 categories (Figure 2F). The ranking of these five categories according to ${T}_{opt}$ is C1 > C2 ≈ C3 > C4 > C5, which is obtained by the onetailed Wilcoxon ranksum 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 activatorinhibitor, 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 ($\tau /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 activatorinhibitor 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).
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 $\tau /T$ and ${T}_{opt}$ for all 474 oscillatory network topologies is 0.72, which is smaller than that of 90 percentiles of $\tau /T$ and ${A}_{opt}$ (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 $\tau /T$ and ${A}_{opt}$ 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 ranksum tests to amplitude average (${A}_{opt}$) 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}_{opt}$ 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 activatorinhibitor, 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 righthand 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:
where $x\in {\mathbb{R}}^{3}$ , $f(\cdot ):{\mathbb{R}}^{3}\to {\mathbb{R}}^{3}$ , $B(\cdot ):{\mathbb{R}}^{3}\to {\mathbb{R}}^{3\times 3}$ and $\xi \left(t\right)\in {\mathbb{R}}^{3}$ is random noise. The unperturbed system $\dot{x}=f\left(x\right)$ has a periodic solution ${x}_{s}\left(t\right)$ (with period $T$). It can be proved that the variance of the phase deviation ${\sigma}^{2}(t)$ satisfies ${\sigma}^{2}(t)=ct$, and $c$ is as follows:
where ${v}_{1}^{T}\left(t\right)$ is the first row of the matrix $V\left(t\right)$. Here, the first column of ${V}^{1}\left(t\right)$ is ${\dot{x}}_{s}\left(t\right)$ , and ${V}^{1}\left(t\right)diag[{\mu}_{1},{\mu}_{2},{\mu}_{3}]V\left(t\right)$ is the state transition matrix for $\dot{w}=A\left(t\right)w\left(t\right)$ where ${\mu}_{i}$ ’s are Floquet exponents and $A\left(t\right)={\frac{\partial f\left(x\right)}{\partial x}}_{x={x}_{s}\left(t\right)}$ (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 $\tau /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 $MT$. Besides, ${v}_{1}\left(t\right)$ becomes $M{v}_{1}\left(\frac{t}{M}\right)$ (see ‘Methods’ for details). For the noise term, we merge $\frac{1}{M}$ with kinetic parameters ${v}_{ij},{\delta}_{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\left(x\right)$ becomes $B\left(x\right)\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\left(x\right)$ becomes $B\left(x\right)\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 $\stackrel{~}{A}/N$, $\stackrel{~}{B}/N$, and $\stackrel{~}{C}/N$, respectively (Figure 4B). This rescaling makes amplitudes of $\stackrel{~}{A}$$N$ times as high as that of $A$, and so do $\stackrel{~}{B}$ and $\stackrel{~}{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 $\stackrel{~}{A}$ , $\stackrel{~}{B}$ , and $\stackrel{~}{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 activatorinhibitor
In previous sections, we have used two and threenode 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 threenode 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 feedbackonly circuit 1 (Figure 5A) by integrating the synthetic RelAIκBα signaling circuit into the yeast MAPK pathway. The nucleartocytoplasmic 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 singlecell 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 singlecell time trajectories for these three circuits, we employ timelapse 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 activatorinhibitor topology.
Such improvement of the oscillation accuracy when adding a repressilator topology to the activatorinhibitor 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}_{CA}$ (the binding affinity of protein C to gene A): (near) zero $1/{K}_{CA}$ indicates the little effect of protein C to protein A, and thus protein A is not affected by protein C, leading to the activatorinhibitor; nonzero $1/{K}_{CA}$ implies the existence of the inhibition from protein C to protein A, resulting in the network with an activatorinhibitor and a repressilator. For a given oscillation parameter set for the activatorinhibitor, we first set $1/{K}_{CA}$ 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}_{CA}$ , 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 activatorinhibitor are usually smaller than those with an additional repressilator, and this gap enlarges with increased $1/{K}_{CA}$ , 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 activatorinhibitor 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 bottomup strategy (Lim et al., 2013; Zhang and Tang, 2019).
Here, we systematically explore the network design principles for accurate oscillation in both two and threenode 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 activatorinhibitor 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 negativefeedback 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 lowrobustness 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 timeaveraging 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 ‘‘switchlike’’ 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 feedforward 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 PotvinTrottier 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 toggleswitch system and excitable system (Lindner et al., 2004; Terebus et al., 2019; Zaks et al., 2005). The possible reason might be the noiseinduced transition between different states. Since our work only studied network topologies whose deterministic model can generate oscillation, we did not count the topologies that cannot oscillate in the deterministic model but begin to oscillate in the stochastic model. Due to the popularity of such topologies, how these topologies buffer noise will be of interest and may lead to the discoveries of new principles.
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; VelizCuba 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 Hilltype rate functions is the low computation cost, and in some biological networks, it leads to consistent results with the model composed of all elementary reactions (Gonze et al., 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 Hilltype reaction rate is based on the quasisteadystate 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 Hilltype form is decomposed into the elementary reactions. The recent development about stochastic algorithms with fast computation makes it feasible to simulate such detailed model for all two and threenode network topologies, for example, algorithms focusing on solving the chemical master equations (Cao et al., 2010; Cao et al., 2016; Munsky and Khammash, 2006; Terebus et al., 2021) and variants of 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 twonode and threenode 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 threenode transcriptional regulatory network:
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}_{ij}$ represents the maximal production rate of protein $i$ regulated by protein $j$, with ${K}_{ij}$ binding affinity. If there exist proteins activating gene i, ${\delta}_{i}$ is zero; if no protein activates gene i, ${\delta}_{i}$ is nonzero and represents the production rate caused by other proteins. ${k}_{basal}$ is the basal production rate. The equations for the twonode 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:
where ${G}_{B},A,B$ denote gene B, protein A, and protein B, respectively. The dissociation rates (${k}_{\mathrm{r}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{s}\mathrm{e}}/{k}_{\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{w}\mathrm{a}\mathrm{r}\mathrm{d}}$) for these six reactions are ${K}_{1},{K}_{2},\cdots ,{K}_{6}$ . Therefore, the fraction of the gene B at the state ${G}_{B}\bullet {B}_{3}$ in equilibrium is given by
Furthermore, we assumed that (2) ${K}_{6}\ll {K}_{5},{K}_{4}$ and ${K}_{3}\ll {K}_{2},{K}_{1}$ , so that this fraction can be rewritten as
We further assumed (3) that only gene B staying at the state ${G}_{B}\bullet {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
where ${v}_{BB}$ 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}_{AB}^{3}$ and ${K}_{BB}^{3}$ , respectively.
Stochastic model in the presence of extrinsic noise
To generate extrinsic noise, we perturb each kinetic parameter ${v}_{ij},{\delta}_{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:
Here, the control parameter $\epsilon $ indicates magnitude of perturbation of kinetic parameters, and large $\epsilon $ represents big perturbation of kinetic parameters. ${\eta}_{ij}$ , ${\xi}_{i}$, ${\zeta}_{i}$ are independent noise terms and all modeled by the Ornstein–Uhlenbeck process:
where ${W}_{t}$ is standard Wiener processes. This equation implies that $z\left(t\right)$ has zero mean and variance $\frac{{\mathrm{\sigma}}^{2}}{2{\tau}^{noise}}$ .
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:
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; VelizCuba 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_{A}X_{A}, and increasing protein A by 1 corresponds to ${k}_{basel}+V\left({\mathrm{\Sigma}}_{i}{v}_{iA}{\left(\frac{{X}_{i}}{{K}_{iA}}\right)}^{3}+{\delta}_{A}{V}^{3}\right)/\left({V}^{3}+{\mathrm{\Sigma}}_{i}{\left(\frac{{X}_{i}}{{K}_{iA}}\right)}^{3}+{\mathrm{\Sigma}}_{i}{\left(\frac{{Y}_{i}}{{K}_{iA}}\right)}^{3}\right)$. 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:
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 ${\eta}_{ij},{\xi}_{j},{\zeta}_{j}$ and used the Euler scheme to solve the dynamics of proteins’ concentrations. To be specific, the noise term $z$ ($z={\eta}_{ij},{\xi}_{j},{\zeta}_{j}$) at $n+1$ time step is determined by the following manner (${\tau}^{noise}=1$):
where $\u2206t$ is the time step, and ${\delta W}_{n}$ obeys the normal distribution with mean zero and variance $\u2206t$. Then, the protein’s concentration is solved by the Euler scheme (taking A as an example):
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:
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 threenode 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 $\tau /\mathrm{T}$ is $1/\mathrm{log}\left(c\right)$ , 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:
where $x\in {\mathbb{R}}^{3}$ , $f(\cdot ):{\mathbb{R}}^{3}\to {\mathbb{R}}^{3}$ , $B(\cdot ):{\mathbb{R}}^{3}\to {\mathbb{R}}^{3\times 3}$ and $\xi \left(t\right)\in {\mathbb{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\left(x\right)$ has a periodic solution ${x}_{s}\left(t\right)$ (with period $T$). Linearizing the noiseperturbed system around ${x}_{s}\left(t\right)$ gives the following system:
where $w\left(t\right)=x\left(t\right)\mathrm{}{x}_{s}\left(t\right)$ , $A\left(t\right)={\frac{\partial f\left(x\right)}{\partial x}}_{x={x}_{s}\left(t\right)}$ is $T$–periodic. From Floquet theory, the state transition matrix $\mathrm{\Phi}(t,s)$ for $\dot{w}=A\left(t\right)w\left(t\right)$ is given by
where $U\left(t\right)$ is $T$periodic nonsingular matrix with columns denoted by ${u}_{i}\left(t\right)$ , and $V\left(t\right)$ with rows denoted by ${v}_{i}^{T}\left(t\right)$ is ${U}^{1}\left(t\right)$. The ${\mu}_{i}$ ’s are Floquet exponents. We can further choose ${\mu}_{1}=0$ and ${u}_{1}\left(t\right)={\dot{x}}_{s}\left(t\right)$. Then corresponding ${v}_{1}\left(t\right)$ will play an important role in calculating the phase noise.
From the nonlinear perturbation analysis, $z\left(t\right)={x}_{s}\left(t+\alpha \left(t\right)\right)+y\left(t\right)$ solves Equation 5 for a small $y\left(t\right)$. The $\alpha \left(t\right)$ and $y\left(t\right)$ are called as phase noise and deviation noise, respectively. It can be proved that the variance of the phase noise $\alpha \left(t\right)$ increases linearly with time $t$, that is,
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 ${\mathit{v}}_{1}\left(\mathit{t}\right)$ when tuning $\mathit{M}$
We first consider the oscillator governed by the following equation:
We still use notations in the previous section to denote the quantities for this system. For example, the solution ${x}_{s}\left(t\right)$ , the period $\mathrm{T}$, the Jacobi $A\left(t\right)$ of $f\left(x\right)$ at the solution ${x}_{s}\left(t\right)$ . We assume the state transition matrix $\mathrm{\Phi}\left(t,s\right)=U\left(t\right)\mathrm{exp}\left(D\left(ts\right)\right)V\left(s\right)$ , which satisfies the first Floquet exponent is 0 and the first column of $U\left(t\right)$ is the time derivative of ${x}_{s}\left(t\right)$ . Then the first row of $V\left(t\right)$ is denoted as ${v}_{1}^{T}\left(t\right)$ , which can be used to calculate the variance of phase noise.
Next, we explore how ${v}_{1}\left(t\right)$ changes when the righthand term is divided by the timescale $M$. By this way, we obtain
It is easy to verify that the system governed by Equation 10 has a periodic solution ${x}_{s}\left(t/M\right)$ with period $MT$. The linearization of this system gives
where $w\left(t\right)=x\left(t\right)\mathrm{}{x}_{s}\left(t/M\right)$ . According to the definition of $\mathrm{\Phi}\left(t,s\right)$ , $\mathrm{\Phi}\left(t,s\right)$ satisfies
So $\mathrm{\Phi}\left(\frac{t}{M},\frac{s}{M}\right)$ satisfies
Therefore, $\mathrm{\Phi}\left(\frac{t}{M},\frac{s}{M}\right)$ is the state transition matrix for $\dot{w}=\frac{1}{M}A\left(\frac{t}{M}\right)w\left(t\right)$ . Since $\mathrm{\Phi}\left(t,s\right)=U\left(t\right)\mathrm{exp}\left(D\left(ts\right)\right)V\left(s\right)$ , we can obtain
where the first term $\frac{1}{\mathrm{M}}$ in the righthand term is to ensure the first column of $\frac{1}{M}U\left(\frac{t}{M}\right)$ is the time derivative of ${x}_{s}\left(\frac{t}{M}\right)$ , that is, $\frac{1}{M}{\dot{x}}_{s}\left(\frac{t}{M}\right)$ . Thus, the first row of $\mathrm{M}\mathrm{V}\left(\frac{\mathrm{t}}{\mathrm{M}}\right)$ is $M{v}_{1}^{T}\left(\frac{\mathrm{t}}{\mathrm{M}}\right)$ , 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 righthand side and perturb the ${v}_{ij},{\delta}_{i},{r}_{i}$ to introduce the extrinsic noise, thus leading to the following equations:
For simplicity, we still use $x$ and $f\left(x\right)$ to denote $\left(A,B,C\right)$ and the terms in the first brackets in the righthand terms, respectively. Thus, the above equations can be rewritten as
where $\Pi \left(t\right)$=$({\eta}_{AA},{\eta}_{BA},{\eta}_{CA,}{\xi}_{A},{\zeta}_{A},{\eta}_{AB},{\eta}_{BB},{\eta}_{CB,}{\xi}_{B},{\zeta}_{B},{\eta}_{AC},{\eta}_{BC},{\eta}_{CC,}{\xi}_{C},{\zeta}_{C})$, and ${B}_{ex}\left(x\right)$ is the matrix whose elements are coefficients of the random noise when $M=1$. Recall that the ‘ ${v}_{1}\left(\mathrm{t}\right)$ ’ in Equation 8 is $M{v}_{1}\left(\frac{\mathrm{t}}{\mathrm{M}}\right)$ for the system $\dot{x}=f\left(x\right)\frac{1}{M}$ , so the slope of variance of phase noise over period for the system with timescale $M$ is
By replacing ${t}^{`}$ with ${t}^{``}M$ and using ${t}^{`}$ to denote ${t}^{``}$ , we obtain
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 righthand side and introduce the cell volume $V$ to incorporate the intrinsic noise, thus leading to the following chemical Langevin equations:
We use $X$ and $f\left(X\right)$ to denote $\left({X}_{A},{X}_{B},{X}_{C}\right)$ and the terms in the first brackets in the righthand terms, respectively. Thus, the above equations can be rewritten as
where $\mathrm{\Lambda}\left(t\right)$ is $(d{W}_{A}/dt,d{W}_{B}/dt,d{W}_{C}/dt)$, and ${B}_{in}\left(X\right)$ is the matrix whose elements are coefficients of the random noise when $M=1$. If we use ${v}_{1}^{in}\left(\mathrm{t}\right)$ to represent the ‘ ${v}_{1}\left(\mathrm{t}\right)$ ’ for the system $\dot{X}=f\left(X\right)$ , the ‘ ${v}_{1}\left(\mathrm{t}\right)$ ’ for the system $\dot{X}=f\left(X\right)\frac{1}{M}$ is $M{v}_{1}^{in}\left(\frac{\mathrm{t}}{\mathrm{M}}\right)$ . So the slope of variance of phase noise over period for the system with cell volume $V$ is
By replacing ${t}^{`}$ with ${t}^{``}M$ and using ${t}^{`}$ to denote ${t}^{``}$ , we obtain
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 $\mathit{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 $\left(\stackrel{~}{A}/N,\stackrel{~}{B}/N,\stackrel{~}{C}/N\right)$ , which allows us to tune the amplitude by varying $N$. After this rescaling, we obtain the following equations for $\stackrel{~}{A}$ , $\stackrel{~}{B}$ and $\stackrel{~}{C}$:
where ${\stackrel{~}{x}}_{i}=\frac{{x}_{i}}{N},\mathrm{}{\stackrel{~}{y}}_{i}=\frac{{y}_{i}}{N}$ . If $N=1$, this equation is the same as Equation 2, so $\stackrel{~}{A}$ , $\stackrel{~}{B}$, and $\stackrel{~}{C}$ show same amplitudes with $A,B$, and $C$, respectively. Nevertheless, if $N\ne 1$, the amplitude of $\stackrel{~}{A}$ , $\stackrel{~}{B}$, or $\stackrel{~}{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}_{ij}$ ’s, and ${\delta}_{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}_{ij},{\delta}_{i},{r}_{i}$ by the same method as mentioned in the section ‘Mathematical modeling’ and obtain a new system described by the following equations:
where ${\eta}_{ij}$ , ${\xi}_{i}$, ${\zeta}_{i}$ are independent noise terms and all modeled by Equation 4 By multiplying $1/N$ to both sides of above equations, we get
Let $\stackrel{~}{\stackrel{~}{A}}=\frac{\stackrel{~}{A}}{N},$$\stackrel{~}{\stackrel{~}{B}}=\frac{\stackrel{~}{B}}{N}$ , and $\stackrel{~}{\stackrel{~}{B}}=\frac{\stackrel{~}{B}}{N}$ , the set of equations for $\stackrel{~}{\stackrel{~}{A}}$ , $\stackrel{~}{\stackrel{~}{B}}$ , and $\stackrel{~}{\stackrel{~}{C}}$ is the same as that in Equation 3, in which $N$ do not appear. So the dynamics of $\stackrel{~}{\stackrel{~}{A}}$ , $\stackrel{~}{\stackrel{~}{B}}$ , or $\stackrel{~}{\stackrel{~}{C}}$ will not change with varied $N$, thus leading to the same accuracy of oscillation when varying $N$. Based on $\stackrel{~}{\stackrel{~}{A}}N=\stackrel{~}{A},$$\stackrel{~}{\stackrel{~}{B}}N=\stackrel{~}{B}$, and $\stackrel{~}{\stackrel{~}{C}}N=\stackrel{~}{C}$ and the fact that the rescaling has no effect on the correlation function, $\stackrel{~}{A}$ shows the same accuracy of oscillation as $\stackrel{~}{\stackrel{~}{A}}$ , and so does $\stackrel{~}{B}$ or $\stackrel{~}{C}$ . Therefore, in the system for $\stackrel{~}{A}$ , $\stackrel{~}{B}$, and $\stackrel{~}{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
where ${X}_{\stackrel{~}{A}}=V\stackrel{~}{A}$ , ${X}_{\stackrel{~}{B}}=V\stackrel{~}{B}$, and ${X}_{\stackrel{~}{C}}=V\stackrel{~}{C}$ , and $V$ is the cell volume. By multiplying $1/N$ to both sides of above equations, we get
Let ${X}_{\stackrel{~}{\stackrel{~}{A}}}=\frac{{X}_{\stackrel{~}{A}}}{N},$${X}_{\stackrel{~}{\stackrel{~}{B}}}=\frac{{X}_{\stackrel{~}{B}}}{N},$ and ${X}_{\stackrel{~}{\stackrel{~}{C}}}=\frac{{X}_{\stackrel{~}{C}}}{N}$ , equations for ${X}_{\stackrel{~}{\stackrel{~}{A}}}$ , ${X}_{\stackrel{~}{\stackrel{~}{B}}}$ , and ${X}_{\stackrel{~}{\stackrel{~}{C}}}$ are
In above equations, $N$ only negatively affects the magnitude of noise term, so the oscillation accuracies of ${X}_{\stackrel{~}{\stackrel{~}{A}}}$ , ${X}_{\stackrel{~}{\stackrel{~}{B}}}$ , and ${X}_{\stackrel{~}{\stackrel{~}{C}}}$ increased with increased $N$. Thus, the oscillation accuracies of ${X}_{\stackrel{~}{A}}$ , ${X}_{\stackrel{~}{B}}$ , and ${X}_{\stackrel{~}{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

The free energy cost of accurate biochemical oscillationsNature Physics 11:772–778.https://doi.org/10.1038/nphys3412

Accurate chemical master equation solution using multifinite buffersMultiscale Modeling & Simulation 14:923–963.https://doi.org/10.1137/15M1034180

Phase noise in oscillators: a unifying theory and numerical methods for characterizationIEEE Transactions on Circuits and Systems I 47:655–674.https://doi.org/10.1109/81.847872

Noiseinduced multistability in chemical systems: discrete versus continuum modelingPhysical Review. E, Statistical, Nonlinear, and Soft Matter Physics 91:042111.https://doi.org/10.1103/PhysRevE.91.042111

Stochastic gene expression in a single cellScience 297:1183–1186.https://doi.org/10.1126/science.1070919

Improved leapsize selection for accelerated stochastic simulationJ Chem Phys 119:8229–8234.https://doi.org/10.1063/1.1613254

Stable oscillations in mathematical models of biological control systemsJournal of Mathematical Biology 6:207–223.https://doi.org/10.1007/BF02547797

Deterministic versus stochastic models for circadian rhythmsJournal of Biological Physics 28:637–653.https://doi.org/10.1023/A:1021286607354

Signaldependent dynamics of transcription factor translocation controls gene expressionNature Structural & Molecular Biology 19:31–39.https://doi.org/10.1038/nsmb.2192

BookNumerical Solution of Stochastic Differential EquationsBerlin: SpringerVerlag.

Effects of noise in excitable systemsPhysics Reports 392:321–424.https://doi.org/10.1016/j.physrep.2003.10.015

Delayinduced degradeandfire oscillations in small genetic circuitsPhysical Review Letters 102:068105.https://doi.org/10.1103/PhysRevLett.102.068105

Robustness of clocks to input noisePhysical Review Letters 121:078101.https://doi.org/10.1103/PhysRevLett.121.078101

Noise control and utility: from regulatory network to spatial patterningScience China Mathematics 63:425–440.https://doi.org/10.1007/s1142501916331

Design principles of biochemical oscillatorsNature Reviews. Molecular Cell Biology 9:981–991.https://doi.org/10.1038/nrm2530

Molecular architecture of the mammalian circadian clockTrends in Cell Biology 24:90–99.https://doi.org/10.1016/j.tcb.2013.07.002

Noise filtering strategies in adaptive biochemical signaling networks: application to E. coli chemotaxisJournal of Statistical Physics 142:1206–1217.https://doi.org/10.1007/s109550110169z

Adaptation with transcriptional regulationScientific Reports 7:42648.https://doi.org/10.1038/srep42648

A synthetic lowfrequency mammalian oscillatorNucleic Acids Research 38:2702–2711.https://doi.org/10.1093/nar/gkq121

Sources of variability in a synthetic gene oscillatorPLOS Computational Biology 11:e1004674.https://doi.org/10.1371/journal.pcbi.1004674

Stochastic dynamics and nonequilibrium thermodynamics of a bistable chemical system: the schlögl model revisitedJournal of the Royal Society, Interface 6:925–940.https://doi.org/10.1098/rsif.2008.0476

A critical quantity for noise attenuation in feedback systemsPLOS Computational Biology 6:e1000764.https://doi.org/10.1371/journal.pcbi.1000764

Bifunctional biochemical networksPhysical Biology 16:016001.https://doi.org/10.1088/14783975/aae74c
Decision letter

Naama BarkaiSenior and Reviewing Editor; Weizmann Institute of Science, Israel
Our editorial process produces two outputs: (i) public reviews designed to be posted alongside the preprint for the benefit of readers; (ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.
Decision letter after peer review:
Thank you for submitting your article "Network design principle for robust oscillatory behaviors with respect to biological noise" for consideration by eLife. Your article has been reviewed by 2 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Naama Barkai as the Senior Editor. The reviewers have opted to remain anonymous.
The reviewers found the work interesting, but reviewer #2 raised some issues that need to be addressed prior to publication. please address all these issues below.
Reviewer #1 (Recommendations for the authors):
The authors may consider to discuss further how the regulation changes under the fixed topology influence the oscillation quality compared with those under different topologies and see which one causes more significant changes.
Reviewer #2 (Recommendations for the authors):
The initial analysis was based on oscillation behavior of ODE models at different parameter values. There have been detailed studies in the past with this approach, and the current results appear to be similar, eg. oscillation appear to coincide with the graph cycle structures in the network. The authors may wish to discuss their findings from ODE analysis with previous results of 'Glass, Leon, and Joel S. Pasternack. "Stable oscillations in mathematical models of biological control systems." Journal of Mathematical Biology 6.3 (1978): 207223.’ and explain the relationship of their findings with these earlier results.
https://doi.org/10.7554/eLife.76188.sa1Author response
Reviewer #1 (Recommendations for the authors):
The authors may consider to discuss further how the regulation changes under the fixed topology influence the oscillation quality compared with those under different topologies and see which one causes more significant changes.
Thank you for the good suggestion. In the revision, we made two new figures (Figure 2—figure supplement 2 and Figure 3—figure supplement 4) and added the following paragraph in the Discussion:
“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 and Figure 3—figure supplement 4): the network topology with high robustness tends to show high dimensionless autocorrelation time when varying regulatory strengths, i.e., accurate oscillation (first ten 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 ten 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 ten 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 lowrobustness ones and then tune the regulation strength, as illustrated in Figure 5C.”
Reviewer #2 (Recommendations for the authors):
The initial analysis was based on oscillation behavior of ODE models at different parameter values. There have been detailed studies in the past with this approach, and the current results appear to be similar, eg. oscillation appear to coincide with the graph cycle structures in the network. The authors may wish to discuss their findings from ODE analysis with previous results of 'Glass, Leon, and Joel S. Pasternack. "Stable oscillations in mathematical models of biological control systems." Journal of Mathematical Biology 6.3 (1978): 207223.’ and explain the relationship of their findings with these earlier results.
Thank you for the good suggestion. The paper you mentioned modeled the feedback inhibition biological network by piecewise linear equations and derived conditions to generate oscillation. The core motif to produce oscillation is the same as that in our manuscript, and we pointed this out in Page 8 by adding the following text: “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).
https://doi.org/10.7554/eLife.76188.sa2Article and author information
Author details
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.
Senior and Reviewing Editor
 Naama Barkai, Weizmann Institute of Science, Israel
Publication history
 Received: December 7, 2021
 Preprint posted: December 23, 2021 (view preprint)
 Accepted: August 22, 2022
 Version of Record published: September 20, 2022 (version 1)
Copyright
© 2022, Qiao, Zhang, Zhao et al.
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.
Metrics

 355
 Page views

 140
 Downloads

 0
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Computational and Systems Biology
 Neuroscience
The locus coeruleus (LC) houses the vast majority of noradrenergic neurons in the brain and regulates many fundamental functions including fight and flight response, attention control, and sleep/wake cycles. While efferent projections of the LC have been extensively investigated, little is known about its local circuit organization. Here, we performed largescale multipatch recordings of noradrenergic neurons in adult mouse LC to profile their morphoelectric properties while simultaneously examining their interactions. LC noradrenergic neurons are diverse and could be classified into two major morphoelectric types. While fast excitatory synaptic transmission among LC noradrenergic neurons was not observed in our preparation, these mature LC neurons connected via gap junction at a rate similar to their early developmental stage and comparable to other brain regions. Most electrical connections form between dendrites and are restricted to narrowly spaced pairs or small clusters of neurons of the same type. In addition, more than two electrically coupled cell pairs were often identified across a cohort of neurons from individual multicell recording sets that followed a chainlike organizational pattern. The assembly of LC noradrenergic neurons thus follows a spatial and cell typespecific wiring principle that may be imposed by a unique chainlike rule.

 Computational and Systems Biology
Computational models starting from large ensembles of evolutionarily related protein sequences capture a representation of protein families and learn constraints associated to protein structure and function. They thus open the possibility for generating novel sequences belonging to protein families. Protein language models trained on multiple sequence alignments, such as MSA Transformer, are highly attractive candidates to this end. We propose and test an iterative method that directly employs the masked language modeling objective to generate sequences using MSA Transformer. We demonstrate that the resulting sequences score as well as natural sequences, for homology, coevolution and structurebased measures. For large protein families, our synthetic sequences have similar or better properties compared to sequences generated by Potts models, including experimentallyvalidated ones. Moreover, for small protein families, our generation method based on MSA Transformer outperforms Potts models. Our method also more accurately reproduces the higherorder statistics and the distribution of sequences in sequence space of natural data than Potts models. MSA Transformer is thus a strong candidate for protein sequence generation and protein design.