PMC:4305980 JSONTXT 3 Projects

Analytical study of robustness of a negative feedback oscillator by multiparameter sensitivity 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 time-delay 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 ever-fluctuating 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 single-parameter 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 PER-TIM 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 Elowitz-Leibler repressilator [22] and Barkai-Leibler 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 time-delay 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: Figure 1 Schematic diagram of the negative feedback oscillator model. xi is the ith molecular species, αi is the decay rate constant, βi is the production rate constant, and Ki is the threshold for turning on/off the production of the target molecular species (i∈{1,2,…,n}). n is the number of molecular species. xi activates xi+1 (i∈{1,2,…,n-1}). xn represses x1. All molecular species are subject to degradation. (1) d x 1 d t = β 1 ⋅ θ ( K n , x n ) - α 1 x 1 d x 2 d t = β 2 ⋅ θ ( x 1 , K 1 ) - α 2 x 2 ⋮ d x n d t = β n ⋅ θ ( x n - 1 , K n - 1 ) - α n x n where xi is the concentration of the ith molecular species (e.g., transcription factor, specifically modified form of protein, metabolite, etc.), αi is the decay rate constant, βi is the production rate constant, and Ki is the threshold for turning on/off the production of the target molecular species (i∈{1,2,…,n}). αi, βi and Ki take positive values. n is the number of molecular species or indicates feedback loop length.  θ is the unit step function given by (2) θ ( a , b ) = 0 ( a < b ) 1 ( a ≥ b ) 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. Ki corresponds to the concentration of a regulator molecule at which the production rate of a target molecule reaches a half maximum in Hill function. xi takes a positive value in the range of (0, βi/αi). The negative feedback oscillator model can produce oscillations when both n≥3 and Ki<βi/αi (i∈{1,…,n}) 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. Figure 2 Example of dynamics of the negative feedback oscillator model. n=3, αi=1, βi=1 and Ki=0.5 (i∈{1,2,3}). Iij are the intervals used to derive the analytical solutions for the period and amplitude (see Text). 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 xi are clock genes or specifically modified forms of clock proteins. Our model can also be found in the field of synthetic biology. Elowitz-Leibler 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 (Iij 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 and the ith species amplitude εi for the oscillatory model of Eq (1) are given by (3) τ = - ∑ i = 1 n ln ( γ i δ i ) α i (4) ε i = C 11 ∏ j = 1 n γ j α 1 / α j - 1   ( i = 1 ) C i 1 ∏ j = i n γ j α i / α j ∏ j = 1 i - 1 δ j α i / α j - 1   ( i > 1 ) where γi and δi are (5) γ i = K i - β i / α i C i 1 (6) δ i = K i C i 2 Cij (i∈{1,…,n}, j∈{1,2}) are the integral constants: (7) C i 1 = K 1 ∏ j = 2 n δ j α 1 / α j - β 1 α 1   ( i = 1 ) K i ∏ j = 1 i - 1 γ j α i / α j ∏ j = i + 1 n δ j α i / α j - β i α i   ( 1 < i < n ) K n ∏ j = 1 n - 1 γ j α n / α j - β n α n   ( i = n ) (8) C i 2 = C 11 ∏ j = 1 n γ j α 1 / α j + β 1 α 1   ( i = 1 ) C i 1 ∏ j = i n γ j α i / α j ∏ j = 1 i - 1 δ j α i / α j + β i α i   ( i > 1 ) Although Cij 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 xi are βi/αi and zero, respectively, the integral constants can be determined without any numerical computations, and the analytical solutions for the period and the ith species amplitude εi are given by (for details, see Methods) (9) τ = - ∑ i = 1 n ln [ ρ i ( 1 - ρ i ) ] α i (10) ε i = β i α i where ρi is the ratio of the threshold Ki to the amplitude βi/αi: (11) ρ i = K i α i β i Multiparameter sensitivity (MPS) represents how fragile system's properties are to fluctuations in kinetic parameters (for details, see Methods). The period MPS Φτ and amplitude MPS Φεi are given by (12) Φ τ = ∑ i = 1 n ln [ ρ i ( 1 - ρ i ) ] α i - 1 - 2 ρ i α i ( 1 - ρ i ) 2 + 2 1 - 2 ρ i α i ( 1 - ρ i ) 2 ∑ i = 1 n ln [ ρ i ( 1 - ρ i ) ] α i 2 (13) Φ ε i = 2 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 βi/αi peak and zero trough). In order to determine Cij, 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 αi, βi and Ki (see Methods). The computed period, amplitude and MPSs are referred to as the semi-analytical solutions because numerical methods were partially used. The semi-analytical 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 semi-analytical and numerical integration solutions. We assigned uniform random values over (0, 1) to αi and βi, and those over (0, βi/αi) to Ki. The semi-analytical solutions were consistent with the numerical integration solutions (Figure 3). Thus, the analytical solutions given by Eqs (3)-(4) were validated. Table 1 Summary of how to obtain numerical integration, semi-analytical, and full analytical solutions Figure 3 Comparison of the semi-analytical and numerical integration solutions . A: period, B: amplitude, C: period MPS, D: amplitude MPS. The values of αi and βi (i∈{1,2,…,n}) were randomized with a range of (0, 1), while the value of Ki (i∈{1,2,…,n}) was randomized with a range of (0, βi/αi). In addition to the correctness of the solutions, the semi-analytical approach is computationally more efficient than the numerical integration approach. When n=3, the semi-analytical 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 βi/α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 semi-analytical solutions. We assigned random values to kinetic parameters as we did for the comparison between the semi-analytical and numerical integration solutions. The results are shown in Figure 4. The full and semi-analytical solutions are consistent especially when the feedback loop is long, because the assumptions (the peak is βi/α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 semi-analytical 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 semi-analytical solutions are consistent. In summary, the analytical solutions given by Eqs (9)-(10) and Eqs (12)-(13) were validated. Figure 4 Comparison of the full analytical and semi-analytical solutions. A: period, B: amplitude, C: period MPS, D: amplitude MPS. The values of αi and βi (i∈{1,2,…,n}) were randomized with a range of (0, 1), while the value of Ki (i∈{1,2,…,n}) was randomized with a range of (0, βi/αi). Figure 5 Examples of dynamics of the negative feedback oscillator model. A: n=5, B: n=7. αi=1, βi=1 and Ki=0.5 (i∈{1,2,…,n}). 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, αi=α, βi=β, Ki=K (i∈{1,2,3}), and ρ=Kα/β. When the period is given as a function of , it reaches the minimal values at ρ=0.7228 (Figure S1AB in Additional file 1). When the period is given as a function of β or K, it reaches the minimal values at (Figure S1CDEF). Note that the negative feedback oscillator model produces oscillations only when K<β/α (ρ<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 ρ=0.5959, implying that the period MPS depends solely on ρ (when n is fixed). Assuming αi=α and ρi=ρ (i∈{1,2,…,n}), Eq (12) reduces to (14) Φ τ = f ( ρ ) n where (15) f ( ρ ) = ln [ ρ ( 1 - ρ ) ] - 1 - 2 ρ 1 - ρ 2 + 2 1 - 2 ρ 1 - ρ 2 ln [ ρ ( 1 - ρ ) ] 2 f(ρ) reaches the minimal value 0.6667 at ρ=0.5959, and limρ→1f(ρ)=∞ (Figure 6). Eq (14) shows that ρ 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, ρ 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. Figure 6 f(ρ) vs. ρ. f(ρ) is given by Eq (15) in Text. Effect of changes in feedback loop length on period and its MPS Assuming ρi=0.5 (Ki=βi/2αi), Eq (9) and Eq (12) reduce to (16) τ = ln 4 ⋅ ∑ i = 1 n α i - 1 (17) Φ τ = ∑ i = 1 n α i - 2 ∑ i = 1 n α i - 1 2 respectively. Therefore, the value of period MPS is constrained by (18) 1 n ≤ Φ τ ≤ 1 The minimal value of the period MPS decreases as n increases. Here, let Ψτbe the sum of all the single-parameter sensitivities of the period. Then, we get an interesting relationship among the single-parameter sensitivities (see Eqs (61)-(63) in Methods): (19) Ψ τ = ∑ i = 1 n S α i τ + S β i τ + S K i τ = ∑ i = 1 n S α i τ = - 1 Note that Sβiτ=SKiτ=0 and Sαiτ<0 when ρ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αiτ=-1/n), which is achieved by equating all αi. Assuming αi=α, Eq (16) and Eq (17) further reduce to (20) τ = n ln 4 α (21) Φ τ = 1 n 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 α. Eq (21) indicates that changes in α 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 'distributed-time delay' mechanism provides the oscillation period with robustness to parameter fluctuations, decreasing the MPS. Figure 7 Distributed time-delay mechanism. Oscillators with two (left) and eight (right) time delays. The right is more robust than the left because of the distributed time-delay mechanism (see Text). In the previous paper [8], we proved that the merge reactions with addition logic contribute to keeping a steady-state 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 steady-state 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 steady-state 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. ρi, the ratio of threshold value to the amplitude, was found to be an important determinant of the period MPS. When ρi is close to one (Ki≈βi/α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 time-delay 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 time-delay 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 high-performance clock. For instance, ρi (=Kiαi/βi) should not be close to one to make the period insensitive to parameter fluctuations (Eqs (14)-(15)). When ρ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 α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 (22) d x 1 d t = β 1 ⋅ θ ( K 3 , x 3 ) - α 1 x 1 d x 2 d t = β 2 ⋅ θ ( x 1 , K 1 ) - α 2 x 2 d x 3 d t = β 3 ⋅ θ ( x 2 , K 2 ) - α 3 x 3 An example of the dynamics is shown in Figure 2. First, we divide a cycle of oscillations into time intervals (Iij in Figure 2). In Ii1, xi increases from its trough to Ki. In Ii2, xi decreases from its peak to Ki. The time taken for Iijis denoted by τ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 I11. In I11, I21 and I31, the first molecular species is produced (θ=1) and thus the differential equation for x1 in Eqs (22) becomes (23) d x 1 d t = β 1 - α 1 x 1 Integrating Eq (23), we obtain (24) x 1 ( t ) = β 1 α 1 + C 11 e - α 1 t Similarly, in I21, I31 and I12, x2 follows (25) x 2 ( t ) = β 2 α 2 + C 21 e - α 2 ( t - τ 11 ) Since t=0 is the starting time of I11, the subtraction of τ11 appears in Eq (25). The same shall apply to Eq (26) and Eqs (28)-(30). In I31, I12 and I22, x3 follows (26) x 3 ( t ) = β 3 α 3 + C 31 e - α 3 ( t - τ 11 - τ 21 ) In I12, I22 and I32, the first molecular species is not produced (θ=0) and thus the differential equation for x1 in Eqs (22) becomes (27) d x 1 d t = - α 1 x 1 Integrating Eq (27), we obtain (28) x 1 ( t ) = C 12 e - α 1 ( t - τ 11 - τ 21 - τ 31 ) Similarly, in I22, I32 and I11, x2 follows (29) x 2 ( t ) = C 22 e - α 2 ( t - τ 11 - τ 21 - τ 31 - τ 12 ) In I32, I11 and I21, x3 follows (30) x 3 ( t ) = C 32 e - α 3 ( t - τ 11 - τ 21 - τ 31 - τ 12 - τ 22 ) Cij(i∈{1,2,3}, j∈{1,2}) are the integral constants. Here, we would like to calculate τij. Since xi increases and reaches to Ki at t= ∑j=1iτj1 (Figure 2), we obtain (31) x i ( ∑ j = 1 i τ j 1 ) = K i By inserting Eqs (24)-(26), Eq (31) becomes (32) β i α i + C i 1 e - α i τ i 1 = K i By solving Eq (32) for τi1, we obtain (33) τ i 1 = - ln γ i α i where (34) γ i = K i - β i / α i C i 1 Since xi decreases and reaches to Ki at t= ∑j=13τj1+ ∑j=1iτj2 (Figure 2), we obtain (35) x i ( ∑ j = 1 3 τ j 1 + ∑ j = 1 i τ j 2 ) = K i By inserting Eqs (28)-(30), Eq (35) becomes (36) C i 2 e - α i τ i 2 = K i By solving Eq (36) for τi2, we obtain (37) τ i 2 = - ln δ i α i where (38) δ i = K i C i 2 By inserting Eqs (33) and (37) into Eqs (24)-(26), we get the peaks of x1, x2 and x3: (39) x 1 ( ∑ i = 1 3 τ i 1 ) = β 1 α 1 + C 11 γ 1 γ 2 α 1 / α 2 γ 3 α 1 / α 3 (40) x 2 ( ∑ i = 1 3 τ i 1 + τ 12 ) = β 2 α 2 + C 21 γ 2 γ 3 α 2 / α 3 δ 1 α 2 / α 1 (41) x 3 ( ∑ i = 1 3 τ i 1 + τ 12 + τ 22 ) = β 3 α 3 + C 31 γ 3 δ 1 α 3 / α 1 δ 2 α 3 / α 2 By inserting Eqs (33) and (37) into Eqs (28)-(30), we get the troughs of x1, x2 and x3: (42) x 1 ( ∑ i = 1 3 τ i 1 + ∑ i = 1 3 τ i 2 ) = K 1 δ 2 α 1 / α 2 δ 3 α 1 / α 3 (43) x 2 ( ∑ i = 1 3 τ i 1 + ∑ i = 1 3 τ i 2 + τ 11 ) = K 2 δ 3 α 2 / α 3 γ 1 α 2 / α 1 (44) x 3 ( ∑ i = 1 3 τ i 1 + ∑ i = 1 3 τ i 2 + τ 11 + τ 21 ) = K 3 γ 1 α 3 / α 1 γ 2 α 3 / α 2 By summing up τij, we obtain the period: (45) τ = ∑ i = 1 3 τ i 1 + τ i 2 = - ∑ i = 1 3 ln ( γ i δ i ) α i The oscillation amplitude is calculated by subtracting the trough from the peak. From Eqs (24) and (39), the amplitude of x1 is (46) ε 1 = x 1 ( ∑ i = 1 3 τ i 1 ) - x 1 ( 0 ) = C 11 ( γ 1 γ 2 α 1 / α 2 γ 3 α 1 / α 3 - 1 ) Similarly, from Eqs (25) and (40), the amplitude of x2 is (47) ε 2 = x 2 ( ∑ i = 1 3 τ i 1 + τ 12 ) - x 2 ( τ 11 ) = C 21 ( γ 2 γ 3 α 2 / α 3 δ 1 α 2 / α 1 - 1 ) From Eqs (26) and (41), the amplitude of x3 is (48) ε 3 = x 3 ( ∑ i = 1 3 τ i 1 + τ 12 + τ 22 ) - x 3 ( τ 11 + τ 21 ) = C 31 ( γ 3 δ 1 α 3 / α 1 δ 2 α 3 / α 2 - 1 ) The general expressions for the period and amplitude are shown as Eqs (3) and (4), respectively. To calculate τ and εi (i∈{1,2,3}), we have to determine Cij (i∈{1,2,3},j∈{1,2}). Cij are determined by solving the following six equations. x1 at the end of I32 equals to that at the starting point of I11. By equating Eq (42) to x1(0) of Eq (24), (49) C 11 = K 1 δ 2 α 1 / α 2 δ 3 α 1 / α 3 - β 1 α 1 x2 at the end of I11 equals to that at the starting point of I21. By equating Eq (43) to x2(τ11) of Eq (25), (50) C 21 = K 2 δ 3 α 2 / α 3 γ 1 α 2 / α 1 - β 2 α 2 x3 at the end of I21 equals to that at the starting point of I31. By equating Eq (44) to x3(τ11+τ21) of Eq (26), (51) C 31 = K 3 γ 1 α 3 / α 1 γ 2 α 3 / α 2 - β 3 α 3 x1 at the end of I31 equals to that at the starting point of I12. By equating Eq (39) to x1(∑i=13τi1) of Eq (28), (52) C 12 = C 11 γ 1 γ 2 α 1 / α 2 γ 3 α 1 / α 3 + β 1 α 1 x2 at the end of I12 equals to that at the starting point of I22 . By equating Eq (40) to x2(∑i=13τi1+τ12) of Eq (29), (53) C 22 = C 21 γ 2 γ 3 α 2 / α 3 δ 1 α 2 / α 1 + β 2 α 2 x3 at the end of I22 equals to that at the starting point of I32. By equating Eq (41) to x3(∑i=13τi1+τ12+τ22) of Eq (30), (54) C 32 = C 31 γ 3 δ 1 α 3 / α 1 δ 2 α 3 / α 2 + β 3 α 3 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 (Cij) 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 θ=1 for Eqs (1), xi moves toward its steady-state value βi/αi, which is given by setting dxi/dt=0. Similarly, assuming θ=0, xi moves toward zero. As shown in Figure 2 and 5, the peak and trough of xi are close to βi/αi and zero, respectively. Therefore, we assume that the peak and trough are βi/αi and zero, respectively, i.e., xi 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 (55) x 1 ( 0 ) = 0 x i ( ∑ j = 1 i - 1 τ j 1 ) = 0   ( i > 1 ) Using Eq (55) and Eqs (24)-(26), the integral constants Ci1 are (56) C i 1 = - β i α i The assumption that the peak is βi/αi gives (57) x 1 ( ∑ j = 1 n τ j 1 ) = β i α i x i ( ∑ j = 1 n τ j 1 + ∑ j = 1 i - 1 τ j 2 ) = β i α i   ( i > 1 ) Using Eqs (57) and Eqs (28)-(30), the integral constants Ci2 are (58) C i 2 = β i α i By inserting Eqs (56) and (58) into Eqs (34) and (38), respectively, the period τ described by Eq (3) (or Eq (45)) becomes (59) τ = - ∑ i = 1 n ln [ ρ i ( 1 - ρ i ) ] α i where ρi=Kiαi/βi. From the assumption that the peak and trough are βi/αi and zero, respectively, the amplitude of xi is (60) ε i = β i α i The relative change of a system property in response to a relative change in a parameter is called single-parameter sensitivity (see the following section). The single-parameter sensitivities for period are (61) S α i τ = α i τ ∂ τ ∂ α i = τ - 1 α i - 1 ln ρ i ( 1 - ρ i ) - 1 - 2 ρ i 1 - ρ i (62) S β i τ = β i τ ∂ τ ∂ β i = τ - 1 α i - 1 1 - 2 ρ i 1 - ρ i (63) S K i τ = K i τ ∂ τ ∂ K i = - τ - 1 α i - 1 1 - 2 ρ i 1 - ρ i MPS is given by summing the squared single-parameter 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 (64) Φ τ = ∑ i = 1 n S α i τ 2 + S β i τ 2 + S K i τ 2 = ∑ i = 1 n ln [ ρ i ( 1 - ρ i ) ] α i - 1 - 2 ρ i α i ( 1 - ρ i ) 2 + 2 1 - 2 ρ i α i ( 1 - ρ i ) 2 ∑ i = 1 n ln [ ρ i ( 1 - ρ i ) ] α i 2 The single-parameter sensitivities of amplitude are (65) S α j ε i = α j ε i ∂ ε i ∂ α j = - 1 ( i = j ) 0 ( i ≠ j ) (66) S β j ε i = β j ε i ∂ ε i ∂ β j = 1 ( i = j ) 0 ( i ≠ j ) (67) S K j ε i = K j ε i ∂ ε i ∂ K j = 0 Amplitude MPS is given by (68) Φ ε i = ∑ j = 1 n S α j ε i 2 + S β j ε i 2 + S K j ε i 2 = 2 Multiparameter sensitivity (MPS) Generally a dynamic model for biochemical networks is described by ordinary differential equations: (69) x · = F ( t , x , p ) where t is time, x is the vector whose elements are the variables for molecular concentrations, p=(p1,p2,…,pm) is the kinetic parameter vector, and m is the number of kinetic parameters. In this study, p=(α1,…,αn,β1,…,βn,K1,…,Kn), 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 single-parameter sensitivity of the output with respect to a change in the ith parameter is given by (70) S p i q = p i q ∂ q ∂ p i = ∂ ln q ∂ ln p i Single-parameter 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 (71) Φ q = ∑ i = 1 m S p i q 2 Although MPS is expressed as the sum of the squared single-parameter 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 Monte-Carlo 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 (αi, βi and Ki in this study) [8]. Eq (70) can be rewritten as (72) S p i q = p i q ( p ) lim Δ p i → 0 q ( p ′ ) - q ( p ) Δ p i = lim Δ p i → 0 ln q ( p ′ ) - ln q ( p ) ln ( p i + Δ p i ) - ln p i where p=(p1,…,pi,…,pm) is the standard parameter vector, p′=(p1,…,pi+Δpi,…,pm) is the perturbed parameter vector, and m is the number of kinetic parameters. In the numerical computation of single-parameter sensitivity, Δpi is not infinitesimal but a small value (Δpi=10-3pi in this study). The numerical computation shown here was used to obtain the numerical integration solutions and semi-analytical solutions (see Table 1). We validated the numerical computation in the previous work [8]. Abbreviation MPS: Multiparameter sensitivity. 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. Supplementary Material Additional file 1 Supplementary Figures.pdf. Figure S1 and S2 are provided. Click here for file