 Research
 Open Access
 Published:
Analytical study of robustness of a negative feedback oscillator by multiparameter sensitivity
BMC Systems Biology volume 8, Article number: S1 (2014)
Abstract
Background
One of the distinctive features of biological oscillators such as circadian clocks and cell cycles is robustness which is the ability to resume reliable operation in the face of different types of perturbations. In the previous study, we proposed multiparameter sensitivity (MPS) as an intelligible measure for robustness to fluctuations in kinetic parameters. Analytical solutions directly connect the mechanisms and kinetic parameters to dynamic properties such as period, amplitude and their associated MPSs. Although negative feedback loops are known as common structures to biological oscillators, the analytical solutions have not been presented for a general model of negative feedback oscillators.
Results
We present the analytical expressions for the period, amplitude and their associated MPSs for a general model of negative feedback oscillators. The analytical solutions are validated by comparing them with numerical solutions. The analytical solutions explicitly show how the dynamic properties depend on the kinetic parameters. The ratio of a threshold to the amplitude has a strong impact on the period MPS. As the ratio approaches to one, the MPS increases, indicating that the period becomes more sensitive to changes in kinetic parameters. We present the first mathematical proof that the distributed timedelay mechanism contributes to making the oscillation period robust to parameter fluctuations. The MPS decreases with an increase in the feedback loop length (i.e., the number of molecular species constituting the feedback loop).
Conclusions
Since a general model of negative feedback oscillators was employed, the results shown in this paper are expected to be true for many of biological oscillators. This study strongly supports that the hypothesis that phosphorylations of clock proteins contribute to the robustness of circadian rhythms. The analytical solutions give synthetic biologists some clues to design gene oscillators with robust and desired period.
Background
Robust oscillations are ubiquitous in biology such as circadian rhythms and cell cycles [1–4]. Robustness is the ability to resume reliable operation in the face of different types of perturbations: environmental and genetic changes, parameter uncertainty, and stochastic fluctuations [5–8]. It is critically important to understand the mechanisms by which biological oscillators robustly work in everfluctuating environments.
To quantify biochemical systems' robustness to fluctuations in all the kinetic parameters, multiparameter sensitivity (MPS) was used [8]. Although MPS is given by the sum of the squared singleparameter sensitivities, it represents how fragile the system's output is when small, random, and simultaneous fluctuations are provided to all kinetic parameters [8]. MPS is mathematically equal to the normalized variance calculated by the Monte Carlo method [9–11]. Use of MPS has revealed that negative feedback loops with multiple phosphorylations produce oscillators robust to parameter uncertainty [8]. The dual feedback model, a simplified version of the Drosophila PERTIM feedback model [12], was found to be the most robust and entrainable among many feedback models with various connection logics [13].
Both numerical and analytical approaches are important for an understanding of design principles underlying robust biological oscillators. A number of numerical studies have been reported for the biological oscillators and many of mathematical models are available from databases such as JWS Online [14], BioModels [15, 16] and BioFNet [17]. Numerical analyses have revealed the mechanisms of how a variety of feedback structures produce robust oscillators [8, 13] and identified critical kinetic parameters that are related to changes in the period and amplitude [9, 11]. Although numerical simulations are useful for a quantitative understanding of how complicated biological oscillators behave, they do not provide explicit information on how the period, amplitude and their robustness depend on kinetic parameters.
On the other hand, analytical solutions directly link the mechanisms and kinetic parameters to dynamic properties such as period and amplitude. However, analytical studies for biological oscillators are scarce compared to numerical counterparts. Most of analytical studies [18–20] focused on whether their models oscillate, but not on the robustness of period and amplitude. Kut et al [21] provided the analytical expressions for the period and amplitude for ElowitzLeibler repressilator [22] and BarkaiLeibler circadian clock [23, 24]. Their work is a great step toward an understanding of how the dynamic properties depend on kinetic parameters. However, the models they analyzed lack generality. Their analytical solutions are not applicable to other biological oscillators. To our knowledge, few analytical solutions for the period and amplitude have been reported. Analytical solutions for a general model of negative feedback loops, which are common structures to biological oscillators [25, 26], greatly contribute to an understanding of design principles underlying robust biological oscillators.
In this paper, we present analytical solutions for the period and amplitude, and their associated MPSs for a general model of negative feedback oscillators. The analytical solutions are in agreement with numerical solutions of their ordinary differential equations. Our analytical or theoretical study of MPSs reveals the mechanisms by which negative feedback loops in biology generate robust oscillations. We present the first mathematical proof that long negative feedback loops make oscillators robust to parameter fluctuations by the distributed timedelay mechanism.
Results and discussion
Negative feedback oscillator model
We consider a general model of negative feedback oscillators with arbitrary loop length, where the feedback is imposed by the last species of the cascade on the first species. Figure 1 shows a schematic diagram of the negative feedback oscillator model. The dynamics is described by a set of differential equations:
where ${x}_{i}$ is the concentration of the i th molecular species (e.g., transcription factor, specifically modified form of protein, metabolite, etc.), ${\alpha}_{i}$ is the decay rate constant, ${\beta}_{i}$ is the production rate constant, and ${K}_{i}$ is the threshold for turning on/off the production of the target molecular species ($i\in \left\{1,2,\dots ,n\right\}$). ${\alpha}_{i}$, ${\beta}_{i}$ and ${K}_{i}$ take positive values. n is the number of molecular species or indicates feedback loop length. $\theta $ is the unit step function given by
where $a$ and $b$ are arbitrary positive values. Introducing the step function makes dynamic models tractable [27]. This on/off behavior arises from a sigmoidal Hill function. ${K}_{i}$ corresponds to the concentration of a regulator molecule at which the production rate of a target molecule reaches a half maximum in Hill function. ${x}_{i}$ takes a positive value in the range of (0, ${\beta}_{i}/{\alpha}_{i}$). The negative feedback oscillator model can produce oscillations when both $n\ge 3$ and ${K}_{i}<{\beta}_{i}/{\alpha}_{i}$ ($i\in \left\{1,\dots ,n\right\}$) are satisfied. It has been proved that this type of negative feedback models cannot generate a sustained oscillation when $n<3$[20]. An example of the dynamics is shown in Figure 2.
Since at least one negative feedback loop is necessary for any biochemical networks to produce oscillations [25], the above model is found everywhere as a network motif in biological oscillatory networks. For example, circadian clock networks can be described by assuming that ${x}_{i}$ are clock genes or specifically modified forms of clock proteins. Our model can also be found in the field of synthetic biology. ElowitzLeibler repressilator [22] and Atkinson's genetic circuit [28] can be considered as slightly modified versions of our model.
Analytical solutions for period, amplitude and their associated MPSs
The period and amplitude were symbolically solved as follows. First, we divided a cycle of oscillations into time intervals (${I}_{ij}$ in Figure 2). Next, we obtained the interval times, peaks and troughs. Finally, we connected all the intervals to calculate the period and amplitude (for details, see Methods). The period $\tau $ and the i th species amplitude ${\epsilon}_{i}$ for the oscillatory model of Eq (1) are given by
where ${\gamma}_{i}$ and ${\delta}_{i}$ are
${C}_{ij}$ ($i\in \left\{1,\dots ,n\right\}$, $j\in \left\{1,2\right\}$) are the integral constants:
Although ${C}_{ij}$ are hard to solve in symbolic form, they can be determined by numerically solving the system of Eqs (7)(8).
Assuming that the peak and trough of ${x}_{i}$ are ${\beta}_{i}/{\alpha}_{i}$ and zero, respectively, the integral constants can be determined without any numerical computations, and the analytical solutions for the period $\tau $ and the i th species amplitude ${\epsilon}_{i}$ are given by (for details, see Methods)
where ${\rho}_{i}$ is the ratio of the threshold ${K}_{i}$ to the amplitude ${\beta}_{i}/{\alpha}_{i}$:
Multiparameter sensitivity (MPS) represents how fragile system's properties are to fluctuations in kinetic parameters (for details, see Methods). The period MPS ${\text{\Phi}}_{\tau}$ and amplitude MPS ${\text{\Phi}}_{{\epsilon}_{i}}$ are given by
Validation of the analytical solutions
First, we like to validate the analytical solutions given by Eqs (3)(4) (which are free of the assumptions of ${\beta}_{i}/{\alpha}_{i}$ peak and zero trough). In order to determine ${C}_{ij}$, fsolve function of MATLAB (The MathWorks, Inc.) was used. The MPSs for period and amplitude were numerically calculated by providing a small perturbation to the values of ${\alpha}_{i}$, ${\beta}_{i}$ and ${K}_{i}$ (see Methods). The computed period, amplitude and MPSs are referred to as the semianalytical solutions because numerical methods were partially used. The semianalytical solutions were compared with the numerical integration solutions. The numerical integration solutions for the period and amplitude were obtained by numerically integrating Eq (1). Table 1 summarizes the methods to obtain the semianalytical and numerical integration solutions. We assigned uniform random values over (0, 1) to ${\alpha}_{i}$ and ${\beta}_{i}$, and those over (0, ${\beta}_{i}/{\alpha}_{i}$) to ${K}_{i}$. The semianalytical solutions were consistent with the numerical integration solutions (Figure 3). Thus, the analytical solutions given by Eqs (3)(4) were validated.
In addition to the correctness of the solutions, the semianalytical approach is computationally more efficient than the numerical integration approach. When $n=3$, the semianalytical approach (fsolve function of MATLAB) required 0.1 sec per parameter set in order to calculate the period, amplitude and their associated MPSs. On the other hand, the numerical integration approach (ode15s function of MATLAB) required 12 sec per parameter set. Although it looks computationally inexpensive to numerically integrate the differential equations of Eq (1), it takes much computational cost because of 'stiffness' that comes from the step function of Eq (2).
Next, we like to validate the analytical solutions given by Eqs (9)(10) and Eqs (12)(13) (which were derived by assuming ${\beta}_{i}/{\alpha}_{i}$ peak and zero trough). The solutions given by Eqs (9)(10) and Eqs (12)(13) are referred to as the full analytical solutions because they are free of any numerical methods (Table 1). We compared the full and semianalytical solutions. We assigned random values to kinetic parameters as we did for the comparison between the semianalytical and numerical integration solutions. The results are shown in Figure 4. The full and semianalytical solutions are consistent especially when the feedback loop is long, because the assumptions (the peak is ${\beta}_{i}/{\alpha}_{i}$ and the trough is zero) are met when the feedback loop is long (Figure 5, and compare it with Figure 2). Even for the short loop ($n=3$), the differences between the full and semianalytical solutions are less than an order of magnitude for most parameter sets. Note that, in Figure 4D, the full analytical solutions for the amplitude MPS are always two, independent of kinetic parameter values (Eq (13)). Although some symbols seem scattered away from the diagonal line in Figure 4D, most of the symbols are so concentrated on the diagonal line that the full and semianalytical solutions are consistent. In summary, the analytical solutions given by Eqs (9)(10) and Eqs (12)(13) were validated.
Effect of changes in kinetic parameters on period and its MPS
Having validated the analytical solutions given by Eqs (9)(10) and Eqs (12)(13), we investigated how the period and period MPS depend on the values of kinetic parameters (we do not investigate the amplitude and amplitude MPS because how they depend on kinetic parameters is clear from Eqs (10) and (13)). Figure S1 shows that how the period depends on kinetic parameters. For simplicity, we assumed $n=3$, ${\alpha}_{i}=\alpha $, ${\beta}_{i}=\beta $, ${K}_{i}=K$ ($i\in \left\{1,2,3\right\}$), and $\rho =K\alpha /\beta $. When the period is given as a function of $\alpha $, it reaches the minimal values at $\rho =0.7228$ (Figure S1AB in Additional file 1). When the period is given as a function of $\beta $ or K, it reaches the minimal values at $\rho =0.5$ (Figure S1CDEF). Note that the negative feedback oscillator model produces oscillations only when $K<\beta /\alpha $ ($\rho <1$). Figure S2 (Additional file 1) shows how the period MPS depends on kinetic parameters. The period MPS always reaches the minimal value 0.2222 at $\rho =0.5959$, implying that the period MPS depends solely on $\rho $ (when n is fixed).
Assuming ${\alpha}_{i}=\alpha $ and ${\rho}_{i}=\rho $ ($i\in \left\{1,2,\dots ,n\right\}$), Eq (12) reduces to
where
$f\left(\rho \right)$ reaches the minimal value 0.6667 at $\rho =0.5959$, and $\underset{\rho \to 1}{\mathsf{\text{lim}}}f\left(\rho \right)=\infty $ (Figure 6). Eq (14) shows that $\rho $ has a significant effect on the period MPS. In order to reduce the period MPS, i.e., to make the period robust to parameter fluctuations, $\rho $ should not be close to one. Eq (14) also shows the feedback loop length $n$ is important. The following section explains how the period MPS depends on the feedback loop length.
Effect of changes in feedback loop length on period and its MPS
Assuming ${\rho}_{i}=0.5$ (${K}_{i}={\beta}_{i}/2{\alpha}_{i}$), Eq (9) and Eq (12) reduce to
respectively. Therefore, the value of period MPS is constrained by
The minimal value of the period MPS decreases as n increases. Here, let ${\Psi}_{\tau}$be the sum of all the singleparameter sensitivities of the period. Then, we get an interesting relationship among the singleparameter sensitivities (see Eqs (61)(63) in Methods):
Note that ${S}_{{\beta}_{i}}^{\tau}={S}_{{K}_{i}}^{\tau}=0$ and ${S}_{{\alpha}_{i}}^{\tau}<0$ when ${\rho}_{i}=0.5$. Eq (19) indicates that all the degradation reactions share the influence on the period. The MPS given by Eq (17) is minimized when all the degradation reactions equally share the influence (${S}_{{\alpha}_{i}}^{\tau}=1/n$), which is achieved by equating all ${\alpha}_{i}$. Assuming ${\alpha}_{i}=\alpha $, Eq (16) and Eq (17) further reduce to
respectively. Eq (20) indicates that the period increases with an increase in the feedback loop length $n$ and/or with a decrease in the decay rate constant $\alpha $. Eq (21) indicates that changes in $\alpha $ do not alter the MPS. On the other hand, the MPS decreases as n increases.
In the negative feedback oscillators, the period is given by the sum of the time delays generated by the reactions belonging to the feedback loop (Figure 7). As the number of the reactions increases (thus, n increases), the time delays can be distributed to more reactions. This 'distributedtime delay' mechanism provides the oscillation period with robustness to parameter fluctuations, decreasing the MPS.
In the previous paper [8], we proved that the merge reactions with addition logic contribute to keeping a steadystate component concentration constant against fluctuations in kinetic parameters. The concentration MPS decreases as the number of merging influx reactions increases (see Appendix A4 in [8]). The period MPS given by Eq (17) resembles the concentration MPS in the merge reactions. Although the oscillation period and steadystate concentration are different properties, our studies suggest that the underlying mechanisms that provide robustness are common to both the cases. The oscillation period is given by the sum of time delays, and the steadystate concentration is determined by the sum of influxes.
Conclusions
We presented the analytical solutions of period, amplitude, and their associated MPSs for a general model of negative feedback oscillators. We validated the analytical solutions by comparing them with numerical solutions. Next, using the analytical solutions, we investigated how changes in kinetic parameters affect the period and its associated MPS. ${\rho}_{i}$, the ratio of threshold value to the amplitude, was found to be an important determinant of the period MPS. When ${\rho}_{i}$ is close to one (${K}_{i}\approx {\beta}_{i}/{\alpha}_{i}$ ), MPS is very large, indicating that the period is very sensitive to fluctuations in kinetic parameters. Finally, we gave the first mathematical proof that long negative feedback loops make oscillations robust to parameter fluctuations by the distributed timedelay mechanism. Since the analytical solutions were derived for a general model of negative feedback loops that are common structures in biological oscillators, the results presented in this paper are expected to be true for many of biological oscillators.
In circadian rhythm networks, clock proteins have multiple phosphorylation sites. It has been reported that the phosphorylation is responsible for sleep disorders [29–31] and temperature compensation [32]. In the previous study, we carried out numerical simulations to predict that the multiple phosphorylations contribute to enhanced robustness of circadian rhythms [8]. The theoretical study in this paper strongly supports our previous result. The multiple phosphorylations form a long negative feedback loop, providing the distributed timedelay mechanism. It is reasonable for organisms to devise the multiple phosphorylation sites in order to make circadian cycles robust.
To date synthetic biologists have focused on whether their artificial gene regulatory circuits oscillate [22, 28, 33–35]. Their next challenge would be to design oscillator circuits with desired and accurate period. Our study demonstrates key factors for designing such a highperformance clock. For instance, ${\rho}_{i}$ ($={K}_{i}{\alpha}_{i}/{\beta}_{i}$) should not be close to one to make the period insensitive to parameter fluctuations (Eqs (14)(15)). When ${\rho}_{i}=0.5959$, the period becomes robust to parameter fluctuations. The robustness can be further increased by lengthening a feedback loop. To design an oscillator with long period, one can employ a long negative feedback loop (large n) and/or slow decay rates (small ${\alpha}_{i}$) as suggested by Eq (9). Our analytical study recommends increasing the feedback loop length instead of decreasing decay rates, achieving a robust oscillator (Eqs (17)(18) and Eq (21)). In summary, the analytical solutions presented in this paper greatly contribute to designing robust gene oscillators with desired period.
Methods
Derivation of analytical solutions for period and amplitude
In this section, we explain how to derive the analytical solutions for the period and amplitude. For simplicity, we assume $n=3$ and Eq (1) becomes
An example of the dynamics is shown in Figure 2. First, we divide a cycle of oscillations into time intervals (${I}_{ij}$ in Figure 2). In ${I}_{i1}$, ${x}_{i}$ increases from its trough to ${K}_{i}$. In ${I}_{i2}$, ${x}_{i}$ decreases from its peak to ${K}_{i}$. The time taken for ${I}_{ij}$is denoted by ${\tau}_{ij}$. Next, we derive analytical expressions for each time interval, and the peaks and troughs of oscillation. Finally, we connect all the time intervals to obtain period. The amplitude is obtained by subtracting the trough from the peak. We set $t=0$ to the starting time of ${I}_{11}$.
In ${I}_{11}$, ${I}_{21}$ and ${I}_{31}$, the first molecular species is produced ($\theta =1$) and thus the differential equation for ${x}_{1}$ in Eqs (22) becomes
Integrating Eq (23), we obtain
Similarly, in ${I}_{21}$, ${I}_{31}$ and ${I}_{12}$, ${x}_{2}$ follows
Since $t=0$ is the starting time of ${I}_{11}$, the subtraction of ${\tau}_{11}$ appears in Eq (25). The same shall apply to Eq (26) and Eqs (28)(30). In ${I}_{31}$, ${I}_{12}$ and ${I}_{22}$, ${x}_{3}$ follows
In ${I}_{12}$, ${I}_{22}$ and ${I}_{32}$, the first molecular species is not produced ($\theta =0$) and thus the differential equation for ${x}_{1}$ in Eqs (22) becomes
Integrating Eq (27), we obtain
Similarly, in ${I}_{22}$, ${I}_{32}$ and ${I}_{11}$, ${x}_{2}$ follows
In ${I}_{32}$, ${I}_{11}$ and ${I}_{21}$, ${x}_{3}$ follows
${C}_{ij}$($i\in \left\{1,2,3\right\}$, $j\in \left\{1,2\right\}$) are the integral constants. Here, we would like to calculate ${\tau}_{ij}$. Since ${x}_{i}$ increases and reaches to ${K}_{i}$ at $t={\displaystyle \sum _{j=1}^{i}}{\tau}_{j1}$ (Figure 2), we obtain
By inserting Eqs (24)(26), Eq (31) becomes
By solving Eq (32) for ${\tau}_{i1}$, we obtain
where
Since ${x}_{i}$ decreases and reaches to ${K}_{i}$ at $t={\displaystyle \sum _{j=1}^{3}}{\tau}_{j1}+{\displaystyle \sum _{j=1}^{i}}{\tau}_{j2}$ (Figure 2), we obtain
By inserting Eqs (28)(30), Eq (35) becomes
By solving Eq (36) for ${\tau}_{i2}$, we obtain
where
By inserting Eqs (33) and (37) into Eqs (24)(26), we get the peaks of ${x}_{1}$, ${x}_{2}$ and ${x}_{3}$:
By inserting Eqs (33) and (37) into Eqs (28)(30), we get the troughs of ${x}_{1}$, ${x}_{2}$ and ${x}_{3}$:
By summing up ${\tau}_{ij}$, we obtain the period:
The oscillation amplitude is calculated by subtracting the trough from the peak. From Eqs (24) and (39), the amplitude of ${x}_{1}$ is
Similarly, from Eqs (25) and (40), the amplitude of ${x}_{2}$ is
From Eqs (26) and (41), the amplitude of ${x}_{3}$ is
The general expressions for the period and amplitude are shown as Eqs (3) and (4), respectively. To calculate $\tau $ and ${\epsilon}_{i}$ ($i\in \left\{1,2,3\right\}$), we have to determine ${C}_{ij}$ ($i\in \left\{1,2,3\right\},j\in \left\{1,2\right\}$). ${C}_{ij}$ are determined by solving the following six equations. ${x}_{1}$ at the end of ${I}_{32}$ equals to that at the starting point of ${I}_{11}$. By equating Eq (42) to ${x}_{1}\left(0\right)$ of Eq (24),
${x}_{2}$ at the end of ${I}_{11}$ equals to that at the starting point of ${I}_{21}$. By equating Eq (43) to ${x}_{2}\left({\tau}_{11}\right)$ of Eq (25),
${x}_{3}$ at the end of ${I}_{21}$ equals to that at the starting point of ${I}_{31}$. By equating Eq (44) to ${x}_{3}\left({\tau}_{11}+{\tau}_{21}\right)$ of Eq (26),
${x}_{1}$ at the end of ${I}_{31}$ equals to that at the starting point of ${I}_{12}$. By equating Eq (39) to ${x}_{1}\left({\displaystyle \sum _{i=1}^{3}}{\tau}_{i1}\right)$ of Eq (28),
${x}_{2}$ at the end of ${I}_{12}$ equals to that at the starting point of ${I}_{22}$ . By equating Eq (40) to ${x}_{2}\left({\displaystyle \sum _{i=1}^{3}}{\tau}_{i1}+{\tau}_{12}\right)$ of Eq (29),
${x}_{3}$ at the end of ${I}_{22}$ equals to that at the starting point of ${I}_{32}$. By equating Eq (41) to ${x}_{3}\left({\displaystyle \sum _{i=1}^{3}}{\tau}_{i1}+{\tau}_{12}+{\tau}_{22}\right)$ of Eq (30),
The general forms of Eqs (49)(51) and Eqs (52)(54) are expressed as Eq (7) and Eq (8), respectively. Since there are six integral constants (${C}_{ij}$) and six equations (Eqs (49)(54)), the integral constants can be calculated. For this purpose, we used fsolve of MATLAB in this study.
Approximation of analytical solutions for period and amplitude
Assuming $\theta =1$ for Eqs (1), ${x}_{i}$ moves toward its steadystate value ${\beta}_{i}/{\alpha}_{i}$, which is given by setting $d{x}_{i}/dt=0$. Similarly, assuming $\theta =0$, ${x}_{i}$ moves toward zero. As shown in Figure 2 and 5, the peak and trough of ${x}_{i}$ are close to ${\beta}_{i}/{\alpha}_{i}$ and zero, respectively. Therefore, we assume that the peak and trough are ${\beta}_{i}/{\alpha}_{i}$ and zero, respectively, i.e., ${x}_{i}$ increases (decreases) to reach its steady state before it begins to decreases (increases). This assumption allows us to symbolically derive the period and amplitude as shown below. The assumption is validated in "Results and discussion."
The assumption that the trough is zero gives
Using Eq (55) and Eqs (24)(26), the integral constants ${C}_{i1}$ are
The assumption that the peak is ${\beta}_{i}/{\alpha}_{i}$ gives
Using Eqs (57) and Eqs (28)(30), the integral constants ${C}_{i2}$ are
By inserting Eqs (56) and (58) into Eqs (34) and (38), respectively, the period $\tau $ described by Eq (3) (or Eq (45)) becomes
where ${\rho}_{i}={K}_{i}{\alpha}_{i}/{\beta}_{i}$. From the assumption that the peak and trough are ${\beta}_{i}/{\alpha}_{i}$ and zero, respectively, the amplitude of ${x}_{i}$ is
The relative change of a system property in response to a relative change in a parameter is called singleparameter sensitivity (see the following section). The singleparameter sensitivities for period are
MPS is given by summing the squared singleparameter sensitivities and represents system's fragility to small, random, and simultaneous fluctuations in all kinetic parameters (see the following section). Period MPS is given by
The singleparameter sensitivities of amplitude are
Amplitude MPS is given by
Multiparameter sensitivity (MPS)
Generally a dynamic model for biochemical networks is described by ordinary differential equations:
where t is time, x is the vector whose elements are the variables for molecular concentrations, $p=\left({p}_{1},{p}_{2},\dots ,{p}_{m}\right)$ is the kinetic parameter vector, and m is the number of kinetic parameters. In this study, $p=\left({\alpha}_{1},\dots ,{\alpha}_{n},{\beta}_{1},\dots ,{\beta}_{n},{K}_{1},\dots ,{K}_{n}\right)$, where n is the number of molecular species. Let q(p) be a given system's output (oscillation period and amplitude in this study) which depends on the kinetic parameter vector. The singleparameter sensitivity of the output with respect to a change in the i th parameter is given by
Singleparameter sensitivities are used to identify parameters influential on the output.
Assuming that the relative change in the output is the linear combination of a change in each parameter, multiparameter sensitivity (MPS) [8, 36–38] is given by
Although MPS is expressed as the sum of the squared singleparameter sensitivities, it represents how fragile the system's output is when small, random, and simultaneous fluctuations are provided to all kinetic parameters [8]. MPS can be used as an indicator of biochemical system's robustness in a cell, where kinetic parameters constantly fluctuate. MPS is mathematically equal to the normalized variance given by the MonteCarlo method [9–11], where all kinetic parameters are simultaneously and randomly deviated from the nominal values. For more details, see our previous work [8].
Numerical computation of multiparameter sensitivity (MPS)
Generally, it is hard to derive the analytical solution for MPS when the given biochemical models are complicated. As a practical solution, the MPS is numerically computed by providing a small perturbation to kinetic parameters (${\alpha}_{i}$, ${\beta}_{i}$ and ${K}_{i}$ in this study) [8]. Eq (70) can be rewritten as
where $p=\left({p}_{1},\dots ,{p}_{i},\dots ,{p}_{m}\right)$ is the standard parameter vector, ${p}^{\prime}=\left({p}_{1},\dots ,{p}_{i}+\Delta {p}_{i},\dots ,{p}_{m}\right)$ is the perturbed parameter vector, and m is the number of kinetic parameters. In the numerical computation of singleparameter sensitivity, $\Delta {p}_{i}$ is not infinitesimal but a small value ($\Delta {p}_{i}=1{0}^{3}{p}_{i}$ in this study). The numerical computation shown here was used to obtain the numerical integration solutions and semianalytical solutions (see Table 1). We validated the numerical computation in the previous work [8].
Abbreviations
 MPS:

Multiparameter sensitivity.
References
 1.
BellPedersen D, Cassone VM, Earnest DJ, Golden SS, Hardin PE, Thomas TL, Zoran MJ: Circadian rhythms from multiple oscillators: lessons from diverse organisms. Nat Rev Genet. 2005, 6 (7): 544556. 10.1038/nrg1633.
 2.
Dunlap JC: Molecular bases for circadian clocks. Cell. 1999, 96 (2): 271290. 10.1016/S00928674(00)805668.
 3.
Panda S, Hogenesch JB, Kay SA: Circadian rhythms from flies to human. Nature. 2002, 417 (6886): 329335. 10.1038/417329a.
 4.
Tyson JJ, Novak B: Regulation of the eukaryotic cell cycle: molecular antagonism, hysteresis, and irreversible transitions. J Theor Biol. 2001, 210 (2): 249263. 10.1006/jtbi.2001.2293.
 5.
Csete ME, Doyle JC: Reverse engineering of biological complexity. Science. 2002, 295 (5560): 16641669. 10.1126/science.1069981.
 6.
ElSamad H, Kurata H, Doyle JC, Gross CA, Khammash M: Surviving heat shock: control strategies for robustness and performance. Proc Natl Acad Sci USA. 2005, 102 (8): 27362741. 10.1073/pnas.0403510102.
 7.
Kurata H, ElSamad H, Iwasaki R, Ohtake H, Doyle JC, Grigorova I, Gross CA, Khammash M: Modulebased analysis of robustness tradeoffs in the heat shock response system. PLoS Comput Biol. 2006, 2 (7): e5910.1371/journal.pcbi.0020059.
 8.
Maeda K, Kurata H: Quasimultiparameter sensitivity measure for robustness analysis of complex biochemical networks. J Theor Biol. 2011, 272 (1): 174186. 10.1016/j.jtbi.2010.12.012.
 9.
Stelling J, Gilles ED, Doyle FJ: Robustness properties of circadian clock architectures. Proc Natl Acad Sci USA. 2004, 101 (36): 1321013215. 10.1073/pnas.0401463101.
 10.
Wagner A: Circuit topology and the evolution of robustness in twogene circadian oscillators. Proc Natl Acad Sci USA. 2005, 102 (33): 1177511780. 10.1073/pnas.0501094102.
 11.
Kurata H, Tanaka T, Ohnishi F: Mathematical identification of critical reactions in the interlocked feedback model. PLoS ONE. 2007, 2 (10): e110310.1371/journal.pone.0001103.
 12.
Leloup JC, Goldbeter A: A model for circadian rhythms in Drosophila incorporating the formation of a complex between the PER and TIM proteins. J Biol Rhythms. 1998, 13 (1): 7087. 10.1177/074873098128999934.
 13.
Maeda K, Kurata H: A symmetric dual feedback system provides a robust and entrainable oscillator. PLoS ONE. 2012, 7 (2): e3048910.1371/journal.pone.0030489.
 14.
Olivier BG, Snoep JL: Webbased kinetic modelling using JWS Online. Bioinformatics. 2004, 20 (13): 21432144. 10.1093/bioinformatics/bth200.
 15.
Le Novere N, Bornstein B, Broicher A, Courtot M, Donizelli M, Dharuri H, Li L, Sauro H, Schilstra M, Shapiro B: BioModels Database: a free, centralized database of curated, published, quantitative kinetic models of biochemical and cellular systems. Nucleic Acids Res. 2006, 34 (Database): D689691.
 16.
Li C, Donizelli M, Rodriguez N, Dharuri H, Endler L, Chelliah V, Li L, He E, Henry A, Stefan MI: BioModels Database: An enhanced, curated and annotated resource for published quantitative kinetic models. BMC Syst Biol. 2010, 4: 9210.1186/17520509492.
 17.
Kurata H, Maeda K, Onaka T, Takata T: BioFNet: biological functional network database for analysis and synthesis of biological systems. Brief Bioinform. 2014, 15 (5): 699709. 10.1093/bib/bbt048.
 18.
Nguyen LK, Kulasiri D: On the functional diversity of dynamical behaviour in genetic and metabolic feedback systems. BMC Syst Biol. 2009, 3: 5110.1186/17520509351.
 19.
Nguyen LK: Regulation of oscillation dynamics in biochemical systems with dual negative feedback loops. J R Soc Interface. 2012, 9 (73): 19982010. 10.1098/rsif.2012.0028.
 20.
Kurosawa G, Mochizuki A, Iwasa Y: Comparative study of circadian clock models, in search of processes promoting oscillation. J Theor Biol. 2002, 216 (2): 193208. 10.1006/jtbi.2002.2546.
 21.
Kut C, Golkhou V, Bader JS: Analytical approximations for the amplitude and period of a relaxation oscillator. BMC Syst Biol. 2009, 3: 610.1186/1752050936.
 22.
Elowitz MB, Leibler S: A synthetic oscillatory network of transcriptional regulators. Nature. 2000, 403 (6767): 335338. 10.1038/35002125.
 23.
Barkai N, Leibler S: Circadian clocks limited by noise. Nature. 2000, 403 (6767): 267268.
 24.
Vilar JM, Kueh HY, Barkai N, Leibler S: Mechanisms of noiseresistance in genetic oscillators. Proc Natl Acad Sci USA. 2002, 99 (9): 59885992. 10.1073/pnas.092133899.
 25.
Novak B, Tyson JJ: Design principles of biochemical oscillators. Nat Rev Mol Cell Biol. 2008, 9 (12): 981991. 10.1038/nrm2530.
 26.
Tyson JJ, Novak B: Functional motifs in biochemical reaction networks. Annu Rev Phys Chem. 2010, 61: 219240. 10.1146/annurev.physchem.012809.103457.
 27.
Alon U: An Introduction to Systems Biology: Design Principles of Biological Circuits. 2006, Chapman and Hall/CRC
 28.
Atkinson MR, Savageau MA, Myers JT, Ninfa AJ: Development of genetic circuitry exhibiting toggle switch or oscillatory behavior in Escherichia coli. Cell. 2003, 113 (5): 597607. 10.1016/S00928674(03)003465.
 29.
Leloup JC: Circadian clocks and phosphorylation: Insights from computational modeling. Cent Eur J Biol. 2009, 4 (3): 290303. 10.2478/s1153500900251.
 30.
Leloup JC, Goldbeter A: Toward a detailed computational model for the mammalian circadian clock. Proc Natl Acad Sci USA. 2003, 100 (12): 70517056. 10.1073/pnas.1132112100.
 31.
Leloup JC, Goldbeter A: Modeling the circadian clock: from molecular mechanism to physiological disorders. Bioessays. 2008, 30 (6): 590600. 10.1002/bies.20762.
 32.
Jolley CC, Ode KL, Ueda HR: A Design Principle for a Posttranslational Biochemical Oscillator. Cell reports. 2012, 2 (4): 938950. 10.1016/j.celrep.2012.09.006.
 33.
Tigges M, MarquezLago TT, Stelling J, Fussenegger M: A tunable synthetic mammalian oscillator. Nature. 2009, 457 (7227): 309312. 10.1038/nature07616.
 34.
Stricker J, Cookson S, Bennett MR, Mather WH, Tsimring LS, Hasty J: A fast, robust and tunable synthetic gene oscillator. Nature. 2008, 456 (7221): 516519. 10.1038/nature07389.
 35.
Fung E, Wong WW, Suen JK, Bulter T, Lee SG, Liao JC: A synthetic genemetabolic oscillator. Nature. 2005, 435 (7038): 118122. 10.1038/nature03508.
 36.
Rosenblum AL, Ghausi MS: Multiparameter Sensitivity in Active RC Networks. IEEE Transactions on Circuit Theory. 1971, 18 (6): 592599.
 37.
Schoeffler J: The Synthesis of Minimum Sensitivity Networks. IEEE Transactions on Circuit Theory. 1964, 11 (2): 271276.
 38.
Leeds JVJ, Ugron GI: Simplified Multiple Parameter Sensitivity Calculation and Continuously Equivalent Networks. IEEE Transactions on Circuit Theory. 1967, 14 (2): 188191.
Acknowledgements
This work was supported by GrantinAid for Scientific Research (B) (25280107) and Research Activity Startup (24800050) from the Ministry of Education, Culture, Sports, Science and Technology of Japan.
Declarations
Publication charges for this work were funded by GrantinAid for Scientific Research (B) (25280107) from the Ministry of Education, Culture, Sports, Science and Technology of Japan.
This article has been published as part of BMC systems Biology Volume 8 Supplement 5, 2014: Proceedings of the 25th International Conference on Genome Informatics (GIW/ISCBAsia): Systems Biology. The full contents of the supplement are available online at http://www.biomedcentral.com/bmcsystbiol/supplements/8/S5.
Author information
Affiliations
Corresponding author
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors' contributions
KM and HK participated in the design of the study. KM performed the mathematical analysis. KM and HK analyzed the data and wrote the paper. Both authors read and approved the final manuscript.
Electronic supplementary material
Rights and permissions
About this article
Cite this article
Maeda, K., Kurata, H. Analytical study of robustness of a negative feedback oscillator by multiparameter sensitivity. BMC Syst Biol 8, S1 (2014). https://doi.org/10.1186/175205098S5S1
Published:
Keywords
 Negative Feedback Oscillator
 Multiparameter Sensitivity (MPS)
 Biological Oscillators
 Singleparameter Sensitivity
 Numerical Integration Solution