On the functional diversity of dynamical behaviour in genetic and metabolic feedback systems

Background Feedback regulation plays crucial roles in the robust control and maintenance of many cellular systems. Negative feedbacks are found to underline both stable and unstable, often oscillatory, behaviours. We explore the dynamical characteristics of systems with single as well as coupled negative feedback loops using a combined approach of analytical and numerical techniques. Particularly, we emphasise how the loop's characterising factors (strength and cooperativity levels) affect system dynamics and how individual loops interact in the coupled-loop systems. Results We develop an analytical bifurcation analysis based on the stability and the Routh- Hurwitz theorem for a common negative feedback system and a variety of its variants. We demonstrate that different combinations of the feedback strengths of individual loops give rise to different dynamical behaviours. Moreover, incorporating more negative feedback loops always tend to enhance system stability. We show that two mechanisms, in addition to the lengthening of pathway, can lower the Hill coefficient to a biologically plausible level required for sustained oscillations. These include loops coupling and end-product utilisation. We find that the degradation rates solely affect the threshold Hill coefficient for sustained oscillation, while the synthesis rates have more significant roles in determining the threshold feedback strength. Unbalancing the degradation rates between the system species is found as a way to improve stability. Conclusion The analytical methods and insights presented in this study demonstrate that reallocation of the feedback loop may or may not make the system more stable; the specific effect is determined by the degradation rates of the newly inhibited molecular species. As the loop moves closer to the end of the pathway, the minimum Hill coefficient for oscillation is reduced. Furthermore, under general (unequal) values of the degradation rates, system extension becomes more stable only when the added species degrades slower than it is being produced; otherwise the system is more prone to oscillation. The coupling of loops significantly increases the richness of dynamical bifurcation characteristics. The likelihood of having oscillatory behaviour is directly determined by the loops' strength: stronger loops always result in smaller oscillatory regions.


Background
Feedback control mechanisms are vital to the robust functioning of gene regulatory and metabolic pathways. They have been extensively researched over the last two decades: we now know more about the topology and functionality of positive and negative feedback in intra-and inter-cellular systems than ever before [1]. For example, positive feedback is essential for the existence of multiple steady states within multi-scale gene regulatory systems [2] and to help prolong and amplify the response to weak signals in intracellular signalling [3]. Operating synchronously, negative feedback is found to help (1) stabilize and maintain the concentration of gene products, (2) maintain the homeostasis of gene expression rates [4], (3) improve the robustness of developing cells [3], and (4) facilitate the sustaining of oscillations of gene transcription rate [2]. Many examples of negative feedback systems exist: (1) regulatory pathways: the main repressor of the SOS regulon in bacteria Escherichia coli, LexA, represses its own production [5]; the Hes1 oscillator represses its own transcription [6]; the p53-Mdm2 network with p53 activates the Mdm2 gene and Mdm2 sequesters p53 [7]; and the tryptophan operon system with multiple negative feedback regulations [8]; (2) metabolic pathways: in the linear mandelate-acetate pathway in Pseudomonas fluorescens the acetate represses seven preceding reactions [9,10]; and, (3) signaling pathways: the NF-kB signalling pathway in which nuclear NF-kB activates production of IkBα which in turn inhibits nuclear import of NF-kB by sequestering it in the cytoplasm [11]; and circadian clocks [12]. These feedback loops orchestrate the molecular fluxes in multi-scale manner so that the organisms survive and thrive in many different environments.
Early work on negative feedback in gene regulatory systems goes back to that of Goodwin [13] who proposed an auto-repressive transcriptional model with inhibition imposed by the gene's own protein product. This initial model has provided a useful framework for later studies of systems involving negative feedback regulation [14][15][16][17][18][19]. In addition, a number of Goodwin-based models of biological oscillators characterised by one or more negative feedback loops have recently been developed, for example of the circadian clocks [20,21]. Nevertheless, our understanding of the dynamic nature of negative feedback regulation is still limited in many aspects. First, the mere presence of a negative feedback loop within a system is insufficient to understand its dynamical behaviour [22,23]. In fact, negative feedback is found to promote both system stability and oscillatory instability; the same feedback, if it is "loose", it may support stability, and if tighter, it may give rise to sustained oscillations [24]. Therefore, to gain deeper insights into the dynamical behaviours of biologically regulated systems, it is necessary for us to understand and characterize the differences among the feedback loops that influence system dynamics by systematically studying the different types of negative feedback loops that occur in these systems.
Two important factors that can be used to characterise a negative feedback loop are the feedback strength and level of binding cooperativity (nonlinearity) between an inhibitor and its regulated molecule [25,26]. In this paper we investigate the effects of these factors on system dynamics. Earlier work has only looked at the effects of changes in the cooperativity levels, but not those of feedback strength [14][15][16][17][18][19]. To coordinate complex and rich interactions within the cell, cellular systems often consist of not just one negative feedback loop, but multiple ones, entangled together. What are the functional advantages of the coupled feedback loops which evolved within the host systems? Although attempts have been made -e.g. the interplay between positive and negative feedback regulation have been shown to provide robustness and reliability to system performance [3,27] -the studies on how the coupled loops affect the molecular dynamics have been very limited. We therefore aim to explore in this paper the dynamical aspects of systems with multiple negative feedback loops in comparison with their single-loop counterparts to understand the possible functional advantages the extra loops may provide.
We consider a commonly encountered motif of the negative feedback systems in which negative feedback is imposed by the last species of the system pathway on the upstream species. We develop mathematical models to analyze for stability and bifurcation, in order to study the behaviour of these systems, confining ourselves to the analytical solutions which allow us to obtain bifurcation points dependent on the feedback strengths and nonlinearity as the parameters. Based on these results, the conditions on the feedback strength or nonlinearity for: no stability; no oscillation; stability enhancement; oscillation enhancement; and guaranteed stability (oscillation) can be established. Our analyses will lead to regimes in the parameter space in which different dynamical behaviours can be identified. In contrast to numerical methods, analytical methods facilitate an analysis of the parameter sensitivities of the system dynamics.
For clarification, here we define the meanings of robustness and stability used in this paper. Robustness is referred to as the ability of a system to maintain its functionality against internal perturbations and environmental variations. Stability, on the other hand, is only concerned with the ability to maintain the system state. Although both are important properties of living systems, robustness is a broader concept than stability, with the emphasis on system functionality rather than system state [28]. A system can preserve its function amid perturbations by actively switching between different (stable and unstable) states [28]. In this study, we focused on the stability aspect of systems.
The remaining structure of this paper is given as follows: the Methods section discusses Hill function and its use for negative feedback modelling, followed by a description of the regulatory motifs studied in the current paper and their mathematical models. The analytical methods developed for stability and bifurcation analyses are outlined in the last subsection. Detailed analyses along with discussion of pertinent results are presented in the section Results and Discussion. In the section Biological Examples, we analyze the Hes1 oscillator in vertebrates and the Tryptophan operon system in bacteria Escherichia coli as examples. Finally, we end the paper with the Summary and Conclusion section.

Modelling negative feedbacks -Hill function
In modelling biochemical systems, the rate of a reaction representing concentration change per unit time can be written as a function of the concentrations of reactants and products. There exist a number of rate laws corresponding to different types of reaction mechanisms: the mass action rate law, the Michaelis-Menten kinetics, and the Hill functions [1,29]. The level of inhibition caused by negative feedback loop due to product X can be described by the Hill function of the following form (another form of the Hill function can also be used to model activation - [1]): where the parameter K i represents the half-saturation constant (i.e. the concentration of X that gives 0.5 ratio repression). It is also commonly referred to as the dissociation constant or binding constant. The parameter n (Hill coefficient) is related to the cooperativity level of the chemical process.
The strength of feedback is inversely proportional to K i : increasing K i lowers the repression level while decreasing K i increases the repression for a given n; therefore we may define FS = 1/K as feedback strength ( Figure 1); and refer to K as the inverse feedback strength indicator. The Hill coefficient n can be interpreted as the sensitivity of the feedback loop. As n becomes larger, the Hill curve becomes more sensitive to change in X in the vicinity of X 0.5 and acts like an on-off switch. When n is very large (infinity), the Hill function resembles the step function ( Figure 1).
In gene regulatory networks, change in FS could be brought about in a number of ways: (1) by mutations that alter the DNA sequence of the binding site of X in the inhibited molecular species' promoter -even alternation of a single DNA base can strengthen or weaken the chemical bonds between X and the DNA -which will subsequently change FS; (2) by change of binding site position within the DNA. The Hill coefficient can be changed, for example, by mutations that alter number of binding sites within the DNA. It has been experimentally shown for bacteria that they can accurately tune these parameters within only several hundred generations for optimal performance when faced with environmental change [1].
Hill function for modelling inhibition Figure 1 Hill function for modelling inhibition. (a) Hill function with increasing n (FS is fixed at 10) (b) Hill function with decreasing FS (n is fixed at 2).

Regulatory Motifs and Models
Our model motifs consist of generic pathways of activation steps (reactions) with arbitrary lengths with single or coupled negative feedback loops imposed by the endproduct of the pathway. Figure 2 shows the schematic diagrams of a few example motifs. We denote a general l-species system with length l by L l ; while a system having single feedback loop on the step k is denoted by (1 ≤ k ≤ l). The system with the coupled loops on the steps k 1 , k 2 ,.., k j is denoted by .
A model system can be described by using a set of l differential equations as follows, where x j , k j , k dj (j = 1,.., l) represent the concentration of species X j , its synthesis rate, and its degradation rate, respectively. And, if feedback loop is present while H(K j , n j , x l ) = 1 if no feedback loop is present for the j th step.
The first species of the pathways X 0 is often assumed to be static, i.e. its concentration is unchanged. In most cases, it represents the gene which is the source of the pathway and activates the whole sequence of reactions. However, we also consider the cases when X 0 is variable by studying the variants of system L l with shorter lengths. The degradation process of model species is assumed to follow first-order kinetics. Degradation parameter k dj is actually an aggregated rate combining transport (or modification) and decay rate of corresponding species. The Goodwin oscilla- tor is a special case and represented by system . We do not consider time delay in these models; readers who are interested in this aspect are referred to the work of Mac-Donald [30]. We first analyse the single-loop systems with three and four variables. The coupled-loops systems are examined next. We then study the generalised, extended systems with arbitrary pathway lengths.

Analysis Methodology
Biological systems display many types of dynamic behaviours including stable steady state, sustained oscillations, and irregular fluctuating dynamics (chaos). Change of system parameters may lead to change of system dynamics. Bifurcation analysis allows one to subdivide the parameter space into qualitatively different regions within each, the system dynamics are homogeneous. Furthermore, the changes in the size and location of resulting regions due to parameters variation can be investigated.

A Summary of Stability Analysis and the Routh-Hurwitz Theorem
The stability analysis of a system consisting of a set of differential equations can be conducted by considering its dynamical behaviour in the neighbourhood of its equilibrium (i.e. steady) state. A steady state is classified locally stable if the system returns to this steady state after a sufficiently small but arbitrary perturbation. Local stability of a steady state can be analysed by linearising the differential equations around the steady state and assessing the eigenvalues of the resulting Jacobian matrix (J) [31]. For a system of differential equations If the real parts of all J's eigenvalues are negative, the steady state is said to be stable, while if any of the real parts are positive, the steady state is unstable (in this case the system oscillates if the imaginary part is nonzero).
Because J's eigenvalues are actually the roots of the following characteristic equation α i s are the coefficients -to assess the signs of J' eigenvalues, we make use of the Routh-Hurwitz theorem [31,32] which states that eigenvalues λ all have negative real parts if where

Bifurcation -a Geometrically Motivated Approach
Our aim is to establish analytical bifurcation points for the feedback strengths, Hill coefficients, and other model parameters of single-loop as well as multiple-loop systems. System stability conditions are first formulated using the Routh-Hurwitz stability criteria outlined above. These conditions are then examined using a geometrically-motivated approach. We demonstrate the method below using a system with length 3. Longer pathway systems are similarly analysed. Consider the case where all three loops are present (Figure 3g), and the equations for this system are: Denote the equilibrium values of the state variables x i , (i = 1, 2, 3) . Steady-state (equilibrium) values of the system variables can be determined by setting the right hand sides of (1) to zeros. This subsequently gives (see Additional file 1) , we have f n x f n x n ,..., , Schematic diagram of feedback motifs analysed in the paper together with their model equations Figure 3 Schematic diagram of feedback motifs analysed in the paper together with their model equations.
where In this case, the characteristic polynomial is cubic Following the Routh-Hurwitz theorem, the system is stable if and only if the following Stability Condition (SC) is satisfied: It is convenient to introduce the following composite variables Variables M 1 , M 2 , M 3 interestingly have the characteristics of activation functions. Working with M 1 , M 2 , M 3 is more straightforward than with x 1 , x 2 and x 3 directly as it spares one from having to deal with the exponential and rational forms in (2). We also have 0 <M 1 , M 2 , M 3 < 1. The equilibrium condition (2) is now simplified to Equations (4) also allow the characteristic coefficients α 1 , α 2 and α 3 to be expressed in terms of only M 1 , M 2 , M 3 and other model parameters, i.e. the synthetic and degradation rates (see equations (9) below for example). Particularly, the conditions (3) and (5) for simpler system motifs with less feedback loops can be easily derived. This means if there exist bounds M 2l , M 2h (based on equation (3)) such that M 2l ≤ M 2 ≤ M 2h , then the stability condition (3) would be equivalent to Condition (8) represents an analytical relationship between the feedback strength indicators K 1 and K 2 of the loops in action, and f 1 (M 2l , K 2 ) and f 1 (M 2h , K 2 ) are the bifurcation points of K 1 .
To determine the bounds M 2l and M 2h , note that (3) can be manipulated to take the form where g is a function whose explicit form depends on the particular system motif. For system we have Because α 1 , α 2 , α 3 > 0, (3) is equivalent to α 1 α 2 -α 3 > 0.
Substituting (9) into this relation we obtain where f(M 2 ) = M 1 and , with coefficients a 0 , a 1 and b 1 are expressions of the system parame- α α α α α f M K K f M K

Single-loop Systems
Because two-species systems are incapable of demonstrating oscillatory dynamics, we only consider systems with three species or more. Here, we present the results for the systems with a single feedback loop. To this end, the three-species systems are first considered. We then examine the four-species systems and investigate potential effects on system functional dynamics as a result of lengthening the pathways.

Three-species Systems
Three negative feedback motifs are possible for the threespecies system where the feedback loop is imposed on the first, second, and the last step of the pathway. These are denoted , and and schematically demonstrated in Figure 3a, b, and 3c, respectively. We found that systems and are both incapable of having oscillatory dynamics, regardless of their parameter values. System , essentially the Goodwin system of length three, possesses both stable and oscillatory dynamics. Switching between these dynamical regimes occurs through a Hopf bifurcation. We present below the analytical condition governing this bifurcation.
The system is stable if and only if the following condition is satisfied (see section 1.1.1 in Additional file 1 for the derivation): where K 1 versus n 1 Manipulating the inequality (11) yields an equivalent condition between the inverse feedback strength indicator K 1 versus the Hill coefficient n 1 and the remaining model parameters (i.e. the synthetic and degradation rates -see section 1.1.4 in Additional file 1 for the derivation), given below where This shows, given other model parameters' values, the existence of a threshold feedback strength (1/K 1thresh ) at which the system loses stability to an oscillatory regime. Based on (13), two-parameter bifurcation diagrams of K 1 against other model parameters can be set up. Figure 5a illustrates on the K 1 vs. n 1 plane, regions of stable and oscillatory dynamics, separated by the K 1thresh curve.
To obtain oscillatory behaviour, the feedback loop must be sufficiently strong (K 1 below the curve). Furthermore, K 1thresh saturates at high n 1 which indicates that there exists a critical value of K 1 or feedback strength above which, the system is guaranteed to have sustained oscillations regardless of the value of the Hill coefficient n 1 (Figure 5a). Denote this critical value K 1crit , we found that K 1crit = 1/A, which is the product of the system synthesis rates divided by the product of the system degradation rates (see section 1.1.5 in Additional file 1 for the derivation). Figure 5a further shows that higher cooperativity level improves the likelihood to observe oscillation, since oscillation is obtained over a wider range of K 1 , i.e. lower n 1 provides more stability. However, this improvement diminishes at high cooperativity level due to the saturation behaviour of K 1thresh .

Parameter's "Ranges of Guaranteed Stability"
Here, we define the Ranges of Guaranteed Stability (RGS) of a model parameter p with respect to model parameter q as all possible values of p that always give a stable system dynamics, subjected to arbitrary variation in q. For instance, as shown above, K 1 > 1/A or (1/A, +∞) is the RGS of K 1 with respect to n 1 .
Because M 1 < 1, equation (11) means that the system is always stable if n 1 ≤ B. This threshold value depends on the degradation rates only. It also yields the RGS of n 1 to equal (0, B], with respect to all model parameters except the degradation rates. Since B ≥ 8 for any arbitrary values of k d1 , k d2 and k d3 , the interval (0,8] becomes the "global" RGS of n 1 , i.e. with respect to all model parameters. Furthermore, for any n 1 ≥ 8, the system can be made oscillating with a proper set of the degradation rates and having a sufficiently weak feedback loop. The number of these sets is found indefinite and shown in section 1.1.3 in Additional file 1.

Effects of turnover parameters
Here, we investigate effect of the synthesis and degradation parameters on the system's bifurcation characteristics. Since A and B are symmetrical expressions, K 1thresh is also symmetrical with respect to the degradation as well as the synthesis rates. This means that all system species equally affect the system's bifurcation characteristics in K n B n An B n K thresh (a) Bifurcation diagram of K 1 against n 1 Figure 5 (a) Bifurcation diagram of K 1 against n 1 . The stable and oscillatory regions are separated by the K 1thesh curve which increasingly approaches the K 1crit line (dashed). B, indicated on the n 1 -axis, is the minimum value of the Hill coefficient at which oscillations are possible (parameter values k 1 = 1, k 2 = 2, k 3 = 1, k d1 = 1.2, k d2 = 3, k d3 = 2 were used for graphing). (b, c) Comparison of the K 1 vs. n 1 bifurcation diagrams for different scenarios. (b) K 1 vs. n 1 bifurcation diagram for the base parameter set, k 1 = 1, k 2 = 2, k 3 = 1, k d1 = 1.2, k d2 = 3, k d3 = 2, (solid); when a synthesis rate k 1 is doubled (dashed); and 10-times increased (dot). Note that the exact same effects are also obtained for changing other synthesis rates. (c) K 1 vs. n 1 bifurcation diagram for the base parameter set (solid) compared to the set in which the degradation rate k d1 is increased 10 times (dashed), and the set when the degradation rates are set identical, k d1 = k d2 = k d3 = 0.5 (dot).
spite of the fact that the feedback loop is only acting on the first reaction of the pathway.
Regarding the synthesis rates, K 1thresh changes proportionally with these parameters. Increase in the production of any of the model species therefore gives rise to a more oscillatory-prone system, indicated by a larger oscillatory region in the two-parameter K 1 vs. n 1 plane (Figure 5b). More interestingly, raising the production rate of any species results in an exactly same K 1thresh curve as raising the production rate of any other species by the same proportion. Bifurcation patterns are therefore conserved under these different changes. This knowledge is potentially useful in many cases. For example, it can facilitate the engineering of synthetic circuits with desirable dynamical behaviour; as one could effectively choose appropriate points to perturb to attain desired dynamical behaviours. It can also help in the process of parameter estimation and optimisation of synthetic circuits.
The degradation rates, on the other hand, have opposite effect on system dynamics. Higher degradation rates tend to reduce K 1thresh , leading to smaller oscillatory region ( Figure 5c) and consequently a more stable system. System stability, therefore, is most likely when model species are rapidly degraded. This is because at large k dj , rapid degradations of the state species significantly weaken the strength of the negative feedback loop that is required for oscillations. In contrast, very slow degradation makes it almost impossible for the system to obtain stability, unless the feedback loop is greatly relaxed with significantly weak inhibition strength (i.e. very high K 1 ).
In examining parameter effects on the threshold value of the Hill coefficient (B), our analysis reveals that comparable degradation rates across model species (k d1 ≈ k d2 ≈ k d3 ) leads to minimum B and thus minimum RGS for n 1 ; whereas if one is many folds greater than another (k di >> k dj , i, j ∈ {1, 2, 3}), B will be high, resulting in a large RGS (see section 1.1.6 in Additional file 1 for the justification). This suggests a way to enhance system stability by unbalancing the degradation rates of molecular species, preferably, towards higher values. Figure 5c compares the bifurcation profiles between a reference parameter set and one with equally low degradation rates; and one with unequal, enhanced degradation rates.

Four-species Systems
These systems can be considered extension of the previous models, which consist of four species. There are a total of four feedback motifs for single-loop systems. Similar to and , systems and are found to be incapable of producing oscillations. Here, we consider and in turn (Figure 3h, i).
Interestingly for , our analysis arrives at the same bifurcation points for n 1 and K 1 as in (11) and (13), however with different expressions of A and B: We found that B ≥ 4 for arbitrary values of the degradation rates. Compared with the three-species system, n 1 's RGS is reduced, supposedly due to the lengthening of system pathway. Furthermore, for n 1 > 4 the system is capable of displaying sustained oscillation for an indefinite number of parameter sets, given a proper selection (see section 1.2 in Additional file 1 for the justification).
Let us now consider system , a variant of the singleloop Goodwin system where inhibition is imposed by the end-product on the second rather than the first step of the pathway (Figure 3i). At first glance, this system design looks like its counterpart . However, there is a fundamental difference in the synthesis of the repressed variables between the two systems. The synthetic term for x 2 in contains x 1 which changes dynamically, having its rate of change defined by the first model equation in Figure 3i; whereas the synthetic term for x 1 in does not contain a varying variable; can thus be considered as a special case of by setting k d1 = 0 in Figure 3i. This difference might give rise to distinct dynamical behaviours between the two systems. Here, we analyse 's bifurcation to identify and investigate possible behavioural discrepancies compared to .
For system , the analytical bifurcation point for feed- where A has the form of (15) while B is reduced to resemble that of the three-species system in equation (12),

Generalised Single-loop Systems Minimum Hill Coefficient for Sustained Oscillation
In fact, for the single-loop Goodwin system with arbitrary length l, the minimum Hill coefficient required for sustained oscillations has been theoretically computed to be although this calculation is done under the stringent assumption of equal degradation rates k d1 = k d2 =...= k dl [15,19]. Figure 6 plots this minimum Hill coefficient against the pathway length l (≥ 3). We observe dramatic reduction of the minimum n 1 at small length (≤ 10) but this reduction becomes insignificant for longer pathway; a saturation trend is instead observed. Our derivation, however, gives us explicit form of the minimum n 1 as analytical expression of the degradation rates.

Effects of System Extension
Comparative study of the systems and allows us to examine dynamical effects resulting from system extension (i.e. the lengthening of system pathway). We found that the extended system with more species is not always more stable. In fact, whether the extended system is more stable or more prone to oscillation is determined by the kinetics (i.e. synthetic and degradation rates) of the added species.
If we denote K 1crit (L) the critical K 1 value of system L, (14) and (15) then give Minimum Hill coefficient (n 1 ) required for sustained oscillations in the Goodwin system with length l Figure 6 Minimum Hill coefficient (n 1 ) required for sustained oscillations in the Goodwin system with length l.
Equation (18) indicates that if the additional species is more quickly degraded than produced (k 4 <k d4 ), K 1crit will be reduced for the extended system. On the other hand, Also, note that if k d4 is small (large) relative to other degradation rates, n 1crit for the extended system becomes greater (smaller) (see (12) and (16)). Figure 7 demonstrates the effects of system extension on the bifurcation characteristics under different scenarios of the added species' kinetics. We conclude that the extended system obtained as a consequence of pathway lengthening becomes more stable only when the added species degrades slower than it is being produced. In this case, the feedback loop must increase its strength to a proportional level, if sustained oscillation is to be obtained (see (18)).
On the other hand, the extended systems are more prone to sustained oscillations if the additional species degrades faster and being created. These observations provide us with useful indications of how regulatory system might tune its feedback strength to achieve certain types of dynamics.

Generalisation of Feedback Strength and Hill coefficient
The equation defining the threshold feedback strength for the systems , and can be extended by induction  (13) with A's generalised form below, and B is a function of only the degradation rates k d1 , k d2 ,..., k dl . B was found for the systems with length 3 and 4 to be neat expressions of the degradation rates. However, for the system with 5 variables or more, B becomes complex but can be derived explicitly. Moreover, B reaches its minimum when all degradation rates are equal, at which B = 1/Cos l (π/l). In addition, the generalised critical value of the feedback strength is A. A system with feedback strength weaker than this value is guaranteed a stable dynamics independent of the Hill coefficient values.
Equation (17) shows that the threshold Hill coefficient of only depends on the degradation rates of the inhibited species and its downstream molecules; and is independent of the upstream species. More interestingly, equations (12) and (17) show that and share the same analytical form for B. Further generalisations, therefore, can be made for variants of the general Goodwin systems: the systems with single negative feedback loop in which any step in the biochemical pathway can be potentially inhibited by the end-product (i.e. systems with 1 <= k <= I - Figure 8a). The threshold K kthresh of is where A is retained its form as in (19) whereas B involves only the degradation rates of downstream species of the repressive targeted species X k .
More importantly, B has similar form as that of the reduced Goodwin system ( Figure S1 in Additional file 1). We confirmed these generalised analytical equations using numerical computations as well in which we estimate K kthesh for variant systems and found that they fit the theoretical form given by .

Effects of Feedback Loop Reallocation
The generalised findings above have important implications concerning the dynamical behaviour of feedback systems. Comparing (19) and (20) reveals that reallocation of the negative feedback loop has no effect on the critical feedback strength (1/K crit ). Regardless of the loop's position, the system's stability is ensured if the feedback strength is weaker than this value.
sional bifurcation diagram K 1 vs. n 1 could either shrink or expand due to loop reallocation.

Non-oscillatory Systems
Because the systems with less than three species are not able to display oscillatory dynamics, the generalised results obtained in the previous section indicate that for We will show in the next section that, adding extra inhibition loops tends to make the system less likely to demonstrate sustained oscillations. This implies that any system with multiple negative feedback loops encompassing one imposing on either the last or second last step ( or with k 1 ,.., k m ∈ {1, 2,.., l}) is also incapable of having oscillatory dynamics (Figure 8b).

Coupled-loop Systems
Let us consider systems with increased complexity in which several negative feedback loops are coupled together. Understanding of composite behaviour of coupled loops has been limited. We aim to establish meaningful connections between the feedback strength of these loops (through the inverse indicator parameters Ks) under certain conditions.

Three-species, Doubled-loop System
Detailed analysis for this system was presented in section "Analysis Methodology" as an example (Figure 3d). The stability condition (10) is ensured if g(M 2 ) intersect the M 1 -axis at a point above point (0,1), indicated by a dot in Figure 4a. This translates a 0 ≥ b 1 or n 1 ≤ B with B as in (12). Compared to the results from the previous section of the three-species, single-loop system, we found that adding a second feedback loop does not affect the RGS of n 1 . On the other hand, if n 1 > B, the line g(M 2 ) must intersect f(M 2 ) within the unit square U and so condition (10) is violated for some M 2 , subsequently destabilising the system (Figure 3a).

K 1 -K 2 Bifurcation Diagram
Range of M 2 satisfying (10) can be determined, indicated by its lower bound M 2l and higher bound M 2h in Figure 4a. For each value of K 2 , we obtain corresponding values for M 2l and M 2h . Using (8), we construct the two-parameter bifurcation diagram with the feedback strength indicators K 1 and K 2 being the axes.
Note that in this case, there exists at most one (unique) intersection point between f(M 2 ) and g(M 2 ) within U, indicating a simple binary separation of the K 1 -K 2 bifurcation diagram into stable and oscillatory regions. A typical K 1 -K 2 bifurcation profile for is illustrated in Figure 9a.
We refer to the feedback loop involving K 1 and K 2 as loop L 1 and L 2 , respectively.
As loop L 1 is relaxed (larger K 1 ), sustained oscillation becomes more difficult to obtain at stronger L 2 (indicated by raised minimum K 2 possible for oscillation). Stability is most likely under weak L 1 coupled with strong L 2 . Oscillations, on the other hand, are most likely at weak L 2 coupled with strong L 1 ; shifting system dynamics to being stable at strong L 1 , however, requires L 2 to be very strong too. Dynamical behaviours are summarised in Table 1 based on combinations of the individual loop's strength. Figure 9a also shows that there exists a critical value for K 1 above which the system must be stable. This means stability is guaranteed if L 1 is sufficiently weak, regardless of the nature of the second loop L 2 . Unlike L 1 , stable as well as oscillatory dynamics can be obtained at any strength of L 2 , given the proper specification of L 1 (Table 1). The critical K 1crit mentioned above was found to have the exact same expression as the K t1hresh in (13) Note that K 1crit is independent of L 2 's specification, whereas it increases with n 1 . As a result, higher n 1 enhances oscillatory behaviour due to the expansion of the oscillatory region. Increase n 1 therefore enables oscillatory exhibition at weaker L 1 (Figure 9b). Moreover, we found that increasing n 2 does also expand the oscillatory region (Figure 9c), enabling oscillatory exhibition at stronger L 2 . Therefore, for a coupled-loop system, raising the Hill coefficient of any loop tends to enhance oscillatory behaviour.

Effects of turnover parameters
The above equation indicates that K 1crit increases proportionally with the synthesis rates. This causes the oscillation region to approximately increase by the corresponding proportion. Oscillatory dynamics is now achievable at higher K 1 given fixed K 2 (Figure 10a). On the other hand, comparable degradation rates (k d1 ≈ k d2 ≈ k d3 ) leads to low B (see section 1.1.6 in Additional file 1) and as a result raises K 1crit . Particularly, K 1crit is maximised when this comparable rate is minimised. Whereas, if these parameters are different by many folds, K 1crit is small and Effects of the synthetic and degradation parameters on the two-parameter K 1 vs. K 2 bifurcation diagram Figure 10 Effects of the synthetic and degradation parameters on the two-parameter K 1 vs. K 2 bifurcation diagram. (a) Comparison of the K 1 vs. K 2 bifurcation diagram between the reference parameter set (k 1 = 1, k 2 = 2, k 3 = 1, k d1 = 1.2, k d2 = 3, k d3 = 2) and when a synthesis rate (k 1 ) is doubled; (b) Comparison of the K 1 vs. K 2 bifurcation diagram for 3 parameter sets: the reference set above (solid); when the degradation rates are all equal = 1.2 (dashed); and when one degradation rate is much greater than another (dot).
so is the corresponding oscillatory region. Figure 10b compares three scenarios in this case. The parameter set with k d1 = k d2 = k d3 = 1.2 obtains the largest oscillatory region while setting k d1 = 4k d2 has it significantly diminished.

Three-species, Doubled-loop System
An alternative doubled-loop system is considered here in which the first and the last pathway steps are inhibited (see Figure 3e for the schematic diagram and model equations). Follow the analytical methodology in section 2; we obtain the following stability condition for the system with and , where coefficients a 0 , a 1 , a 2 and b 1 are given in the supplementary material (section 4.3, Additional file 1).

K 1 -K 3 Bifurcation Diagram
Similar analysis was carried out on the two-dimensional

Three-species, Multiple-loop System
In this section, we consider the system structure which incorporates all three feedback loops imposing on all pathway steps (see Figure 3g for the circuit diagram and model equations). The analysis becomes more complicated, as a result. For this multiple-loop system, the stability condition is given by Due to space restriction, we give the explicit forms of function f and g in the supplementary material (section 4.4, Additional file 1). K 1 was derived as a function of the remaining model parameters: Following the similar methodology laid out above, we were able to compute bifurcation diagrams for any pair of feedback strengths (K 1 vs. K 2 , K 1 vs. K 3 , and K 2 vs. K 3 ). In all the cases, it is found that having extra third loop always increases the stability of the system, illustrated by expansion of the stability region on bifurcation diagrams. The likelihood of obtaining oscillatory dynamics is directly controlled by strengths of the loops in effect, with stronger feedback loops always result in smaller oscillatory region. Figure 11b compares the bifurcation profiles of the doubled-loop and three-loop systems.

Four-species, Coupled-loop System
To further understand the effects of loops coupling in negative feedback systems, we examine a system consisting of four species with two loops in action, denoted (see Figure 3l for the circuit diagram and model equations).
Our derivation arrives at the stability condition for this system, given below In this case f(M 2 ) still have the usual form with A as in equation (15). However, we have where the coefficients a 0 , a 1 and a 2 are given in the supplementary material (section 4.5, Additional file 1).
In a similar fashion, to construct bifurcation diagram of the loops' feedback strengths, we need to determine the ranges of M 2 that satisfies (22) via examining the curves in (23) and (24) on the M 1 -M 2 coordinate. System stability occurs over the ranges of M 2 such that f(M 2 ) lies below g(M 2 ), while oscillatory dynamics reigns over the remaining ranges.
As expected, system extension greatly complicates the analysis due to the increased number of parameters and the increased complexity of g(M 2 ).

Loops Coupling lowers Hill coefficients for Oscillations
To find the RGS for the Hill coefficients n 1 and n 2 , we determine on the two-parameter Hill coefficients n 1 -n 2 plane regions that give rise to system stability regardless of the feedback loops' strengths and other parameters. These regions can be referred to as the Regions of Guaranteed Stability for n 1 and n 2 .
From Figure 12a we can see that condition (22) will always be satisfied if g(M 2 ) always lies above the unit square U and therefore above f(M 2 ) for all K 2 . The RGS for n 1 , n 2 can thus be found by solving the following system of inequalities,  Put , which only depends on the degradation parameters. We obtain the RGS as containing those points (n 2 , n 1 ) on the n 1 -n 2 plane such that where B 1 and B 2 have the same form as in (16) and (17), respectively. Specifically, This region is presented in Figure 12b as the marked area bounded by the axes (with B 1 and B 2 ). Figure 12b shows the RGS region for three different sets of the degradation rates. Interestingly, we found that for arbitrary values of degradation parameters, the corresponding RGS always contain the 8 by 4 triangular, indicating that any (n 1 , n 2 ) inside this triangular guarantees system stability regardless of all other model parameters including the feedback strengths. On the other hand, for (n 1 , n 2 ) outside the tri-angular, we can always find a set of k di (i = 1,..,4) so that its RGS does not contain (n 1 , n 2 ), giving rise to unstable system equilibrium.
Recall that in the cases of single-loop systems considered before, there exist lower bound conditions for the Hill coefficients if oscillatory dynamics is to be obtained. For example, n 1 must > 4 for , and n 2 must > 8 for . Feedback loops coupling, however, effectively removes these constraints for the Hill coefficients. In fact, sustained oscillation is now achievable for any value of n 1 (n 2 ) given proper choice of n 2 (n 1 ). As a result, sustained oscillation can occur at much more biologically plausible values of n 1 , n 2 ; e.g. (n 1 , n 2 ) = (3, 3) or (2,4), indicated by the dots in Figure 12b.
It is important to note that the RGS for (n 1 , n 2 ) solely depends on only the degradation rates. Variation on these rates affects its size and location. B 2 is maximised if among the degradation parameters, one is many folds greater than another (k di >> k dj with i, j ∈ {2,3,4}). Similarly, B 1 is maximised if k di >> k dj . On the other hand, B 1 , B 2 are lowest when k di >> k dj with k d1 ≈ k d2 ≈ k d3 ≈ k d4 . Therefore, reducing any degradation rate to extreme low or high level will expand the RGS, resulting in system stability for wider range of Hill coefficients. Oscillation, consequently, is enhanced when the degradation rates are kept comparable between the system species.

Loops Coupling Generates 13 different K 1 -K 2 Bifurcation Patterns
If (n 1 , n 2 ) lies outside the RGS region, g(M 2 ) must cross U and therefore must intersect with f(M 2 ) at least once at some value of the feedback strength indicator K 2 . Unlike the previously considered systems where only one intersection point is detected, the number of intersection points in this case could be up to three. This provides a rich variety of different bifurcation profiles for the system. In fact, we identify a total of 13 distinct patterns of bifurcation on the K 1 vs. K 2 bifurcation plane; each pattern for one choice of the Hill coefficients. These 13 patterns are displayed in Figure 13.
These bifurcation diagrams differ in their characteristics: the shapes of the stable and oscillatory regions. For example, Figure 13a shows a simple bifurcation profile: for each and every value of K 1 , the system displays oscillations over a range of K 2 with a lower bound but no upper bound. Figure 13c shows a similar feature but the range of K 2 for oscillations is now bounded by both the lower and upper bounds; moreover, this range only exists for a certain range of K 1 . Figure 13k displays even more intricate bifurcation characteristics: as the parameter K 1 moves up the vertical axis, the corresponding set of values of K 2 for oscillations continually changes with no, one bounded range, two bounded ranges, and one unbounded range. This indicates the complexity between the feedback loops' strengths in contributing to shaping up the dynamics of the system, as a whole.

End-product Utilisation
In the preceding model systems, reduction of model species was assumed to occur only through degradative processes (decay and modification). However, species' reduction could also occur via other mechanisms. We consider here the scenario where the end-product (inhibitor), besides being degraded, is due to be consumed by the cell for synthesising other cellular components. Notable examples are of common amino acid biosynthesis pathways in which the amino acid (pathway's end-product) is utilised by the cell for protein synthesis [33][34][35]. Earlier work [16,35] suggests that the change of the endproduct in this manner has important effects on the dynamical stability of the system. However, these work were numerical; there are no analytical analyses of this effect.
Although the inclusion of g complicates our analysis, we were able to obtain the stability condition for the system in simple form similar to (11) (see section 3 in Additional file 1 for detailed derivation) with B g is given by (note that A and B are as in (15) and (16) of system ): As shown in Additional file 1, section 3, for the system to have equilibrium, the utilisation rate g must not exceed a critical value

End-product utilisation enables oscillations at any Hill coefficient
With (27), it is easy to check that B g < 1 for all n 1 and g.
Hence, for any n 1 , stability condition (3.24) can be breached by choosing K 1 sufficiently large (see Additional file 1, section 1.1.3), and so the system is destabilised (oscillatory). Moreover, this is true for arbitrary value of g > 0. The interesting implication here is that, unlike system where sustained oscillation is only attainable for certain n 1 (n 1 must be greater than 4), the inclusion of g, even small, has enabled the system to attain oscillation at any n 1 . End-product utilisation therefore allows oscillatory dynamics at low cooperativity level. This is demonstrated in Figure 14b where bifurcation diagrams on the K 1 -n 1 plane are compared for system with (thick line) and without (thin line) end-product utilisation. The bifurcation diagrams were constructed based on the threshold K 1 which we calculated to be (see Additional file 1, section 3 for the derivation): This threshold approaches a critical value as n 1 increases; K 1crit is given by (1-G)/A which is smaller than K 1crit of (1/A). This indicates in order to achieve oscillation; the feedback strength of the system with end-product utilisation must generally be stronger (Figure 14b).

Effects of g on systems dynamics
We expect that change in g would bring change in the dynamical characteristics of the system. We found that this change comes about in an interesting way. As g (G) increases, K 1crit reduces, causing shrinking of oscillatory region, especially at high n 1 (Figure 14c). Therefore, higher utilisation of the end-product generally requires .
stronger feedback loop if oscillation is to be obtained. However, when G exceeds 0.5, the oscillatory region changes its shape significantly resembling an L shape, indicated by the crossed area in Figure 14c (Figure 14c compares when G = 0.8 against G < 0.5). Now at low n 1 , oscillatory region is greatly expanded, enabling oscillation at a much wider range of K 1 .
We also present on Figure 14d the two-parameter K 1 vs. G bifurcation diagrams for n 1 = 2, 4 and 6. In this case, higher cooperativity generally gives rise to larger oscillatory region, consequently promoting oscillation. We observe an intermediate value of G (and so g) for which oscillation is most likely (widest range of K 1 ) while at low and high G the system tends to be more stable. This is consistent for all three plotted n 1 . As n 1 increases, this intermediate g moves further left in its spectrum (between 0 and g c ) To sum up, we showed that end-product utilisation enhances sustained oscillation at low cooperativity level while it enhances stability at high cooperativity. However, it is important to note that raising the utilisation level does not always further these enhancements. In fact, there exists an intermediate rate for utilisation at which sustained oscillation is detected most likely, while less likely at other rates.

Biological Examples
Restricted by the paper's scope, we only present here two biological examples from the literature where useful insights can be readily obtained by applying the analysis outlined in this paper. These are (1) the Hes1 oscillator which plays an important role during somitogenesis of vertebrates and (2) the Tryptophan operon system responsible for regulatory production of tryptophan in Escherichia coli. The former system exemplifies single-loop system while the later is an instance of multiple-loop system.

The Hes1 Oscillator
A wide range of cellular phenomena have their activities centred on oscillations [36,37]. One such notable example is vertebrate somitogenesis. This is a developmental process in which the vertebrate embryo becomes segmented by the regular sequential assignment of mesodermal cells to discrete blocks [38]. Experimental evidence reveals the basic helix-loop-helix (bHLH) transcription factor Hes1 as an important cyclic gene driving this oscillations [6,39,40]. These studies showed that the oscillatory expression of the bHLH factor Hes1 is regulated by a direct negative feedback loop whereby Hes1 represses the transcription of its own coding gene (Figure 15a) A few models have been developed for this network [6,[41][42][43][44]. Here we base our analysis on a model suggested by Zeiser et al. [43] which consists of four ordinary differential equations involving Hes1 mRNA and protein and incorporates their transportation processes between the nucleus and cytoplasm. Using notations in the original paper, the model equations are given below, Here m 1 and m 2 represents the concentration of Hes1 mRNA before and after being transported from the nucleus to the cytoplasm, respectively; while p 1 and p 2 are the concentration of Hes1 protein before and after being transported from the cytoplasm to the nucleus, respectively. Equation (29a) describes the synthesis of mRNA in the nucleus. The mRNA is then transported into the cytoplasm, described by (29b). Translation into protein is specified by (29c) while ( Table 2. It is easy to see that model (29) is just a special case of system model analysed earlier with parameters adapted and presented in Table 3. Consequently, the threshold Hill coefficient is computed based on equation (16) using only the degradation parameters gives a value of about 5.6. This means for the Hes1 oscillator to be oscillating at all, the Hill coefficient must be greater than 5.6. The Hill coefficient used by Zeiser et al., 6.2 is quite close to this minimum value (Table 3 and Figure 16b). We constructed the two-parameter K 1 -n 1 bifurcation diagram in Figure   16a. The threshold feedback strength can also be readily calculated from equation (15), as K 1thesh = 1/A = 3296.
Thus, for the Hes1 system to be an oscillator, the necessary condition for K 1 is that K 1 < 3296, given other parameters'  Table 3. Moreover based on Figure 16a, at n 1 = 6.2, viable range of feedback strength for the Hes1 oscillator is 0 <K 1 < 207.6. Zeiser et al. therefore used a quite small K 1 (i.e. 10). Interestingly, we find that variation of the feedback strength (K 1 ) has little effect on the oscillation period. Figure 16c compares the temporal change of Hes1 protein concentration, p 2 in (29), for the parameter set in Table 2 and when feedback strength is 100-fold stronger (K 1 = 0.1), and 10-fold weaker (K 1 = 100). So in fact, constraining the oscillation period to be 120 mins can gives rise to many more suitable parameter sets other than one in Table 3.

The Tryptophan Operon System
The tryptophan operon system in E. coli controls the production of tryptophan amino acid inside the cell. Key molecular processes include transcription, translation and synthesis of tryptophan. To regulate these processes, the tryptophan operon utilises three negative feedback mechanisms: transcriptional repression, attenuation, and enzyme inhibition [8,45].
The transcription process is initiated as RNA polymerase binds to the promoter. However, when the activated form of repressor which is induced by the attachment of two tryptophan molecules become abundant, it will bind to the operator site and block RNA polymerase from binding to the promoter, thereby, repressing transcription and forming the first feedback loop. Furthermore, transcription can also be attenuated depending on the level of intracellular tryptophan and is controlled by the leader region sitting between the operator and the genes ( Figure  15b). This attenuation makes up the second feedback loop. The tryptophan operon consists of five structural genes positioned consecutively after the leader region. These genes code for five polypeptides that make up enzyme molecules in the form of tetramers, which in turn catalyse the synthesis of tryptophan from chorismates [8,29,45,46].   Hill coefficient values (> 50) and at weak or strong feedback loops. This suggests that the tryptophan system is extremely stable.
The highly stable property of the tryptophan system is probably underlined by the fact that it is regulated by multiple feedback loops in concert. In addition, the system's degradation rates k d1 , k d2 and (k d3 +g) are significantly different to each other as shown in Table 4. This disparity in value of the degradation rates, as discussed earlier, greatly enhances system stability. On the other hand, by adjusting k d1 , k d2 , k d3 and g so that k d1 ≈ k d2 ≈ k d3 +g, oscillatory behaviour can now be observed at much lower Hill coefficient and at appropriate feedback strength of the loops. For example, setting k d1 = k d2 = k d3 +g = 15 can give rise to oscillatory dynamics with n 2 , n 3 as in Table 4 and n 1 as low as 8.5.

Summary and conclusion
Previous studies [14][15][16][17][18][19] have looked mainly at the effect of cooperativity level on system dynamics, while largely neglecting the effects of feedback strength. Furthermore, most of these analyses were carried out numerically; those with analytical approaches were however often done under great simplification for model equations such as stringently assuming that all degradation parameters are identical [15,19]. The major contributions from our study are summarised and discussed below.

Threshold feedback strength
For single-loop systems where inhibition is fedback by the end-product on the first reaction step, i.e. the original Goodwin system, it was found that oscillatory behaviour is only obtainable if the feedback loop is sufficiently strong. Otherwise, the system is stable and achieves steady state. Switching between these dynamics occurs through a Hopf bifurcation. We derived an explicit, analytical form for the feedback strength's bifurcation point which can be straightforwardly computed if the other model parameters are known. Interestingly, this threshold strength was found to follow a saturation trend and approaches a critical level as the Hill coefficient increases. We further showed that this critical feedback strength equals the ratio of the product of the degradation rates and the product of the synthesis rates. So for a system with feedback strength weaker than this critical level, system stability is guaranteed regardless of how high the Hill coefficient is.
Studying the two-parameter bifurcation diagram with the feedback strength and the Hill coefficient as parameters revealed that as the Hill coefficient is raised, sustained oscillation can be obtained over a wider range of feedback The Hes1 oscillator   [15] and Goldbeter [49].

Threshold Hill coefficient
Assuming identical degradation rates, previous studies of the Goodwin system have identified a threshold value for the Hill coefficient lower than which, sustained oscillation is unachievable. Under this over-simplifying assumption, the threshold cooperativity level was found only dependent of the length of the system pathway. However, such a similar threshold Hill coefficient has not been analytically determined for systems with general, arbitrary degradation rates. We explicitly derived this threshold for the three and four-species systems, which turns out to be rather simple, symmetrical functions of only the degradation rates and are not influenced by the synthesis rates. More importantly, the threshold reaches its minimum when all the degradation rates are equal. The threshold Hill coefficient obtained under identical degradation rates, therefore, provides a lower bound for that of the general case. We further showed that for the systems with Hill coefficient exceeding the threshold value, it can always oscillate with a properly chosen set of parameters' values. In fact, the number of parameter sets giving rise to sustained oscillations is indefinite. For the systems with more species, the threshold Hill coefficient is also explicitly derivable; however, the form of the resulting function becomes significantly more complex. Nevertheless, the above results also apply for longer systems which we confirmed using numerical simulations.

Effects of parameters variation
Parameter sensitivity analysis revealed interesting effects of parameters (the synthetic and degradation rates) variation on the dynamical characteristics of the system. Because of the symmetry in the expressions involving the threshold feedback strength and Hill coefficient, the individual model species equally characterise the system bifurcation profiles despite the fact that the feedback loop is only acting on the first reaction step. Specifically, increasing the synthesis rate of any model species by the same proportion results in a proportionally larger oscillatory region and hence in a system which is more likely to oscillate. In contrast to this simple linear relationship between the synthesis rates and the bifurcation profiles, the degradation rates affect system dynamics in a more intricate manner. We found that system stability is most likely when the model species are rapidly degraded while slow degradation only leads to stability if the feedback strength is significantly weak. We further showed that having comparable degradation rates between the model species promotes oscillations, whereas stability is promoted if one rate is significantly larger than another. This sug-gests a way to enhance system stability by unbalancing the degradation rates; preferably, towards high levels. These results are particularly helpful for the engineering of synthetic circuits with desirable dynamical behaviour, as well as for parameter estimation and optimisation.

Feedback reallocation
For single-loop systems, reallocation of the feedback loop to inhibit a reaction step further downstream may or may not make the system more stable. Interestingly, the specific effect is determined only by the degradation rates of the model species downstream of the newly inhibited species. The dynamical properties of the new system closely resemble those of the Goodwin system with reduced length, which equals to the number of species downstream of the inhibited species. Therefore, as the loop moves closer towards the end of the pathway, the minimum Hill coefficient for oscillation is reduced. In addition, we found that feedback reallocation does not influence the critical feedback strength discussed above. This means that for a system possessing a loop weaker than this strength, its stability is ensured regardless of the loop's position and the cooperativity level.

System extension
It has been known that lengthening the system by increasing the number of reaction steps, i.e, increasing number of model species, reduces the cooperativity necessarily required for sustained oscillations [15,19]. The implication here was that system extension enhances oscillations. However, this result was demonstrated under the assumption of identical degradation rates. When this assumption is relaxed, we found that the extended system is not always more stable. More importantly, whether it is more stable or not is attributed to the kinetics of the added species: more stable only when the added species degrades slower than it is being produced; and more oscillationprone otherwise.

Non-oscillatory systems
It has been known that the systems with two species are incapable of exhibiting sustained oscillations, regardless of the feedback strength and the Hill coefficient value [14,15]. Our analysis further showed that those systems with arbitrary lengths which possess a single loop inhibiting either the last step or second-last step of the pathways is also incapable of obtaining oscillatory dynamics. In addition, multiple-looped systems which include at least one loop inhibiting either the last or second last step of the pathway is also incapable of demonstrating oscillatory behaviours.

Effects of end-product utilisation
We also investigated the situation when the pathway's end-product is used up by the cells, which is common in many metabolic pathways. Most interestingly, we showed analytically that end-product utilisation enables oscillatory dynamics at any Hill coefficient value. More specifically, end-product utilisation enhances sustained oscillation at low cooperativity level but enhances stability at high cooperativity level. It is important to note that raising the utilisation level does not always further these enhancements. In fact, there exists an intermediate rate for utilisation at which sustained oscillation is most likely to be detected, while being less likely at other utilisation rates.

Effects of loops coupling
Since cellular systems are complex and often consist of multiple, interlocked feedback loops. Understanding of how the loops act together in giving rise to the system dynamics is absolutely crucial. Designs with interlinked positive and negative feedback loop have been shown to exhibit performance advantages over simple negative feedback loops, such as the ability to easily tune frequency of oscillators, improved robustness and reliability, even under noisy environments [22,27,50]. Multiple-negativefeedback-loop designs have also been shown to enhance system robustness and generates developmental constancy [8,27,47,51]. We obtained in this study a number of results which further our understanding into the dynamics of coupled-loop systems. We discuss these below.
First, coupled loops effectively enable oscillations at lower, more biologically plausible Hill coefficient value. For example, a four-species single loop Goodwin system requires the Hill coefficient (n 1 ) to be at least 4 for oscillations. Its variant design with the loop reallocated to impose on the second pathway step requires the Hill coefficient (n 2 ) to be at least 8 for oscillations. However, a system with both of these loops in effect can achieve oscillations at practically any Hill coefficient value for one loop, given proper choice of the Hill coefficient for the other loop. Oscillations, therefore, are possible at more biologically plausible Hill coefficients, for example at (n 1 , n 2 ) = (3, 3) or (2,4).
Reduction of the Hill coefficient for oscillations is often only suggested via pathway lengthening by previous studies. In this study, loops coupling and end-product utilisation (discussed above) were shown as the two additional mechanisms where this reduction can be obtained without increasing the number of system variables.
Secondly, coupled-loop systems were also shown to exhibit much greater complexity and more diverse behaviours compared to their single-loop counterparts. We showed that, by having just two loops performing cooperatively, the four-species system demonstrates a rich diver-sity of dynamical characteristics. For example, we detected a total of up to13 different bifurcation patterns between the feedback strengths. This enhancement in behavioral complexity and diversity might be the reason why evolution has driven some systems to acquire multiple feedback regulations as it will increase the chance of organisms' survival when facing fluctuating environments.
Thirdly, we found that different combinations of feedback strengths of individual loops give rise to different dynamical regimes. For three species with double loops acting on the first and second steps, stability is most probable when a weak first loop is coupled with a strong second loop. Oscillations, on the other hand, are most likely if a weak second loop is coupled with a strong first loop.
If oscillations are to be obtained with a strong first loop, the second loop must also be significantly strong.
Fourthly, we found a threshold strength for the first loop. If the loop is weaker than this threshold, the system is always stable regardless of the strength of the second loop. This threshold strength turns out to be independent of the second loop's specification (its strength and cooperativity level). On the contrary, at any strength of the second loop, stable as well as oscillatory dynamics are obtainable given a proper choice of the first loop's strength. By further considering the coupled-loop system consisting of the loops on the first and the third reaction step, we discovered that the location of the additional loop has no influence on the threshold strength of the first feedback loop.
Finally, examining the system with all three loops in action showed that incorporating the extra third loop always enhances system stability. The likelihood of having oscillatory behaviour is directly determined by the loops' strength: stronger loops always result in smaller oscillatory region.
We demonstrate the practicality of our analysis by including a brief investigation of two example systems: the Hes1 oscillator and the Tryptophan operon system. The former system represents a single-loop system while the latter represents one with multiple negative feedback loops coupled together. Because of the abundant number of biological systems regulated by negative feedback loops (and many can be represented under simplifying assumptions by one of the motifs considered here) the methods developed in this study may prove useful in gaining better understanding of their dynamical behaviours.

Authors' contributions
LKN devised the work, carried out the mathematical research and implemented the numerical simulations with guidance from DK. LKN and DK wrote the paper.
Both authors have read and approved the final version of the manuscript.