{'fraction_non_alphanumeric': 0.0264572808245582, 'fraction_numerical': 0.00720271065392497, 'mean_word_length': 4.856463878326996, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Design principles of noise-induced oscillation in biochemical reaction networks: II. coupled positive and negative feedback loops

According to the chemical reaction network theory, the topology of a certain class of chemical reaction networks, regardless of the kinetic details, sets a limit on the dynamical properties that a particular network can potentially admit; the structure of a network predetermines the dynamic capacity of the network. We note that stochastic fluctuations can possibly confer a new dynamical capability to a network. Thus, it is of tremendous value to understand and be able to control the landscape of stochastic dynamical behaviors of a biochemical reaction network as a function of network architecture. Here we investigate such a case where stochastic fluctuations can give rise to the new capability of noise-induced oscillation in a subset of biochemical reaction networks, the networks with only three biochemical species whose reactions are governed by mass action kinetics and with the coupling of positive and negative feedback loops. We model the networks with the master equations and approximate them, using the linear noise approximation. For each network, we read the signal-to-noise ratio value, an indicator of amplified and coherent noise-induced oscillation, off from the analytically derived power spectra. We classify the networks into three performance groups based on the average values of the signal-to-noise ratio and the robustness. We identify the common network architecture among the networks belonging to the same performance group, from which we learn that the coupling of negative and positive feedback loops generally enhance the noise-induced oscillation performance better than the negative feedback loops alone. The performance of networks also depends on the relative size of the positive and negative feedback loops; the networks with the bigger positive and smaller negative feedbacks are much worse oscillators than the networks with only negative feedback loops.

prevalent in cellular biology partly because of the living organisms' adaptation or readiness to a periodic or abrupt environmental change. The few examples of the cellular biochemical oscillations are cell cycle, circadian rhythm, NF-κB, p53, developmental clock, neural rhythms, and hormone [12,13]. To quantitatively understand such periodic phenomena, many theoretical models have been put forth and extensively studied. The previous studies revealed a few important requirements for biochemical oscillators [14]: i) negative feedback loops with sufficient time delay, ii) non-linearity of kinetic laws, and iii) appropriate balancing of synthesis and degradation rates of chemical species. One of the above three restrictions is the structral condition of the underlying networks: negative feedback loops with time delay. A negative feedback loop is the cyclic pathway consisting of the odd number of inhibitory edges. The time delay can be realized in the networks in two different manners: One by having the explicit time delay in biochemical interaction and another by having positive feedback loops in addition to the negative feedback loops. Our primary questions in this paper are how the noise relaxes or tightens the above three conditions for biochemical oscillations. Particularly, we are very keen to enact the new requirements for stochastic biochemical oscillators and to compare the stochastic oscillatory behaviors between networks with negative feedback loops alone and the networks with coupled positive and negative feedback loops.
In our previous work [15], we investigated the requirements for stochastic biochemical oscillators and compared the networks with only negative feedback loops. The networks consist of three biochemical species governed by mass action kinetics and are allowed to have only negative feedback loops. We proved that the negative feedback loops are required for stochastic oscillation in these small-sized networks with mass action kinetics. Then, we numerically demonstrated that all the networks have one positive and locally stable steady state and the stochastic fluctuations enable all of them to exhibit prominent stochastic oscillations, i.e., coherence resonance or noise-induced oscillation, in the various biologically feasible parameter ranges [16][17][18][19][20][21]. Stochastic fluctuations confer an additional dynamical capacity of stochastic oscillation to the group of biochemical reaction networks with negative feedback loops.
The few examples are mitiotic trigger in Xenopus [26], Ca 2+ spikes/oscillations [27], circadian clock [28], galactose-signaling network in yeast [29], and p53-Mdm2 oscillator [30]. The addition of positive feedback loops to the biochemical oscillators can confer some functional and performance advantages [22][23][24][25]. The networks with the positive-negative interlinked feedbacks can have the inreased frequency-tunable range and the enhanced robustness of biochemical oscillations [22]. Additionally, the positive-negative interlinked feedback loops can give rise to a variety of dynamical behaviors such as monostability, bistability, exitability, and oscillations in the change of relative strength of the feedback loops [23] and even make the cellular signal responses more desirable in noise environment [25]. It seems that the nature must have forced the living organisms to evolve to favor the networks with interlinked postive and negative feedback loops because of their many functional advantages.
In this present paper, we add one more favorable point to the long list of the advantage of networks with positive and negative interlinked feedbacks: networks with positive-negative coupled feedbacks are much better noise-induced oscillators than networks with only negative feedback loops in various biologically feasible parameter ranges. Again, we consider small-sized networks with mass action kinetics and two types of repression: repression by proteolysis and repression by transcriptional control. We generate the exhaustive list of all possible network structures with interlnked feedback loops, resulting in the sixty-three networks with the differential coupling of positive and negative feedback loops, and model each network with a chemical master equation which is approximated to a linear Fokker-Planck equation through van Kampen's system size expansion. Firstly, stochastic fluctuations indeed confer a new dynamical capacity of stochastic oscillation to all the networks with interlinked feedback loops. Secondly, implementing a Kmedoids clustering algorithm, we group the sixty-three networks into three performance groups based on the average values of signal-to-noise ratio and robustness and identify the common network architecture among the networks belonging to the same performance group. We learn that the coupling of negative and positive feedback loops (PNFBL) generally enhance the noise-induced oscillation performance better than the negative feedback loops (NFBL) alone. However, the performance of PNFBL networks depends on the size of the positive feedback loops (PFBL) relative to that of the NFBL in the networks; the performance of the networks with the bigger PFBL than NFBL is worse than that of the networks with only NFBL. We single out two networks which stand out in their performance of noise-induced oscillation in all four different parameter ranges and for two different models of repression. Thirdly, we also elucidate the machanisms of noise-induced oscillation and find that for most networks, the natural frequency sets the low bound of the resonant frequency of the noise-induced oscillation. This corroborates the present understanding of the noise-induced oscillation mechanisms; the imaginary part of a complex eigenvalue from the Jacobian matrix of a linearized network is responsible to generate an inherent rotation in the network and the noise enhances such a rotation. However, we also find that a few networks with purely real eigenvalues can admit a very prominent and amplified noise-induced oscillation.
The paper is organized in the following way. In the result section, we discuss the K-medoids clustering of sixty-three networks into three performance groups, the classification of networks by network structure, the identification of a common network structural signature per performance group, the distributions of signal-to-noise ratio from real and complex eigenvalues of Jacobian matrices, and finally the correlation between the values of signal-to-noise ratio and the spiralness and proximity. In model and method section introduced all the networks under our consideration, the derivation of the power spectrum from the chemical master equation based on the van Kampen system size expansion, parameter sampling method, definition of signal-to-noise ratio, prominence and robustness, the eigenvalue calculation, and the linear stability of the networks.   prominence robustness FIG. 1: Classification of sixty-three biochemical reaction networks into three performance groups, based on the values of prominence and robustness and using K-medrid clustering algorithm. In top row (a)-(d), the repression is modeled by proteolysis. The proteolysis rate is sampled from the following four intervals: a) (10 −7 , 10 −5 ), b) (10 −5 , 10 −3 ), c) (10 −3 , 10 −1 ) and d) (10 −7 , 10 −1 ). In bottom row, repression is modeled by transcriptional control. The transcriptional repression rate is sampled from the following four intervals: e) (10 −3 , 10 −1 ), f) (10 −1 , 10), g) (10, 10 3 ), and h) (10 −3 , 10 3 ). The values of all the other kinetic rate constants are sampled from a preset respective biologically feasible interval for a)-c) and e)-g) whereas all of them are sampled from the same interval for d) and h).

10
The exhaustive enumeration of all the networks with three nodes and with at least one negative feedback loop reveals that there are sixty-three biochemical reaction networks. Each of sixty-three networks has one positive and locally stable steady state in the chosen parameter ranges. If a network admits an unstable steady state, then that set of parameter values is resampled until the network admits the stable fixed point. The statics of resampling cases is provided in the supplimentary information. Sixty-three biochemical reaction networks are classified, based on their values of prominence and robustness. Prominence is a measure of the coherence and amplification of a stochastic oscillation, defined as the averaged maximum signal-to-noise ratio (SNR). Robustness is the fraction of the sample points that admit the value of the maximum SNR greater than one. The maximum SNR means that the largest SNR value among three chemical species. See the method section for the detailed definition of prominence and robustnss and how to calculate them. To classify the networks, we use a machine-learning classification tool, K-medoids clustering. The kinetic rate constant values are sampled in two different manner, biologically and non-biologically. For biologically feasible parameter values, all the kinetic rate constants are sampled from a preset biologically feasible range in Table  IV. Particularly, we consider three different intervals of the repression strength corresponding to weak, intermediate, and strong repression. We do not know even the approximate values of repression strength and want to see the dependence of noise-induced oscillation performance on repression strength. For non-biological parameter values, all the kinetic rate constants are sampled from the same "non-biological" range. This is to see the overall behavior of the networks across a parameter space. Fig. 1 shows the classification of the sixty-three networks into three performance groups. The performance of networks is based on the values of prominece and robustness of stochastic oscillation. We classify the biochemical reaction networks into three clusters, using K-medoids clustering. The K-medroid performs robustly and the choice of K=3 is well justified by the occurrence of a dramatic drop at K=3 in the graph of sum distance (error) versus the numebr of clusters over almost all parameter ranges. The graph of sum distance versus the number of clusters for Fig. 1 is provided in the Fig. V in the supplimentary information. Just like other searching algorithms, the K-metroid clustering algorithm can be trapped in local mimina. To prevent such a trap, we start from many different initial conditions and optimize the performance of the K-metriod clustering algorithm. Because of techanical issues, we could not evaluate one network (network 38) for the case of repression by proteolysis and five networks (networks 36-39, 61) for the case of repression by transcriptional control. We label three clusters of networks as follows: red as the best performing networks, green as the medium performers, and blue as the worst performers. In Fig. 2 and 3, we will indentify the common network structural properties among the networks belonging the same color group. The clustering also enables us to identify the all-season best performer of stochastic oscillation, network 31. This network 31 belong to the red group for both repression models and across four different sampling intervals of repression strength. The Network 31 consists of two three-dimensional NFBL plus three two-dimensional PFBL. One interesting observation is that the clustering is mainly done by the prominence values and does not strongly depend on the robustness. We do notice the decreasing trend of the prominence of most of the networks as the repression strength increases for both repression models in Fig. 1. For example, As the repression strength increases from a) to c), the prominence value of network 31 decreases respectively. This pattern is most conspicuously noticed with the networks in green performance group.
In Fig. 2, the sixty-three networks are classified into eight different architectural classes by network structural properties as provided in the Tables I and II. This newtork structural classification is based on the number and the size of positive and negative feedback loops. Those eight architectural classes of the sixty-three networks can form three architectural groups in a coarse-grained manner. The architectural group I includes all the networks with only negative feedback loops, Networks 1-13, whereas the architectural group II (classes II through IV) contains the networks with the coupling of smaller postive feedback loops and larger negative feedback loops, Networks 14-39. Finally the architectural group III (classes V through VIII) consists of the networks with the coupling of larger positive feedback loops and smaller negative feedback loops and the networks with linear chains, Networks 40-63. The networks in the architectural group I dominantly belong to blue performance group for the case of repression by proteolysis whereas for repression by transcription control all three colors are mixed. A majority of the networks in the architectural group II belong to green performance group for repression by proteolysis whereas for transcriptional repression, they belong to either green or red performance group with exception of a very few blue performance networks. The networks in the architectural group III, for both cases of repression by proteolysis and transcription control, mostly belong to blue or black performance group. As far as the case of repression by proteolysis is concerned, the blue performance arises from two architectural groups (or four architectural classes): group I (class I) and group III (classes V through VIII). The clear message is that the networks either with exclusively negative feedback loops or with the larger-sized positive feedback loops and the smaller-sized negative feedback loops can admit the noiseinduced oscillation in the chosen biologically feasible parameter ranges, but their oscillations are neither sufficiently well amplified nor coherent. In other words, the networks with the coupling of smaller-sized positive feedback loops and larger-sized negative feedback loops can admit quite well amplified and coherent stochastic oscillations. Thus, the recommended network strucure for biochemical oscillators is one of the networks belonging to the architectural group II.
To see the clear relationship among the individual network structures, the parameter sampling range, and the performance of networks altogether, we represent each network not with the absolute value of average max SNR but with the performance group to which it belongs. Then, we average the performance colors of each network over the four different parameter sampling ranges for two different repression models, as presented in the rightmost subfigure of Fig. 2. Since the average performance scores are quite noisy, it is very hard to appreciate any correlation between the individual networks and their average performance scores. But, it is easy to see that the average performance scores and the network architectural groups (represented with three different background darknesses in Fig. 2(c) are closely correlated for both repression models (blue and red lines in Fig. 2(c)). In other words, the networks belonging to their respective architectural group demonstrate the similar performance, regardless of the detailed repression models.
In Fig. 3, we identify the common architectural properties among the networks that are classified into the same performance group as discussed in Figure 1. Each of the networks belonging to the same performance group are rotated until they are all aligned such that the Hamming distance among the networks are minimized with respect to a reference network. We then calculate the average edge value, defined as η ij where ij denotes the directed  Tables I and II. While color indicates the absence of data due to technical issues. Black color denotes that all sample points yield SN R ≤ 1. Red, green or blue color is assigned to each network in each parameter sampling range and exactly corresponds to the classification from the K-medoids clustering in Figure 1. In the right graph is the average performance value for each network: integer numbers, 0,1,2,3 are assigned to four colors, black, blue, green, and red, respectively. Red line indicates the average performance value from (a) repression by proteolysis whereas blue line denotes the performace value from (b) repression by transcriptional control.
edge from a node X to a node Y as presented in Fig On the one hand, the networks belonging to the blue performance group share the following structural properties: one three-dimensional positive feedback loop coupled with many two-dimensional negative feedback loops. In other words, The blue networks are characerized with the strong presence of multiple two-dimensional (small) negative feedback loops and three-dimensional (large) positive feedback loop. We conclude that the networks with a larger positive feedback loop coupled with small negative feedback loops are likely to belong to the blue (worst) performance group.
On the other hand, the networks belonging to the green performance group are characterized by the strong presence of three-dimensional (large) negative feedback loop and one two-dimensional (small) positive feedback loop. In the networks belonging to the Red performance group, we see the same strong presence of three-dimensional negative feedback loop as wellas the stronger presence of multiple two-dimensional positive feedback loops. The best network topology for highly amplified and coherent stochastic oscillators is a large negative feedback loop coupled with many small positive feedback loops. We demonstrate the existence of the network signature for each performance group consistently across two different repression models and over four different intervals of the model parameters. This is a remarkable numerical evidence of the relationship between network architecture and network performance. In Fig. 4, we present the distributions of the resonant frequencies of stochastic oscillations from the sixty-three networks which are colored according to their performance group as defined in Fig. 1. The resonant frequency is defined as a non-zero frequency at which the power spectrum peaks. The distribution is plotted with a subset of power spectra which has a peak at non-zero frequency, i.e., a fraction of 10 6 sample points that yield maximum SNR≥ 1. Since any noise-free network has a single stable fixed point, the existence of a resonance frequency is purely due to noise-induced effect. For the case of repression by transcriptional control (e)-(h), the distributions of resonant frequencies are substantially overlapped over almost all of the networks, independent of their performance and network architectures. For the case of repression by proteolysis (a)-(d), the distributions are not as homogeneous as for the case of repression by transcriptional control. The blue-colored networks tend to have more homogeneous distribution than the green-colored networks. For both repression models, the peaks of resonant frequency distributions are shifted to the right as the repression strength increases. This is a clear signal that the noise-induced resonant frequency is positively correlated with the repression strength. The range of resonant frequencies falls within the biologically relevant range, between hours (corresponding to 10 −4 Hz) to years (corresponding to 10 −8 Hz).
In Fig. 5 are discussed the origins and mechanisms of the stochastic oscillations accompanied with the large maximum SNR values. Calculating the discriminant of a Jacobian matrix with a set of randomly sampled parameter values for a chosen network, we accurately determine if all three eigenvalues of the 3 x 3 Jacobian matrix are real or the mixture of real and complex conjugate pairs. We repeat the discrimimant calculation with the 10 6 sets of randomly sampled parameter values from four different parameter ranges for the sixty-three networks with two different repression models. Each subfigure in Fig. 5 presents three histograms of the maximum SNR values: one from the samples yielding all real eigenvalues, another from the samples resulting in complex conjugate pairs, and the last from the total samples. The distributions from only networks 31 and 44 are presentd in the main text, but the distributions from the rest of the networks are provided in the Figs. 9 and 10 in the supplimentary information. Note that the log of the maximum SNR values less than 0 do exist, but they are not presented here. For the case of network 31 and almost all of the networks in the network architectural group II (networks 14-39), the complex eigenvalues make a dominant contribution to the distribution of maximum SNR values. Particularly the high values of maximum SNR come exclusively from the samples yielding the complex eigenvalues. The imaginary part of the complex eigenvalues is indicative of a rotational flow in a deterministic system and the noise tends to amplify that rotational motion, resulting in the maximally amplified and coherent oscillation. So, this high correlation between the existence of complex eigenvalues and the high values of maximum SNR is well understood. The further discussion of the complex eigenvalue cases will follow in Fig. 6. A significant contribution of the real eigenvalues to the maxum SNR values is detected from network 44. All the networks belonging to the network architectural group III (networks 40-63) benefit the similar huge contribution from the real eigenvalues as shown in Fig. 9 in the supplimentary information. It is worthwhile to note that only a very few networks belonging to this third architectural group are in green performance group, but a majority of the networks are in blue performance group. This seems to suggest that the real eigenvalues do not elicit as large SNR values as the complex eigenvalues do. However, the network 44 belongs to the green performance group and its intermediate performance is exclusively due to purely the real eigenvalues, leaving us puzzled. Finally, the networks from the first network architectural group I (networks 1-13) have the very similar distributions from both real and complex eigenvalues. The above obervations are primarily based on the proteolysis repression model as shown in Fig. 9 in the supplimenary information. For the case of repression by transcriptional control, all of sixty-three networks dominantly produce the complex eigenvalues across four different parameter sampling ranges. The network architecture seems to determine the ratio between twos distributions from real and complex eigenvalues and to affect the dynamical behaviors accordingly.
The Fig. 6 shows that the resonant frequencies are positively correlated with the imaginary part of the complex eigenvalues, Imλ, of the Jacobian matrix of the linearized system. As seen in top panel of Fig. 6, the Imλ plays a role of the lower bound for the resonance frequency ω and the noise increases the values of resonant frequencies. As briefly discussed in value of Imλ is very large. However, as the value of Imλ gets smaller, the resonant frequency gets larger than the angular speed of the deterministic stable spiral. The noise seems to push the system along the deterministic rotational flow rather than to go against it. This tendency is consistently observed in three networks, networks 31, 35, and 59 and under two repression models as presented in the Fig. 11 in the supplimentary information. Also, as shown in the lower panel of Fig. 6, the stochastic oscillation is primarily driven by two causes: the complex eigenvalues of the Jacobian matrix being very close to the imaginary axes and having a larger imaginary part. Those two causes are quantified by two numeric values of proximity and spiralness. The maximum SNR values which are represented by color in heat maps are highly correlated with those two values of proximity as well as spiralness. Proximity is defined as the ratio of the magnitude of noise (B along an eigen-direction to the real negative part of the complex eigenvalues (Reλ) whereas spiralness is defined as the ratio of the angular speed of a deterministic stable spiral (Imλ) to the negative real part (Reλ) which can be thought as an attractive force to the fixed point . Thus the proximity measures the likelihood of making the fixed point stable and pushed the system away from its stable fixed point. The spiralness measures the coherence and amplification of the stochastic oscillations because Imλ is very closely related to the resonant frequency and Reλ is inversely proportional to the amplitude of the stochastic oscillations, i.e., the peak amplitude of a power spectrum. In the Fig. 12 in the supplimentary information, one representative network is selected from each performance group: network 59 for blue, network 35 for green, and network 31 for red performance group. The Fig. 12 clearly shows that for both repression models, the networks belonging to the better performance group have the larger values of proximity and spiralness overall and the larger values of maximum SNR.

III. DISCUSSION AND CONCLUSIONS
This paper is the first exensive comparative study pertaining to a stochastic dynamical behavior from stochastic biochemical reaction networks. Most importantly, we numerically demonstrate the strong correlation between the stochastic behavior and the network architectural properties, namely the coupling patterns of positive and negative feedback loops. We investigate noise-induced oscillation in the networks with only three biochemical species whose reactions are governed by mass action kinetics and with the coupling of positive and negative feedback loops. Modeling a set of so many stochastic biochemical reaction networks by using linear noise approximation and reading the signal-tonoise ratio values from the analytically derived power spectra, we show that all the networks with coupled positive and negative feedbacks are capable to admit the noise-induced oscillations in biologically feasible parameter ranges. Also, using a K-metroid clustering algorithm, we group the sixty-four networks into three performance groups and identify the common network architecture among the networks belonging to the same performance group. We learn that the coupling of negative and positive feedback loops (PNFBL) generally enhance the noise-induced oscillation performance better than the negative feedback loops (NFBL) alone. However, the performance of PNFBL networks depends on the size of the positive feedback loops (PFBL) relative to that of the NFBL in the networks; the performance of the networks with the bigger PFBL than NFBL is worse than that of the networks with only NFBL.
As shown in the table V in the supplimentary information, we realize that a few networks can generate an unstable fixed point in a very small fraction of the parameter space and this dynamical instability affects the performance of those networks tremendously. The noise-induced oscillation can arise with a significantly high signal-to-noise ratio when the dynamical instability is nearby, whose effect is known as the noisy precursors of the nonlinear instabilities ??. For the case of proteolysis, all the networks except three networks, networks 10, 31, and 55, have one locally stable fixed point across all parameter ranges. The networks 10, 31, and 55 generate an unstable fixed point with up to 9 pecentage of random sampling from each of four parameter ranges. The networks 10 and 31 outperform the other networks across all parameter ranges whereas the network 55 is in the blue performance group despite having the instabilities inside the parameter sampling ranges. However, the networks from the network architectural group II (networks 14-39) except network 31 don't exhibit any dynamical instability, but are still capable of generating the noise-induced oscillation with the extraordinarily large values of maximum SNR, i.e., up to maxSN R ∼ 10 1 0 in the biologically feasible ranges. For the case of transcriptional control, we find that the presence of instabilities of the networks within the chosen parameter range is nicely correlated with the performance of the networks almost without exception. The sample points leading to instabilities goes up to 9.5% for a few networks, but stay much less than 1%.

A. Biochemical reaction models
We consider biochemical reactions in which the chemical species get synthesised and degraded. The chemical species get synthesised constitutively and their synthesis can be enhanced by another chemical species, called activation. Most often, the activation occurs in a gene regulation where a protein functioning as a transcription factor binds to the promoter of a target gene and enhances the activity of the gene, increasing the production rate of the target protein. The chemical species get degraded spontaneously and they can be negatively regulated by another chemical species, called repression. In this paper, we consider two repression models: repression by proteolysis and repression by transcription control. Repression by proteolysis can be found in protein-protein leading to proteolysis such as ubiquitin-mediated protein-degradation. Repression by transcription control can be often found in gene reguation where a protein functioning as a repressor binds to a promoter of a target gene and shut off the activity of that gene, decreasing the synthesis of the target protein.
In this paper, we consider an ensemble of three-node directed graphs with the following constraints: (a) There can exist at most one directed edge from one node to another. (b) When a directed edge is present, it can be either inhibition or activation. In graphical representation, X → Y is for biochemical species X activating species Y , and X Y is for species X inhibiting species Y . (c) We only consider networks in which all three nodes are part of some loop. Thus networks that have nodes with only in-coming edges or out-going edges or are isolated are not considered. (d) Self-directed edges are not allowed. The exhaustive list of three-node directed graphs under our consideration is presented in the Tables I and II. Table I is for networks with only negative feedback loops while Table II is for networks with coupled positive and negative feedback loops. Only topologically distintive networks are allowed. If two networks are identical after rotation and/or mirroring, then the two networks are topologically identical. All the graphs can have at most 6 directed edges and the topology of each graph is determined by a distinctive arrangement of the 6 directed edges. Thus, we represent each graph with a string of 6 characters. Each of the 6 characters can be either "A" for activation, "I" inhibition, or "0" when an edge is absent. The first three characters indicate three directed edges going counterclockwise, one from species X to species Y , another from species Y to species Z, and the last one from species Z to species X. The last three characters indicate three directed edges going clockwise, one from species X to species Z, another from species Z to species Y , and then the last one from species Y to species X. For example see Fig. 7 for network number 7. This network consist of three counter-clockwise edges staring from X going to Y (activation), Y to Z (inhibition), and Z to X (activation). There are two clockwise edges in this network one going from X to Z (inhibition) and Y to X (inhibition), and there is no edge going from Z to Y . Therefore, we can represent network 7 with a string of 'AIAI0I'. We represent the rest of the networks in our ensemble by applying the same rules. The networks in our ensemble can be categorized into different groups depending on their underlying topological characteristics such as: 1) cyclic negative feedback loops (NFBLs); 2) linear negative feedback loops; 3) cyclic positivenegative interlinked feedback loops (PNFBLs); and, 4) linear positive-negative interlinked feedback loops. The linear feedback loops are made up of only two component feedback loops, which forms a linear chain of three nodes lying in a row. The cyclic feedback loops always involve three component feedbacks as a backbone structure, which form a triangular closed loop. In Table I, we classified networks 1-10 into type I and type II. Type I NFBLs have "AIA000" as the backbone structure while Type II NFBLs have "III000" as the backbone structure. All Type I and Type II NFBLs can be created by adding other links to the backbone structures "AIA000" and "III000", respectively. The cyclic PNFBLs can be constructed by adding positive feedback loops in an exhaustive way to the different types of the cyclic NFBL backbone networks. They are termed as the cyclic PNFBLs Type 1 and Type 2 when they are obtained by adding positive feedback loops to NFBLs Type 1 and Type 2, respectively.
In Tables I and II, we further indicate the number and types of feedback loops present in the networks in the form "nNpP", where n, p = 1, 2, 3. Here, n indicate the number of negative feedback loops (N) while p indicate the number of positive feedback loops (P).

B. Master Equations
In order to study the stochastic dynamics of a network, we start by writing down the chemical master equation that describes the time evolution of the states of the network of chemical reactants. We denote the state of the network as a vectorS = (X, Y, Z) that represent the current populations of the three chemical species of the network. We also denote the joint probability distribution of the state of the network at a particular time t as P (S; t) and the transition rate from stateS to stateS as T (S|S ). This allows us to formally write the master equation governing the time evolution of P (S; t) as To mathematically model the transition rates among the states, we make the following assumptions: (a) any given node is constitutively synthesized, which depends on the system size Ω, (b) any given node is degraded at a rate which depends on its current abundance, (c) activation is proportional to the abundance of the activating species, (d) for inhibitory interaction, we consider two models as discussed in Sec. IV A, repression by proteolysis and repression by transcription control; and (e) the corresponding stochastic chemical reaction is a Markovian process. Based on these assumptions, we can specify all the allowed transitions among different states. In Table III, we show the transition rates of synthesis and degradation associated with all the possible chemical reactions that can occur on a single node X. In Table III, we have shown the transition rates for both repression by proteolysis and repression by transcritpion control discussed in Sec. IV A. These transition rules combined with the particular network topology fully determine the mathematical formulation of the chemical master equations. The master equation can be rewritten by using the step operator notation in the following way:

Reactions
Transition rate where the step operators for species X, E + X and E − X , are defined by where f (X, Y, Z) is an analytic function. Operators E ± Y and E ± Z are similarly defined. Thus, the operation of the step operator is to either raise or lower the number of a node by 1.

C. System size expansion
The form of the chemical master equation does not allow the analytic solution in general due to the nonlinearity. One possible way to deal with this problem is to use the system size expansion method introduced by van Kampen. It allows us to expand the master equation systematically in terms of the system size Ω. Thereby, we can obtain the approximate solutions to the leading orders of the large system size.
Here, we sketch the standard procedure of the system size expansion briefly following Ref. [31]. Interested readers should refer to Ref. [31] for more detail. The main idea of the system size expansion assumes that the abundance of chemical species can be decomposed into its macroscopic mean value and the fluctuation around the mean, both of which depends on the system size Ω. For example, we can write the number of chemical species S i as where the index i runs over all the chemical species, i.e., S i = X, Y, Z. The first term Ωφ i describes the macroscopic mean value of S i . The second term Ω 1/2 ξ i is the gaussian fluctuations of S i around its mean. Then, we can rewrite the joint probability distribution P ({S i }; t) in terms of the set of stochastic variables of ξ i instead of S i . Upon such change of variables, the time derivative of the joint probability distribution transforms as Accordingly, the step operators are also rewritten in terms of the new variable.
Plugging Eq. (4) through Eq. (6) in Eq. (2), we can expand the master equation systematically in the order of the system size Ω. To the leading order of Ω expansion (Ω 1/2 ), we obtain a set of the macroscopic deterministic equations.
In the next leading order of Ω expansion (Ω 0 ), we obtain the linear Fokker-Plank equation.
where C(ξ i ) ≡ j J ij ξ j . J ij denotes the (i, j) th element of the Jacobian matrix J of the macroscopic rate equations and B ij denotes the (i, j) th element of the noise covariance matrix B.

D. Power spectrum
We can recast the linear Fokker-Plank equation into the mathematically equivalent form of Langevin equations.
where η i (t)η j (t ) = B ij δ(t−t ). The matrix elements J ij and B ij are evaluated at the fixed points of the macroscopic deterministic equations. In order to examine the existence of the stochastic ocillations, we can study the power spectrum of the stochastic variable ξ i (t). Namely, whereξ i (ω) denotes the Fourier transform of ξ i (t) and Φ ij ≡ −iωδ ij − J ij . For our case, J and B are 3 by 3 matrices. The explicit expression of the power spectrum for species i, in terms of matrix elements of J and B, is given by where the coefficients in the numerator are:

E. Sampling methods and root calculation
The networks we consider in this work are stable in the positive real domain in the parameter space, i.e., {k i } ∈ R + where {k i } are the kinetic rate constants shown in Table III. In order to demonstrate the existence of the stochastic amplified oscillations, we sample 10 6 sets of kinetic rate constants uniformly and randomly from the logarithmically scaled intervals by using the unconstrained Monte Carlo sampling technique. For both models that incorporate repression by proteolysis and repression by tracription control, we sample kinetic rate constants from biologically feasible range and uniform range. For the biologically feasible range, we sample constitutive synthesis rates from the interval [10 −7 , 10 −5 ), and enhanced synthesis and degradation rates from the interval [10 −6 , 10 −4 ). Further, repression rates were sampled from three separate intervals. For repression by proteolysis, [   We need to know the roots of the deterministic system to calculate the power spectra using Eq. (10). For some networks the roots of the deterministic system can be obtained analytically. For networks for which the steady state of the deterministic system can not be obtained analytically, we used the bisection method. We left out few networks for which it becomes impossilbe to get steady state values either by analytical expressions or by using bisection method, e.g., when two or more equations are coupled. Since the system size expansion can only be used for the system that has a single stable steady state, for proteolysis, we check the coefficients of the polynomial equations to makes sure that we have a single positive root for a given network for a given set of randomly sampled kinetic rate constants. We resample when we encounter either a negative or a very large root. This is because negative and very large steady states are biologically unfeasible. Thus by resampling we exclude such parts of domian. Maximum limit was set to a very high value (10 15 ). After obtaining a single positive state, we check for the stability of the dynamical system numerically as described in Sec. IV I. If the system is unstable, we resample k i 's from the parameter space. When the root calculation fails in the bisection method, we set SN R values to zero and do not resample.

F. SNR calculation
The SNR (signal to noise ratio) is commonly defined by the power spectrum peak height over its relative width. For species i it is thus where P i (ω o ) and ω o denote the peak height and the peak frequency, respectively, for species i. The ∆ω defines the so-called full width at half maximum (FWHM). Our calculation of the power spectrum is based on the analytical expressions discussed in section IV D. It is, however, not straightforward to obtain the analytic expression for SNR because SNR calculation requires the solution of fourth order polynomial equation given by dP i (z)/dz = 0 where z = ω 2 and the determination of ∆ω around ω o is also non trivial. Therefore, we use a numerical algorithm to calculate the SNR in the following way. First, we numerically identify the occurrence of the maxima, if it exist, of P i (ω) in Eq. (10) for a given set of kinetic rate constants and fixed points of the system. We locate the peak frequency at ω = ω o for which P i (ω o ) is the maximum value. Next, if the peak frequency ω o exists, we search for the two neighboring frequencies, i.e., ω 1 (< ω o ) and ω 2 (> ω o ) that satisfy P i (ω 1 ) = P i (ω 2 ) = P i (ω o )/2. Then, we can determine ∆ω from the difference of two frequencies, i.e., ∆ω = |ω 2 − ω 1 |.

G. Robustness vs Prominence
In this work, the robustness of a networks is defined as the size of the domain of the parameter space that generates stochastic amplified oscillations, i.e., SNR≥ 1, whereas, the prominence is defined as the average value of SNR within such domain for that network. It is to be noted that we only consider {k i } domain in which we do not have multiple roots, negative roots, very large roots and unstable roots. As discussed in Sec. IV E, when the bisection method fails, we set SNR values to zero and do not resample. This gives us lower bounds on robustness. We measure the size of the domain of the parameter space that generates stochastic amplified oscillations by counting the data sets that give SNR ≥ 1 out of total sampling points we accept. For species i, we calculate the average value of SNR within the domain by integrating the probability distribution of SNR, i.e., The step function Θ(x) is inserted to restrict the integral range only for SNR≥ 1. In numerical simulations, ( Here j is the data point and D i is the counts of the data sets with SN R i ≥ 1 for species i in the network.

H. Eigenvalue calculation
The dimensions of Jacobian matrices for all the networks considered in this parper are 3-by-3 and the corresponding third order characteristic polynomials can have only real coefficients. The three eigenvalues of the Jacobian matrix can be either all real or a combination of a real eigenvalue and a complex conjugate pair. Since all the networks we consider exhibit stable deterministic dynamics in the parameter range of our interest, the real parts of all three eigenvalues are necessarily negative. We used cubic polynomial formulae for obtaining the eigenvalues of 3-by-3 Jacobian matrices [32]. We reject any samples for which numerical accuracy is a problem, e.g., for some Jacobian matrices the numerical variables become so small that the calculation using double precision is not adequate and may need higher precision. Apart from looking at the imaginary part of eigenvaules, we can use discriminant of the matrix to determine if the eigenvalues are complex or real. Discriminant < 0 means that we have a complex conjugate pair, otherwise, all three eigevalues are real.

I. Linear stability
The proof of the stability for the deterministic system amounts to the examination of Re(λ) < 0, where λ denotes the eigenvalues of the Jacobian matrix, which are the roots of the cubic characteristic polynomials, i.e., det(J − Iλ) = 0: Here T r(J) representing the trace of the Jacobian, while A ij is the cofactor of the i th row and j th column. For all the eigenvalues to have Re(λ) < 0, the Routh-Hurwitz conditions for a third order characteristic equation, T r(J) < 0, det(J) < 0, det(J) − (A 11 + A 22 + A 33 )T r(J) > 0 should be satisfied. We used Routh-Hurwitz condition to numerically check the stability condition.

J. K-medoids clustering
We used k-medoids algorithm to cluster the data points. K-medoids is a partitioning technique of clustering a data set into a given number of clusters. We chose this method over K-means clustering because it is more robust to noise and outliers. We used the squared Euclidean distance as our distance measuer. A medoid of a cluster is the most centrally located data point of the cluster. The most common realisation of K-medoid clustering is the Partitioning Around Medoids algorithm which is also the method we used in this paper. Since the attainment of the minima is not guaranteed, we made several trials with different initial conditions and chose the clustering that gave least squared Euclidean distance.

K. Calculation of Hamming distance for the identification of common network architecture
We use the Hamming distance measure to find the similarity between any two given networks. To find the Hamming distance between two given networks, in both the networks we assign the value of +1 to an edge when it is activation, −1 to an edge when it is repression and 0 when the edge is absent. We then overlap the two networks and calculate the Hamming distance as the sum of the absolute values of the differences of the values of the edges that point in the same direction. While calculating the Hamming distance between two given networks, we also consider the mirror images of the networks, and the rotations of the networks and their mirror images and choose the Hamming distance which is least after considering these transformations.
To calculate the common network architecture of each cluster, one-by-one, we consider each network in the cluester as a reference network of the cluster. We then find the sums of the Hamming distances of the reference networks from the rest of the networks in its cluster. We choose the minimum of these sums of Hamming distances and the configuration corresponding to this minimal value. Note that the configuration may contin the rotations and mirror images of the networks because we also consider these transformations while calculating the Hamming distance. The common architecture is the average architecture over all the networks corresponding to this configuration. In the common architecture, a positive value of an edge means that the edge is activation while a negative value means the edge is repression. The magnitude of the average edge values indicate how common the edge is within the cluester. Magnitudes close to 1 means the edge is very common inside the cluester while a value close to 0 means the edge is uncommon.