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–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–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.

Methods

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/Kas 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].

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 end-product 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 oscillator 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 MacDonald [30].

Despite their simplicity, these models can readily be used to model real biological systems. For instance, a generic model with three variables can be interpreted as follows: the first equation as the synthesis of nucleic mRNA, the second equation describing transportation of mRNA to the cytoplasm, while the last equation explains translation of mRNA into protein. Extended four-variable system can be interpreted by including a fourth equation describing transportation of protein back to the nucleus.

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:

(1)

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)

(2)

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:

(3)

It is convenient to introduce the following composite variables

(4)

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

(5)

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. For example, setting M_{i} = 0 for some index i gives rise to a system structure lacking the corresponding feedback loop, e.g. M_{1} = 0 gives , M_{1} = M_{2} = 0 gives . Equations (4) lead to:

(6)

Combined with (5) and (6), each of M_{1}, M_{2} and M_{3} can be expressed as functions of the others involving only K_{i}s and n_{i}s. For example, assume M_{3} = 0 for simplicity (system ), we have

Here, M_{1} is a strictly decreasing function of M_{2} over (0, 1). Substitute this into , we obtain:

(7)

as a function of K_{2} and M_{2}. Moreover, K_{1} is strictly increasing with M_{2} since the derivatives with respect to M_{2} of the terms inside the brackets in (7) are positive over (0, 1) (see section 4.1 in Additional file 1).

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

(8)

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

(9)

Because α_{1}, α_{2}, α_{3} > 0, (3) is equivalent to α_{1}α_{2} - α_{3} > 0. Substituting (9) into this relation we obtain

(10)

where f(M_{2}) = M_{1} and , with coefficients a_{0}, a_{1} and b_{1} are expressions of the system parameters and given in the supplementary material (section 4.2, Additional file 1).

Note that (10) is to take different forms depending on specific motifs of feedback loops. Next, we analyse the inequality (10) for the system using a geometrical approach in which the curves f(M_{2}) and g(M_{2}) are drawn on the two dimensional M_{1}-M_{2} coordinate plane (Figure 4a). f(M_{2}) is a strictly decreasing curve contained within the unit square U = {(0,0); (1,0); (1,1); (0,1)} (indicated with dashed boundary in Figure 4a), g(M_{2}) is a straight line with a positive slope. As 0 <M_{1}, M_{2} < 1, the analysis is constrained within U only. Range of M_{2} satisfying (10) can be determined along with its lower and upper bounds M_{2l} and M_{2h}, illustrated in Figure 4a.

Results and discussion

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 three-species 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):

(11)

where

(12)

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

(13)

where

(14)

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 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:

(15)

(16)

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 single-loop 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 feedback strength K_{2} is

where A has the form of (15) while B is reduced to resemble that of the three-species system in equation (12),

(17)

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

(18)

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, K_{1crit} for is raised if k_{4} > k_{d4} and unchanged if k_{4} = k_{d4}. 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 to a system with an arbitrary number of reaction steps. The threshold feedback strength of the Goodwin system with the general length l ( – schematic diagram in Figure 2), for example, can be expressed as in equation (13) with A's generalised form below,

(19)

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

(20)

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.

However, loop reallocation may or may not make the system more stable. This is because the threshold Hill coefficient values B in (19) and (20) are different. The ratios of B for and based on equations (16) and (17) is given by

For simplification, let k_{d2} = k_{d3} = k_{d4} = 1. We observe that B() >> B() when k_{d1} << 1; B() <B() when k_{d1} ≈ 1 but B() approaches B() when k_{d1} >> 1. This means B could become either larger or smaller depending on the relative magnitude of k_{d1} compared to k_{d2}, k_{d3}, k_{d4}. Consequently, the oscillatory region on the two-dimensional 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 the systems with arbitrary length which possess a single loop inhibiting either the last step or second-last step of the pathways ( or ), oscillatory dynamics is also not feasible. This conclusion is in line with analytical results obtained for the systems , , , previously.

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 K s) 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.

K_{1}'s Critical Value

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 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

(21)

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 M_{1}-M_{3} coordinate to determine the ranges of M_{3} satisfying condition (21). Figure 4b shows that, in this case, g(M_{3}) is a concaved-up parabola instead of a straight line like g(M_{2}) in the previously considered system . Nevertheless, there exists still at most one intersection point between f(M_{3}) and g(M_{3}), resulting in a similar bifurcation pattern for the K_{1} vs. K_{3} diagram as in (Figure 4a). Moreover, the critical value for K_{1} discussed previously is found to have the same form here. This again confirms that incorporation of additional feedback loop does not affect K_{1crit}, regardless of the location of the added loop. Moreover, the first loop's Hill coefficient (n_{1}) also has its RGS unchanged: RGS = (0, B]. The two-parameter K_{1} vs. K_{3} bifurcation diagrams were constructed based on the following relation:

Here, we compare the bifurcation profiles between two doubled-loop systems and . To facilitate this, we impose K_{2} = K_{3} and n_{2} = n_{3}. The line g(M_{2}) is superimposed on the M_{1}-M_{3} plane (by setting M_{2} ≡ M_{3}) and indicated by the thin line in Figure 4b. Note that that g(M_{2}) and g(M_{3}) meets on the M_{1}-axis. Moreover, the slope at this point for g(M_{3}) is always steeper than g(M_{2}), suggesting a higher lower bound for M_{3} for stability. The implication is: given the same set of parameter values, adding loop L_{3} results in a larger stability region (a smaller oscillatory region) than adding loop L_{2}, therefore better enhance system stability. On the other hand, is more likely to exhibit oscillatory dynamics than .

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

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,

(25)

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

(26)

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 triangular, 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–35]. Earlier work [16, 35] suggests that the change of the end-product 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

(27)

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):

(28)

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–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,

(29a)

(29b)

(29c)

(29d)

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 (29d) represents transportation of the protein into the nucleus where it represses its own transcription. Parameters b and a denote the decay and modification rates for mRNA respectively; while d and e are used for the protein. To simplify the analysis, Zeiser et al. assumed the decay rates (b and d) as being identical for both forms of the mRNA and protein. By fixing b = 0.028, d = 0.031 and under condition of having oscillation period about 120 min, all determined in Hirata et al. [6], Zeiser et al. estimated a set of fitting parameters, displayed in 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' values in 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].

Anthranilate synthase (AS) is the enzyme catalysing the first reaction step in the tryptophan synthesis pathway. The pathway end product tryptophan is fedback to inhibit anthranilate synthase activity if tryptophan level is high. Enzyme inhibition therefore forms the third negative feedback loop in the tryptophan operon system.

(30a)

(30b)

(30c)

We set up a simple three-species model for the tryptophan system as in equations (30). The state variables are the mRNA (M), the AS enzyme (E) and the tryptophan amino acid (T). Each negative feedback loop is modelled using a Hill function; O_{t} is the static total operon concentration; k_{1}, k_{2}, k_{3} represent transcription rate, translation rate and tryptophan synthesis rate, respectively; k_{d1}, k_{d2}, k_{d3} are the degradation rates (aggregated parameters combining the decay rate and dilution rate due to cell growth). Consumption of tryptophan for protein synthesis is simply assumed to follow first order kinetics and represented by the last term of equation (30c). Parameter values are adapted from [47, 48] and tabulated in Table 4.

It can be seen that model (30) is a case of the general multiple-looped considered in the section "Coupled-loop systems". We applied the analysis to the tryptophan operon system for the parameter set in Table 4. Furthermore, for each parameter p in Table 4, we varied it over a wide range around its nominal value (from p × 10^{-2} to p × 10^{2}) while fixing the other parameters. Under all these scenarios, system stability was always obtained. Given the nominal values of the synthetic and degradation rates, the system failed to demonstrate oscillations even at very high 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–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 strength; suggesting that higher cooperativity level tends to enhance the probability of demonstrating sustained oscillations. This result is in line with previous results from Tyson & Othmer [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 suggests 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 oscillation-prone 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-negative-feedback-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 diversity 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.

References

Alon U: An Introduction to Systems Biology: Design Principles of Biological Circuits. 2007, CRC Press

Thomas R, Thieffry D, Kaufman M: Dynamical behaviour of biological regulatory networks – I. Biological role of feedback loops and practical use of the concept of the loop-characteristic state. Bull Math Biol. 1995, 57 (2): 247-276.

Schnarr M, Oertel-Buchheit P, Kazmaier M, Granger-Schnarr M: DNA binding properties of the LexA repressor. Biochimie. 1991, 73 (4): 423-431. 10.1016/0300-9084(91)90109-E

Nguyen LK, Kulasiri D: On multiple regulatory mechanisms in the tryptophan operon system in Escherichia coli: in silico study of perturbation dynamics. In Silico Biol. 2008, 8: 485-510.

Segel LA, : Mathematical models in molecular and cellular biology. 1980, Cambridge [Eng.]: Cambridge University Press

Hoffmann A, Levchenko A, Scott ML, Baltimore D: The IkappaB-NF-kappaB signaling module: temporal control and selective gene activation. Science. 2002, 298 (5596): 1241-1245. 10.1126/science.1071914

Ishiura M, Kutsuna S, Aoki S, Iwasaki H, Andersson CR, Tanabe A, Golden SS, Johnson CH, Kondo T: Expression of a gene cluster kaiABC as a circadian feedback process in cyanobacteria. Science. 1998, 281 (5382): 1519-1523. 10.1126/science.281.5382.1519

Griffith JS: Mathematics of cellular control processes. I. Negative feedback to one gene. J Theor Biol. 1968, 20 (2): 202-208. 10.1016/0022-5193(68)90189-6

Song H, Smolen P, Av-Ron E, Baxter DA, Byrne JH: Dynamics of a minimal model of interlocked positive and negative feedback loops of transcriptional regulation by cAMP-response element binding proteins. Biophys J. 2007, 92 (10): 3407-3424. 10.1529/biophysj.106.096891

Tao Y, Jia Y, Dewey TG: Stochastic fluctuations in gene expression far from equilibrium: Omega expansion and linear noise approximation. J Chem Phys. 2005, 122 (12): 124108- 10.1063/1.1870874

Maheshri N, O'Shea EK: Living with noisy genes: how cells function reliably with inherent variability in gene expression. Annual review of biophysics and biomolecular structure. 2007, 36: 413-434. 10.1146/annurev.biophys.36.040306.132705.

Kulasiri D, Nguyen LK, Samarasinghe S, Xie Z: A Review of Systems Biology Perspective on Genetic Regulatory Networks with Examples. Current Bioinformatics. 2008, 3 (3): 197-225. 10.2174/157489308785909214.

MacDonald N: Time lag in a model of a biochemical reaction sequence with end product inhibition. J Theor Biol. 1977, 67 (3): 549-556. 10.1016/0022-5193(77)90056-X

Yanofsky C: The different roles of tryptophan transfer RNA in regulating trp operon expression in E. coli versus B. subtilis. Trends Genet. 2004, 20 (8): 367-374. 10.1016/j.tig.2004.06.007

Chassagnole C, Fell DA, Rais B, Kudla B, Mazat JP: Control of the threonine-synthesis pathway in Escherichia coli: a theoretical and experimental approach. Biochem J. 2001, 356 (Pt 2): 433-444. 10.1042/0264-6021:3560433

Tyson JJ: Periodic enzyme synthesis and oscillatory repression: why is the period of oscillation close to the cell cycle time?. J Theor Biol. 1983, 103 (2): 313-328. 10.1016/0022-5193(83)90031-0

Takebayashi K, Sasai Y, Sakai Y, Watanabe T, Nakanishi S, Kageyama R: Structure, chromosomal locus, and promoter analysis of the gene encoding the mouse helix-loop-helix factor HES-1. Negative autoregulation through the multiple N box elements. J Biol Chem. 1994, 269 (7): 5150-5156.

Jensen MH, Sneppen K, Tiana G: Sustained oscillations and time delays in gene expression of protein Hes1. FEBS Lett. 2003, 541 (1–3): 176-177. 10.1016/S0014-5793(03)00279-5

Santillan M, Mackey MC: Dynamic regulation of the tryptophan operon: a modeling study and comparison with experimental data. Proc Natl Acad Sci USA. 2001, 98 (4): 1364-1369. 10.1073/pnas.98.4.1364

Bhartiya S, Chaudhary N, Venkatesh KV, Doyle FJ: Multiple feedback loop design in the tryptophan regulatory network of Escherichia coli suggests a paradigm for robust regulation of processes in series. J R Soc Interface. 2006, 3 (8): 383-391. 10.1098/rsif.2005.0103

Goldbeter A: A model for circadian oscillations in the Drosophila period protein (PER). Proc Biol Sci. 1995, 261 (1362): 319-324. 10.1098/rspb.1995.0153

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.

Additional file 1: Supplementary Information. This file consists of three parts. Section 1 presents the mathematical derivations for the results involving the single-loop systems. Section 2 presents the mathematical derivations for the results involving the coupled-loop systems. Section 3 presents the mathematical derivations for the results involving the system with endproduct utilisation. Section 4 gives the explicit expressions of the coefficients of functions f and g discussed in the main text, and some intermediate derivation steps. The supplementary figure S1 is given in section 5. (DOC 358 KB)

Authors’ original submitted files for images

Below are the links to the authors’ original submitted files for images.

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Nguyen, L.K., Kulasiri, D. On the functional diversity of dynamical behaviour in genetic and metabolic feedback systems.
BMC Syst Biol3, 51 (2009). https://doi.org/10.1186/1752-0509-3-51