## Abstract

Feedback controls are central to cellular regulation. Negative-feedback mechanisms are well known to underline oscillatory dynamics. However, the presence of multiple negative-feedback mechanisms is common in oscillatory cellular systems, raising intriguing questions of how they cooperate to regulate oscillations. In this work, we studied the dynamical properties of a set of general biochemical motifs with dual, nested negative-feedback structures. We showed analytically and then confirmed numerically that, in these motifs, each negative-feedback loop exhibits distinctly different oscillation-controlling functions. The longer, outer feedback loop was found to promote oscillations, whereas the short, inner loop suppresses and can even eliminate oscillations. We found that the position of the inner loop within the coupled motifs affects its repression strength towards oscillatory dynamics. Bifurcation analysis indicated that emergence of oscillations may be a strict parametric requirement and thus evolutionarily tricky. Investigation of the quantitative features of oscillations (i.e. frequency, amplitude and mean value) revealed that coupling negative feedback provides robust tuning of the oscillation dynamics. Finally, we demonstrated that the mitogen-activated protein kinase (MAPK) cascades also display properties seen in the general nested feedback motifs. The findings and implications in this study provide novel understanding of biochemical negative-feedback regulation in a mixed wiring context.

## 1. Introduction

Feedback regulations are arguably the most common control mechanisms employed by cellular systems, ranging from metabolic to transcriptional to signalling levels [1,2]. Feedback provides cells with self-adjusting capability essential for their survival in fluctuating environments. Failure of proper feedback regulation can lead to diseased states of the cells [2]. A large body of research has thus been dedicated to studying feedback regulation and the implications at cellular levels. While positive feedback is critical in giving rise to multiple stable steady states, underlining switch-like behaviours and prolonging signals [3], negative feedback is essential in maintaining homeostasis and generating oscillatory dynamics [1,4].

Despite concrete progress, much is still to be understood about the regulation and function of feedback mechanisms in many cellular systems. Advances in molecular techniques are providing better understanding of the mechanistic picture of the cells, which show that feedback mechanisms often appear tangled up and act cooperatively rather than in isolation. Examples of cellular systems in which several negative-feedback loops intertwine to coordinate complex dynamics are ubiquitous. At least three different negative-feedback mechanisms are found to regulate the oscillatory nucleo-cytoplasmic shuttling of the transcriptional nuclear factor kappaB [5]. More than two negative-feedback loops, acting at different layers of the cascade, were also described in the extracellular signal-regulated kinase (ERK) pathway [6], which may underline the oscillations recently observed in certain cell types [7,8]. Furthermore, multiple negative-feedback structures are prominent in the synthesis pathways of many common amino acids, where the amino acids regulate their own production by concurrently inhibiting multiple upstream reaction steps [9,10]. For instance, the tryptophan operon has evolved seemingly redundant negative-feedback loops, but which appear to be providing distinct dynamical and noise-related functions [11–13]. The cooperative feedback structures exemplified above often underline intricate connections between molecular components within cellular pathways as well as crosstalks between them. Breakdown of such concerted operation could lead to malfunctions in the harbouring system and results in diseased state of the cell. Although studies looking at linked-feedback structures, e.g. mixed positive and negative feedbacks, have been carried out to some extent, their number remains limited [14–16]. Structures of mixed negative feedbacks have received far less attention, creating gaps in our understanding concerning their dynamical behaviour.

In the current work, we set out to fill some gaps by studying in detail the dynamical properties of various general dual-negative-feedback motifs with nested structures, using combined analytical and computational approaches. Schematic of the studied motifs and the single-feedback loop counterpart are given in figure 1. Each motif consists of two negative-feedback loops acting simultaneously at different layers of control. All motifs share the outer feedback loop but differ in the second inner loops, three of them are auto-repression feedback loops. We aimed to understand in these coupled-feedback contexts, how the individual loops contribute to the regulation of oscillations’ existence and of the basic quantitative characteristics including the period (frequency), amplitude and average (mean) value (electronic supplementary material, figure S1*a*). Of particular interest, we asked if the coupled-feedback structures provide any enhanced robustness to the control of oscillation. To this end, we derived analytical conditions based on which parameter regimes leading to distinct dynamical characteristics can be obtained. Crossing between these regimes corresponds to qualitative change of the system's behaviours including fixed-point steady state and sustained oscillations, which are the prominent dynamics in the considered motifs. As a specific example, similar questions were investigated in the well-conserved mitogen-activated protein kinase (MAPK) signalling cascades where a nested negative-feedback structure is at work [6,17].

## 2. Models and methods

### 2.1. Feedback motifs and mathematical description

Early work on modelling feedbacks in gene-related systems started with that of Goodwin [18]. The first repressive model of transcriptional regulation where a gene's expression is inhibited by its own protein product was proposed. This simple model has provided a practical framework for a plethora of subsequent studies exploring the dynamics of negative-feedback genetic systems [19–24]. Here, we extended the Goodwin framework to describe feedback systems with multiple negative mechanisms. A multi-loop system is described by the following differential equations:
2.1where *m* is the length of the system, [*S*_{i}], *k*_{i} and *k*_{di} (*i* = 1 … *m*) denote concentrations of the molecular species *S*_{i} and its rates of synthesis and degradation. Hill kinetics are used to model feedback inhibition [25], given by
2.2

In equation (2.1), there is a subindex set *I*_{k} (*k* = 1 … *m*), which describes the inhibition of *S*_{k} via different feedback loops by the species denoted by the subindex *k*. For instance, *I*_{1} = {2, 3} if *S*_{1} is repressed by *S*_{2} and *S*_{3}. Figure 1 displays the scheme of the dual-negative-feedback motifs of length three considered in this study. The ordinary differential equations (ODE) describing their dynamics following the same notations as in equation (2.1) are given in §1 of the electronic supplementary material. The input *S*_{0}, which can be a gene, ligand or signalling protein that triggers sequential activation of a transcriptional or signalling cascade, was included in parameter *k*_{1} in (2.1) for convenience. Degradation and synthesis are assumed to follow first-order kinetics. Since the half-saturating constant *K* is inversely correlated with the inhibiting Hill function, we refer to *K* as the reciprocal feedback strength of the respective loop [25,26]. Biologically, *K* represents the binding constant between the repressor and the repressed target [26]. A small value of *K* indicates that the feedback is strong while large *K* implies weak inhibition. On the other hand, the Hill coefficient *n* modulates the degree of function sigmoidity and signifies the level of nonlinearity exhibited by the feedback loop.

### 2.2 . Analytical analysis of dynamics based on the Routh–Hurwitz theorem

Biochemical systems can exhibit diverse types of dynamical behaviour including fixed point, mono- or bi-stable steady states, oscillations (limit cycles) and even irregular fluctuation (chaos) [27,28]. Sustained oscillations and fixed-point steady states are often the prominent dynamics displayed by negative-feedback systems which are separated by Hopf bifurcation [21,24]. Four motifs in figure 1*a,d*–*f* can be categorized as monotone cyclic negative-feedback systems and thus have been proved to only exhibit fixed point or sustained oscillations [21]. Linear stability analysis hence provides a convenient way to detect oscillations in these systems. For the remaining dual-feedback motifs that are non-monotone cyclic (figure 1*b,c*), further methods are required to validate oscillations, which we will discuss below. We focus here on how the transition between these two dynamics is governed for the considered motifs.

Dynamical insights of a system can be obtained by linearizing the governing differential equations around the equilibrium state and probing the eigenvalues of the resulting Jacobian matrix, *J*, given by [29]:
2.3

If the real parts of all *J*'s eigenvalues are negative, the system's equilibrium is deemed stable and reaches a fixed-point steady state, whereas if any of the real parts are positive, system stability is lost and in the motifs considered here, the system moves into a limit cycle regime where it oscillates. The eigenvalues are usually obtained numerically and their real part values are assessed. To facilitate analytical treatment, we instead employ the Routh–Hurwitz theorem [29,30] which states that the eigenvalues, essentially the roots of the characteristic polynomial |*J* − *λI*| = *λ*^{n} + *α*_{1}*λ*^{n}^{−1} + … +*α*_{n}_{−1}*λ* + *α*_{n} = 0, all possess negative real parts if and only if the following condition holds true
2.4

where Δ_{k}(*k* = 1 … *n*) are the determinants of the properly constructed matrices

Equation (2.4), therefore, represents the sufficient condition for oscillations in the dual-feedback motifs while it is only a necessary condition in motifs *F*_{31,32} and *F*_{31,21}. To verify that (2.4) also constitute the sufficient condition for oscillations in these systems, we employed analysis based on the Floquet theory that globally assesses the stability of the obtained periodic solutions [31–33]. The mathematical and computational background of the Floquet analysis is given in the electronic supplementary material. This method was numerically implemented in Mathematica v. 8.0 (code can be provided on request). In the following sections, we will present first the analytical analysis based on (2.4) for the single-loop system *F*_{31} followed by the exemplified dual-loop motif *F*_{31,32}. Detailed derivations for the remaining motifs are provided in the electronic supplementary material.

## 3. Results

### 3.1. Dynamical analysis of the single-loop system F_{31}

We begin by examining the single negative-feedback motif *F*_{31}, where the first species in the cascade *S*_{1} is inhibited by the last species *S*_{3}. The ODE model of this motif is given in §1, electronic supplementary material. Systems of this type have been previously studied assuming equal degradation rates [24,34], yielding informative results but are limited in their biological implication owing to the oversimplified assumption. Our following analysis assumes no such restricting assumption and was aimed to facilitate an extension to motifs of multiple-feedback loops. According to the Routh–Hurwitz theorem, *F*_{31} reaches the steady state if and only if
3.1where *α*_{1}, *α*_{2}, *α*_{3} are the coefficients of the system's characteristic polynomial *f*(*λ*) = *λ*^{3} +*α*_{1}*λ*^{2} + *α*_{2}*λ* + *α*_{3}. Hereafter, we denote *x*_{i} as the steady-state values (when the system reaches equilibrium) of *S*_{i} (*i* = 1, 2, 3). These steady-state variables of the system can be obtained by setting the right-hand side (RHS) of equation (2.1) to zero which leads to
3.2
3.3

Using equation (3.2), the coefficients *α*_{1}, *α*_{2}, *α*_{3} are given by
3.4

This introduction of a composite variable *M*_{1} (0 < *M*_{1} < 1), as will soon become apparent, allows us to ignore possible complications involving *x*_{3} and the related rational functions such as those in equations (3.2) and (3.3), which greatly facilitates our treatment. Since *α*_{1}, *α*_{2} > 0, equation (3.1) means that oscillatory dynamics arise if *α*_{1}*α*_{2} − *α*_{3} < 0. Substituting α_{1}, α_{2}, α_{3} into this condition gives
3.5

Rearranging *x*_{3} in terms of *M*_{1}, *K*_{1} and *n*_{1} from equation (3.4) and substitute *x*_{3} into equation (3.3), we have:
*K*_{1} can now be expressed as *K*_{1} = *F*(*M*_{1}), where is a function of *M*_{1}.

Because d*F*(*M*_{1})/d*M*_{1} < 0 for *M*_{1}∈[0,1], *F*(*M*_{1}) strictly decreases with *M*_{1}; the oscillations condition, equation (3.5), is transformed into the following equivalent form:
3.6where *f* and *g* involve the synthesis and degradation parameters given in equations (3.3) and (3.5).

The symbolic, parametric condition equation (3.6) governs the bifurcation points separating the dynamics of fixed-point steady state and limit cycles of the single-loop system. It is immediately clear that the system can only oscillate if the feedback loop is sufficiently strong (small enough *K*_{1}). Moreover, the Hill coefficient *n*_{1} must be greater than *g* for any chance of oscillations to realize. The latter observation is consistent with the previous derivation by Griffith [20].

### 3.2. Dynamical analysis of the dual-feedback systems

Demonstrating this, we describe here treatment of the dual-feedback system *F*_{31,32}, while only the main results are reported for other motifs whose detailed derivations are given in the electronic supplementary material. Motif *F*_{31,32} contains a second, inner-feedback loop, where the cascade product *S*_{3} also inhibits species *S*_{2}. The governing ODEs are presented in §1, electronic supplementary material. This newly incorporated loop is also modelled by Hill kinetic function in which the inverse of the parameter *K*_{2} describes the inhibition strength, and *n*_{2} measures the nonlinearity level of the feedback regulation. The steady-state solution of the system variables can be obtained by setting the RHS of the differential equations to zero which leads to
3.7

Introducing and from equation (3.7) we obtain
3.8where *f* is given in equation (3.3). The coefficient *α*_{1}, *α*_{2}, *α*_{3} of the characteristic polynomial in this case can be expressed as
3.9where *α* *= g*/*n*_{1},*β* = (*n*_{2}/*n*_{1})(*k*_{d2} + *k*^{d3})/*k*_{d1}.

Since *x*_{3} can be written as , substitute this into equation (3.8) and rearrangement yields
3.10where *F*(*M*_{1}) was defined previously. Since *F*(*M*_{1}) strictly decreases with *M*_{1}∈[0,1], equation (3.10) becomes
3.11

It is convenient here to report similar conditions obtained for the remaining dual-feedback motifs (see electronic supplementary material for derivation). 3.12 3.13 3.14 3.15

### 3.3. Addition of the inner-feedback loops reduces and potentially eliminates oscillatory dynamics

It can be observed that equations (3.11–3.15) share a common form despite the difference in the coefficients of the *M*_{2}-related terms inside function *F*. Moreover, an extra term (1−*M*_{2}) appears on the RHS of these equations that does not appear in (3.10) derived earlier for *F*_{31}. It turned out that this additional term and those inside the function *F* brings about important implications in how the inner negative-feedback loops affect system dynamics in the dual-feedback motifs, as shown below.

First, since the composite variable *M*_{2} belongs to the [0,1] interval, the following holds true for all dual-feedback motifs:
3.16

Note that the form of *β* are motif-specific and *γ* = 0 for *F*_{31,21} and *F*_{31,32}. Equation (3.16) indicates that *K*_{1} < *F*(*g*/*n*_{1}) is a necessary condition for oscillations in the dual-feedback motifs which, in comparison with condition equation (3.6), implies that addition of any of the negative-feedback loops: *S*_{3} repressing *S*_{2}, *S*_{2} repressing *S*_{1}, or auto-repression of *S*_{1}, *S*_{2} or *S*_{3} would not produce oscillations if the *F*_{31} does not oscillate in the first place. Equation (3.16) also suggests that the dual-feedback structures reduce the effective upper-bound threshold of *K*_{1} possible for oscillations, thereby shrinking the oscillations domain compared with the single-loop motif.

Secondly, consider the case when the second-feedback loop is weak. In this scenario, *K*_{2} is large that makes *M*_{2} ≈ 0. The RHS of equations (3.11–3.15) then become *F*(*α*). The condition for oscillations in the coupled-loop motif *F*_{31,32} thus approaches that of the single-loop motif *F*_{31} in the limiting case of weak feedback inhibition by the inner loop. On the other hand, when the inner loop is strong, *K*_{2} is large resulting in *M*_{2} ≈ 1. The RHS of equations (3.11–3.15) then become null which would lead to unsatisfactory conditions for oscillations if *K*_{1} > 0. An important and rather unexpected consequence is that the inner-feedback loop, when sufficiently strong could completely eliminate oscillations by overriding the oscillatory dynamics triggered by the outer feedback loop. To verify these predictions, we performed numerical simulations of all the motifs with arbitrary parameter sets, which indeed showed that pre-existing oscillations in the single-loop system *F*_{31} can be abolished when incorporating a sufficiently strong second-feedback loops; whereas oscillations are maintained for weak inhibition levels (electronic supplementary material, figure S2).

### 3.4. Dynamics partitioning revealed differential oscillations-reducing efficiency by the dual-feedback motifs

Although equation (3.11–3.15) allowed us to draw new insights into the effects of the nested-feedback structures on the appearance likelihood of oscillations, the presence of *M*_{2} in these equations complicates the derivation of conditions that involve only the model parameters. Below, we further refine these equations to achieve this aim.

Again, we use the system *F*_{31,32} as an illustrative case. Since *x*_{3} can also be written as , equation (3.8) enables us to express *M*_{1} in terms of *M*_{2} as follows
3.17where is a function of *M*_{2}. Note that function *F*_{1} is similar to *F* defined in §3.1 except for the parameter *n*_{2}. Equation (3.10) can then be rewritten as
3.18where *F*_{2}(*M*_{2}) = *F*(1−*K*_{2}/*F*_{1}(*M*_{2}))(1−*M*_{2}).

Since *F*_{2}(*M*_{2}) strictly increases with *M*_{2} owing to d*F*_{2}(*M*_{2})/d*M*_{2} > 0 for *M*_{2}∈[0,1] (see the electronic supplementary material,), derivation of the parameters-only condition that governs oscillatory dynamics for the considered motif is simplified to determination of values of the variable *M*_{2} satisfying equation (3.10) or the equivalent form
3.19

where *G*(*M*_{2}) = 1−*K*_{2}/*F*_{1}(*M*_{2}).

Such ranges of *M*_{2} cannot be obtained analytically but can be quickly obtained by numerically solving equation (3.19), and can be conveniently visualized on the *M*_{1} versus *M*_{2} coordinate where intersection of the functions *G*(*M*_{2}) and *α* + *βM*_{2} are shown (figure 2*a*). Plugging the obtained numerical value of *M*_{2} into equation (3.18) allows us to readily compute the value of *K*_{1} that result in oscillations given any parameterization of the remaining model parameters. Bifurcation plots that partition the parameter space (i.e. *K*_{1} versus other parameters) can therefore be constructed and based on these, we will be able to perceive how changes in certain parameters affect system dynamics and bring about oscillations. It should be noted that the form of the function *F*_{2} is slightly different across the motifs, as detailed in the electronic supplementary material.

Figure 2*b* compares the partition of systems dynamics (fixed-point steady state and oscillations) on the *K*_{1} versus *K*_{2} coordinate for all five dual-negative-feedback motifs. Consistently in all these motifs, a second loop with strong inhibition, signified by small *K*_{2}, tends to bring the system out of the oscillation domain. In fact, we observed an increasing effect where stronger second loops result in smaller size of the oscillation domain. Since the same value of *n*_{2} was used for all motifs in figure 2*b*, it reveals that the ability to diminish oscillations is weakened as the second loop moves upstream of the cascade. Indeed, the addition of *S*_{3}'s auto-inhibition loop gives the smallest oscillations domain, whereas it is largest with auto-inhibition of *S*_{1}. Furthermore, the figure shows, as predicted, that *K*_{1} approaches a threshold value (*F*(*g*/*n*_{1})) at large *K*_{2}, which is independent of *n*_{2}.

### 3.5. Dependence of oscillations on the degradation rates and optimization of oscillation's appearance

The rates of degradation have been suggested to play important roles in influencing the existence of oscillations in feedback systems [35–37]. Symmetry among elements within a feedback loop was predicted to increase the occurrence likelihood of oscillations [35]. It is, therefore, interesting to study how the degradation rates affect systems dynamics in our dual-feedback context. Figure 3 shows on the *K*_{1} versus *k*_{di} bifurcation diagrams (*i* = 1 … 3) that the threshold *K*_{1} depends biphasically on the degradation rates in all motifs. The bell-shaped oscillation domains suggest that oscillatory behaviour arises only over a bounded range of the degradation rates. The three-dimensional bifurcation plot (figure 3*d*) of *K*_{1} in response to simultaneous changes in pairs of the degradation rates further shows a cone-shaped oscillation domain. The oscillation domain of all the motifs projected in the space of the degradation rates would therefore form a restricted volume, whose size is expectedly diminished as the second-feedback loop's inhibition intensity increases. The ability to oscillate is only realized when the rates of degradation are tuned within this enclosed domain, conferring oscillatory dynamics much less probable in the considered motifs compared with steady-state equilibrium dynamics. The biphasic dependence of the threshold *K*_{1} on the rates of degradation also implies an optimal strength of the *S*_{3}-repressing-*S*_{1} feedback inhibition at which the appearance of systems oscillations is maximized. Importantly, this optimal state is achieved at intermediate levels of the turnover rates of the components within the loop.

Next, we ask the question, what distribution of the degradation rates would lead to maximal oscillations given a fixed steady-state level of the output *S*_{3}? This comparative criterion is realistic since cellular pathways typically maintain a functional range of key proteins. For the single-loop motif *F*_{31}, constant *x*_{3} can be brought about by fixing the product of the rates of degradation (see equation (3.3)). We found that in this case, oscillations are most likely under conditions of equal degradation rates, which agrees with previous prediction [35] (see proof in §2.1, electronic supplementary material). For the dual-feedback systems, *F*_{31,32} and *F*_{31,33}, constant steady-state *S*_{3} levels can be similarly guaranteed by fixing the product *k*_{d1} *k*_{d2} *k*_{d3}. We proved analytically that oscillations are most likely for *F*_{31,32} when *k*_{d2} = *k*_{d3} at any given *k*_{d1}. On the other hand, *F*_{31,33} is most prone to oscillations when *k*_{d1} = *k*_{d2} at a given *k*_{d3} (§2.2, 2.3, electronic supplementary material). To verify these conclusions, e.g. for S_{31,32}, we generated numerous parameter sets that differ only in *k*_{d2}, *k*_{d3} but keep their product fixed. The corresponding value of the threshold *K*_{1} was then computed for each set (electronic supplementary material, figure S2*c*) which indeed shows that the threshold *K*_{1} is maximized at comparable levels of *k*_{d2} and *k*_{d3}. Interestingly, similar analysis shows that *K*_{1} is also maximized at comparable {*k*_{d1}, *k*_{d2}} (electronic supplementary material, figure S2*a*) and {*k*_{d1}, *k*_{d3}} (electronic supplementary material, figure S2*b*), suggesting that oscillations in *F*_{31,32} are most probable under equal degradation rates of all species (*k*_{d1} = *k*_{d2} = *k*_{d3}). This is also the case for the system *F*_{31,33}, as shown in electronic supplementary material, figure S2*d*–*f*. For the dual-feedback motifs *F*_{31,21}, *F*_{31,22} and *F*_{31,11}, constant steady-state *S*_{3} is brought about either by fixing {*k*_{d3}, *k*_{d1}*k*_{d2}}, {*k*_{d3}, *k*_{d1}*k*_{d2}} or {*k*_{d1}, *k*_{d2}*k*_{d3}}. We found that oscillations are most likely for *F*_{31,21}, *F*_{31,22} and *F*_{31,11}, when *k*_{d1} = *k*_{d2}, *k*_{d1} = *k*_{d2} and *k*_{d2} = *k*_{d3}. Analytical proofs of these predictions are shown in §2.4 and 2.5 of the electronic supplementary material, and also confirmed by computational simulations (figure 4 and electronic supplementary material, S3). These results suggest that in systems with nested-feedback structure, symmetry among the elements of the loops also promotes systems oscillations.

### 3.6. Higher nonlinearity of both feedback loops enhances oscillations

Negative-feedback regulation requires a sufficient degree of nonlinearity for oscillations to arise [38]. This is evident from the condition *n*_{1} ≥ *g* needed for oscillations in the single-loop and in the tested dual-feedback systems. Since in the dual-feedback motifs, the feedback strength of the individual loops brought about opposing effects on systems oscillations, we asked if the loops’ nonlinearity levels behave in the same way. Figure 5 shows that this is not the case. In all the tested dual-feedback motifs, increasing the Hill coefficients *n*_{1} and *n*_{2} both enlarged the oscillations domain projected here on the *K*_{1} versus *K*_{2} coordinate, leading to enhanced oscillations. However, the underlying mechanisms of this enhancement are different. While higher *n*_{1} supports oscillations by raising the threshold *K*_{1}, allowing oscillations possible at an even weaker *S*_{3}-repressing-*S*_{1} loop strength, higher *n*_{2} enables oscillations at even stronger intensity of the inner loops. Taken together, the result implies that although addition of stronger inner loops diminished oscillations, raising its nonlinearity level tolerated this effect and thus promoted oscillations instead.

### 3.7. Regulation of quantitative properties of oscillation dynamics in the dual-feedback structures

An oscillatory dynamic bears a number of basic quantitative characteristics such as its frequency (period), amplitude (difference between peak and trough) and the average value (mean; electronic supplementary material, figure S1*a*). Recent studies report that these metrics may play important roles in directing cellular responses. In yeast, the level of extracellular calcium was found to modulate the nuclear translocation frequency of transcriptional factor Crz1 and this frequency directly controls expression of multiple target genes [39]. In calcium and endocrine oscillations, the frequency was also thought to contain key functional signal [4,40,41]. In other cases, the amplitude plays key physiological roles. Allada *et al*. [42] showed that a clock gene mutant in *Drosophila melanogaster* that produces persistent oscillations but with reduced amplitude led to behavioural malfunction. Previous computational work on mixed positive-negative feedback oscillators revealed such coupling enabled robust tuning of the frequency at stable amplitude [16]. It is therefore of our interest to investigate how the three metrics (period, amplitude and mean) of the oscillations dynamics are controlled in the mixed negative-feedback motifs.

First, to understand how the individual loops affect the oscillation pattern in the dual-feedback motifs, the feedback strength of either loop was varied and the period, amplitude and mean value of the oscillatory level of the species S_{3} were computed (see §4.1, electronic supplementary material, for computational procedure). The parameter values used in our calculation was based on the physiological values obtained for the endogenous circadian clock system in *Neurospora crassa* [37]. The molecular basis of this system is described in electronic supplementary material, figure S1*b*, which follows the general motif *F*_{31}. Our computational result showed that the period increases with stronger feedback of the outer loop, while decreases towards zero (i.e. loss of oscillations) at higher feedback of the nested loop (figure 6*a*). This was observed for all the dual-feedback motifs which is consistent with our previous finding that the second-feedback loops can eliminate oscillations when sufficiently strong. We observed a biphasic dependence of the amplitude on the outer loop, indicating existence of an optimal strength a maximal amplitude (figure 6*b*). On the other hand, the inner loops decrease the amplitude across all motifs as expected. Also as expected, both loops decrease the average value of the oscillatory *S*_{3} concentration (figure 6*c*).

Next, we simultaneously varied the strength levels of both the outer and inner loops and recorded the period, amplitude and mean value. Plotting one metric against another in parametric planes revealed interesting observations (figure 6*d–f*). We found that in all the dual-feedback motifs, altering the short inner loop allowed the system to widely modulate the period or amplitude while keeping the mean value at near-constant levels, indicated by the almost perfect horizontal lines in figure 6*e*,*f*. On the other hand, varying the long outer loop allowed the systems to modulate the period while keeping the amplitude relatively stable. This effect is increased when the systems are well within the oscillating regime, indicated by the lines inside the box in figure 6*d*. Furthermore, higher inner loop seemed to make the tunability more pronounced (lower lines, figure 6*d*). These observations suggest that having coupled negative feedbacks provides the systems with robust ways to adjust different features of the oscillation dynamics, an ability similar to that of the mixed positive–negative motifs [16]. Such robustness may be the key in enabling the systems to decouple functional controls regulated by different features of the oscillation dynamics. To see to what extent the above findings are dependent on parameterization, we carried out similar analysis for other random parameter sets generated from the basal set (electronic supplementary material, figure S5*a* and §3). We found that the findings still hold true for these sets, suggesting that they are likely consequences of the motifs’ topology.

### 3.8. The mitogen-activated protein kinase signalling cascade: an example

The MAPK signalling cascades are among the best-studied signalling systems in cells, which play crucial roles in many cellular processes including proliferation and survival [6,17]. Kholodenko [43] predicted that sustained oscillations can occur in these cascades under physiological parametric conditions. This prediction has been backed up by recent experimental reports [8,44]. The MAPK cascade model analysed by Kholodenko [43] however contained only a single negative-feedback loop caused by inhibition of the active form of the upstream GTPase (e.g. RasGTP) by the activated MAPK protein, while infact multiple negative-feedback loops are at work [6]. To further understand oscillation dynamics of the MAPK cascades in multi-loop context, we analysed an extended MAPK model. We considered an additional negative feedback brought about by inhibition of the MAPK kinase kinase (MKK) by the active MAPK protein (ppMAPK, figure 7*a*). For instance, in the extracellular signal-regulated kinases (ERK) pathway, an exemplary MAPK cascade, the second negative feedback is caused by phosphorylation of Raf-1 on inhibitory sites by active ERK [6,17].

The ODE and parameter values of the extended Kholodenko model are given in electronic supplementary material, tables S1 and S2. The feedback strength of the two loops can be adjusted by varying parameters *K*_{I1} and *K*_{I2}. Our numerical stability analysis showed that these feedback loops affect oscillations differently (figure 7*b*). Higher strength of the outer loop tends to move the system into the oscillatory domain, thereby enhancing oscillations, while a stronger inner loop tends to move the system out of the oscillatory domain, thereby weakening oscillations (figure 7*b*). Similar to the general dual-feedback motifs, a sufficiently strong inner loop can destroy oscillations in the MAPK cascades. However, unlike in the general motifs, where the outer loop can always rescue oscillations, in the MAPK cascades, there exists a barrier strength of the inner loop over which the system fails to oscillate regardless of how strong the outer loop is. This difference may be due to the different wiring and level of complexity in the MAPK model. We further computed the period, amplitude and mean value of the oscillatory ppMAPK levels in response to increasing strength of the outer and inner loop (electronic supplementary material, figure S6). Interestingly, both loops reduced the oscillations period and mean value. The outer loop also exhibited a bi-phasic modulation of the amplitude as expected, while the inner loop consistently decreased the amplitude.

## 4. Discussion

Understanding the emergence of functional behaviour from the dynamical controls of cellular networks is of fundamental importance in biology. It is, therefore, critical to dissect how certain system dynamics are brought about and how they are regulated. Oscillations are ubiquitously observed in cells and have been shown to play important roles in a broad array of cellular processes [45,46]. The oscillatory behaviour produced can be time-dependent, as in glycolytic oscillations [47], the circadian rhythms [48], the cell cycles [49], or space-dependent, as in calcium waves [50,51] and neuronal oscillations [52]. In all of these examples, oscillations are driven by some sort of negative-feedback mechanism. However, in many oscillating systems, their circuits contain more than one negative-feedback loop which begs the question how they might cooperate to control oscillations. In the current work, we studied the emergence and control of oscillation dynamics in the context of coupled, nested negative feedback biochemical motifs. By comparative examination of the dual-loop and single-loop structures, we aimed to dissect the relevant roles of the individual feedback loops. To this end, we employed combined analytical and numerical methods.

Most interestingly, we found that the individual-feedback loops regulate oscillations in opposing manners in all of the general dual-feedback motifs. The long, outer loop promotes oscillations whereas the short, inner loops strictly inhibit oscillations. Such oscillations-inhibiting effects, when sufficiently strong, can dominate the oscillation-generating effect of the outer loop and completely destroy oscillations. Coupling multiple negative feedbacks together thus can suppress rather than enhance oscillations as found previously [53]. Note that other ways of negative-feedback coupling such as mutual inhibition can also inhibit oscillations and promote bistability instead [3], as they behave like a positive feedback. Such wiring is, however, structurally different to the nested-feedback motifs considered here, highlighting that different ways of coupling may yield similar system responses yet through distinct mechanisms. Our analytical proof further showed that the above observation is parameter independent and rather a feature of the network topology. Comparison of the oscillations domain among the dual-feedback motifs under the same parametric specification of the inner-feedback loop revealed that as it moves further downstream of the cascade, it represses oscillations more strongly. Auto-inhibition of the final component of the cascade, therefore, reduces oscillations most efficiently.

Having showed that the shorter loop of the dual-feedback structures can abolish oscillations driven by the longer loop, it is useful to understand this in an intuitive way. Negative feedback requires sufficient delay and nonlinearity to give rise to oscillations [38]. Our simulations show that removing the outer loop from the dual-feedback motifs rendered them incapable of exhibiting oscillations regardless of parameter choice (data not shown). The failure of these short loops to maintain oscillations on their own may be due to the lack of delay, a cause that is likely to also underline their ability to suppress oscillations in the coupled structures. To test this idea, we constructed variant motifs of the system *F*_{31,11} which extend the auto-inhibition loop of *S*_{1} by one or two intermediate components (electronic supplementary material, figure S1*c*–*d*). Addition of one intermediate component still reduced oscillations; while having two intermediates introduced enough delay to prevent the extended loop from having an oscillations-reducing property. One important consequence then follows regarding the modelling of feedback systems in general. Models of biochemical systems are in many cases greatly simplified. Quite often, multiple activation steps occurring via intermediate components are described as a single step for convenience without explicitly taking into account the time delay. We, however, showed that description of a negative feedback with one or more intermediate components could result in markedly different effects of the loop towards system dynamics. Time delay, if significant, thus should be considered in the model for faithful representation of the system in reality.

It also follows from our results that motifs formed by combining the outer feedback loop and any combination of the secondary inner loops would behave similarly to the dual-feedback motifs in their ability to suppress oscillations. The opposing roles of the coupled feedbacks suggest intriguing implications regarding the evolution of negative-feedback structures. In principle, changes in system dynamics such as loss of oscillatory behaviour can be brought about by mutations that break the oscillation-generating feedback loop, or alternatively by the emergence of new feedback loops with oscillation-abolishing property. Under certain scenarios, the latter strategy may be more advantageous than the former, for example, when the newly evolved loops provide additional benefit to the cells. Since measurements like period, amplitude and mean value are intrinsic features of an oscillatory dynamics, we went on to examine how these metrics are regulated by the feedbacks interplay in the linked-feedback motifs. We found that the inner loop enabled flexible modulation of the period or amplitude, while keeping the mean value at near-constant levels. More interestingly, in presence of strong inner loop, varying the outer loop allowed the systems to significantly change the period, while keeping the amplitude relatively stable. Such features were not possible in the single-loop motif. Our findings suggest that having nested-feedback structures provides robust controls of oscillation dynamics and its basic quantitative characteristics. Such robustness may have underlined the evolution of some biochemical oscillators possessing these topological structures. To further our investigation in a more biologically realistic system, we asked similar questions in the MAPK signalling cascades that were known to regulate important cellular processes and contain multiple negative-feedback loops in a nested structure. To this end, we extended an existing mathematical model for the dual-feedback MAPK cascade. Our analysis showed that despite the fact that the MAPK model is much more intricate, its negative feedbacks also exhibit opposing regulation of oscillations as in the general motifs. It would be interesting to test this prediction experimentally in the future.

Parameter-sensitivity analysis is of great importance in understanding biological networks. The derivation of the analytical conditions that govern the occurrence of oscillations have allowed us to construct useful bifurcation plots based on accurate knowledge of when switching between different dynamical regimes occurs. Such knowledge is potentially useful in experimental design. It can also aid in the engineering of synthetic circuits with desirable targeted dynamical behaviour, since one could effectively choose appropriate points for perturbation to attain desired dynamics. In particular, we observed that oscillations only occur over limited ranges of degradation rates. The fact that each of these rates must be tuned within its viable range for oscillations to occur suggests that oscillatory behaviour is a strict parametric requirement. Its emergence cannot be guaranteed by the presence of negative feedback alone but depends greatly on the evolutionary tuning of the stability of the systems components, which may imply why oscillations are not observed in all systems with negative-feedback regulations. We went on to examine which distribution of the degradation rates would maximize oscillation occurrence, given they are required to maintain a constant steady-state level of the cascade output. We showed that for the single- as well as the dual-feedback motifs, symmetry among the elements within the feedback loops tends to promote oscillations, consistent with previous reports [35].

In this paper, we employed deterministic modelling as the approach of choice and thus molecular fluctuations (noise) were neglected. However, these fluctuations are inherent features of biochemical networks and may be important in functioning of the cells [54]. Negative feedbacks have been known to suppress noise [55,56], understanding the behaviour of noise in nested feedback structures is an interesting topic. In our recent work [12], we explored how noise is propagated and managed in the prokaryotic tryptophan operon system of *E. coli*, a model system governed by three nested negative-feedback loops. Interestingly, it was found that the individual loops in this system do not all reduce noise, but some actually enhance noise levels of the system output. This work was, however, only concerned with noise when the studied system is outside the oscillating domain. Therefore, we further examined how noise might affect oscillations in the general dual-feedback motifs considered in the current work. Stochastic simulations showed that noise can extend the Hopf bifurcation point in both single- (*F*_{31}) and dual-feedback systems (*F*_{31,32}), producing oscillations at parameter values that does not lead to oscillations in a deterministic setting. Comparing the bifurcation points in *F*_{31} and *F*_{31,32} at different levels of the inner loop showed that the points were not altered, suggesting that addition of the second loop in the nested arrangement had no effect on the stochastic bifurcation points. Further understanding concerning the interplay between noise and oscillations in coupled-feedback settings should be among the focuses of future systems biology studies.

Recent progress in systems biology relies a great deal on numerical analysis and simulations. However, two important bottlenecks are identified with this approach: the lack of reliable parameter values and the limited exploring capability of numerical simulations. The former hurdle can be overcome by development of increasingly advanced experimental techniques; while one way to help reduce the second limitation is to develop novel analytical methods for quantitative analysis of cellular systems. Analytical tools can surmount brute-force simulations and boost the capacity to explore systems dynamics over much wider parameter space. The findings and implications in this study, as a result of such approach, have provided novel understanding of negative-feedback regulation in coupling context.

## Acknowledgements

Systems Biology Ireland would like to acknowledge Science Foundation Ireland's funding award–06/CE/B1129. We thank Prof. Walter Kolch and Boris Kholodenko for stimulating discussion.

- Received January 13, 2012.
- Accepted February 23, 2012.

- This journal is © 2012 The Royal Society