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. $x_i$ is the i th molecular species, ${\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\}$). n is the number of molecular species. $x_i$ activates $x_{i+1}$ ($i\in \left\{1,2,\dots ,n-1\right\}$). $x_n$ represses $x_1$. All molecular species are subject to degradation.

$$\begin{array}{c}\frac{d{x}_1}{dt}={\beta }_1\cdot \theta \left(K_n,x_n\right)-{\alpha }_1{x}_1\\ \frac{d{x}_2}{dt}={\beta }_2\cdot \theta \left(x_1,K_1\right)-{\alpha }_2{x}_2\\ \begin{array}{}\end{array}⋮\\ \frac{d{x}_n}{dt}={\beta }_n\cdot \theta \left(x_{n-1},K_{n-1}\right)-{\alpha }_n{x}_n\end{array}$$

(1)

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

$$\theta \left(a,b\right)=\left\{\begin{array}{c}\hfill \begin{array}{cc}\hfill 0\hfill & \hfill \left(a<b\right)\hfill \end{array}\hfill \\ \hfill \begin{array}{cc}\hfill 1\hfill & \hfill \left(a\ge b\right)\hfill \end{array}\hfill \end{array}\right.$$

(2)

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.

Figure 2

Example of dynamics of the negative feedback oscillator model. $n=3$, ${\alpha }_i=1$, ${\beta }_i=1$ and $K_i=0.5$ ($i\in \left\{1,2,3\right\}$). $I_{ij}$ 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 $x_i$ 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 ($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

$$\tau =-{\displaystyle \sum _{i=1}^n}\frac{\mathsf{\text{ln}}\left({\gamma }_i{\delta }_i\right)}{{\alpha }_i}$$

(3)

$${\epsilon }_i=\left\{\begin{array}{c}\hfill C_{11}\left({\displaystyle \prod _{j=1}^n}{{\gamma }_j}^{{\alpha }_1/{\alpha }_j}-1\right)\left(i=1\right)\hfill \\ \hfill C_{i1}\left({\displaystyle \prod _{j=i}^n}{{\gamma }_j}^{{\alpha }_i/{\alpha }_j}{\displaystyle \prod _{j=1}^{i-1}}{{\delta }_j}^{{\alpha }_i/{\alpha }_j}-1\right)\left(i>1\right)\hfill \end{array}\right.$$

(4)

where ${\gamma }_i$ and ${\delta }_i$ are

$${\gamma }_i=\frac{K_i-{\beta }_i/{\alpha }_i}{C_{i1}}$$

(5)

$${\delta }_i=\frac{K_i}{C_{i2}}$$

(6)

$C_{ij}$ ($i\in \left\{1,\dots ,n\right\}$, $j\in \left\{1,2\right\}$) are the integral constants:

$$C_{i1}=\left\{\begin{array}{c}\hfill K_1{\displaystyle \prod _{j=2}^n}{{\delta }_j}^{{\alpha }_1/{\alpha }_j}-\frac{{\beta }_1}{{\alpha }_1}\left(i=1\right)\hfill \\ \hfill K_i{\displaystyle \prod _{j=1}^{i-1}}{{\gamma }_j}^{{\alpha }_i/{\alpha }_j}{\displaystyle \prod _{j=i+1}^n}{{\delta }_j}^{{\alpha }_i/{\alpha }_j}-\frac{{\beta }_i}{{\alpha }_i}\left(1<i<n\right)\hfill \\ \hfill K_n{\displaystyle \prod _{j=1}^{n-1}}{{\gamma }_j}^{{\alpha }_n/{\alpha }_j}-\frac{{\beta }_n}{{\alpha }_n}\left(i=n\right)\hfill \end{array}\right.$$

(7)

$$C_{i2}=\left\{\begin{array}{c}\hfill C_{11}{\displaystyle \prod _{j=1}^n}{{\gamma }_j}^{{\alpha }_1/{\alpha }_j}+\frac{{\beta }_1}{{\alpha }_1}\left(i=1\right)\hfill \\ \hfill C_{i1}{\displaystyle \prod _{j=i}^n}{{\gamma }_j}^{{\alpha }_i/{\alpha }_j}{\displaystyle \prod _{j=1}^{i-1}}{{\delta }_j}^{{\alpha }_i/{\alpha }_j}+\frac{{\beta }_i}{{\alpha }_i}\left(i>1\right)\hfill \end{array}\right.$$

(8)

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)

$$\tau =-{\displaystyle \sum _{i=1}^n}\frac{\mathsf{\text{ln}}\left[{\rho }_i\left(1-{\rho }_i\right)\right]}{{\alpha }_i}$$

(9)

$${\epsilon }_i=\frac{{\beta }_i}{{\alpha }_i}$$

(10)

where ${\rho }_i$ is the ratio of the threshold $K_i$ to the amplitude ${\beta }_i/{\alpha }_i$:

$${\rho }_i=\frac{K_i{\alpha }_i}{{\beta }_i}$$

(11)

Multiparameter sensitivity (MPS) represents how fragile system's properties are to fluctuations in kinetic parameters (for details, see Methods). The period MPS ${\text{Φ}}_{\tau }$ and amplitude MPS ${\text{Φ}}_{{\epsilon }_i}$ are given by

$${\text{Φ}}_{\tau }=\frac{{\displaystyle \sum _{i=1}^n}\left[{\left(\frac{\mathsf{\text{ln}}\left[{\rho }_i\left(1-{\rho }_i\right)\right]}{{\alpha }_i}-\frac{1-2{\rho }_i}{{\alpha }_i\left(1-{\rho }_i\right)}\right)}^2+2{\left(\frac{1-2{\rho }_i}{{\alpha }_i\left(1-{\rho }_i\right)}\right)}^2\right]}{{\left[{\displaystyle \sum _{i=1}^n}\frac{\mathsf{\text{ln}}\left[{\rho }_i\left(1-{\rho }_i\right)\right]}{{\alpha }_i}\right]}^2}$$

(12)

$${\text{Φ}}_{{\epsilon }_i}=2$$

(13)

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 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 ${\alpha }_i$ and ${\beta }_i$, and those over (0, ${\beta }_i/{\alpha }_i$) to $K_i$. 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 ${\alpha }_i$ and ${\beta }_i$ ($i\in \left\{1,2,\dots ,n\right\}$) were randomized with a range of (0, 1), while the value of $K_i$ ($i\in \left\{1,2,\dots ,n\right\}$) was randomized with a range of (0, ${\beta }_i/{\alpha }_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 ${\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 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 ${\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 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 ${\alpha }_i$ and ${\beta }_i$ ($i\in \left\{1,2,\dots ,n\right\}$) were randomized with a range of (0, 1), while the value of $K_i$ ($i\in \left\{1,2,\dots ,n\right\}$) was randomized with a range of (0, ${\beta }_i/{\alpha }_i$).

Figure 5

Examples of dynamics of the negative feedback oscillator model. A: $n=5$, B: $n=7$. ${\alpha }_i=1$, ${\beta }_i=1$ and $K_i=0.5$ ($i\in \left\{1,2,\dots ,n\right\}$).

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

$${\text{Φ}}_{\tau }=\frac{f\left(\rho \right)}{n}$$

(14)

where

$$f\left(\rho \right)=\frac{{\left(\mathsf{\text{ln}}\left[\rho \left(1-\rho \right)\right]-\frac{1-2\rho }{1-\rho }\right)}^2+2{\left(\frac{1-2\rho }{1-\rho }\right)}^2}{{\left(\mathsf{\text{ln}}\left[\rho \left(1-\rho \right)\right]\right)}^2}$$

(15)

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

Figure 6

$f\left(\rho \right)$ vs. $\rho$. $f\left(\rho \right)$ is given by Eq (15) in Text.

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

$$\tau =\mathsf{\text{ln}}4\cdot {\displaystyle \sum _{i=1}^n}{{\alpha }_i}^{-1}$$

(16)

$${\text{Φ}}_{\tau }=\frac{{\displaystyle \sum _{i=1}^n}{{\alpha }_i}^{-2}}{{\left({\displaystyle \sum _{i=1}^n}{{\alpha }_i}^{-1}\right)}^2}$$

(17)

respectively. Therefore, the value of period MPS is constrained by

$$\frac{1}{n}\le {\text{Φ}}_{\tau }\le 1$$

(18)

The minimal value of the period MPS decreases as n increases. Here, let ${\Psi }_{\tau }$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):

$$\begin{array}{cc}{\Psi }_{\tau }\hfill & ={\displaystyle \sum _{i=1}^n}\left(S_{{\alpha }_i}^{\tau }+S_{{\beta }_i}^{\tau }+S_{K_i}^{\tau }\right)\hfill \\ =\sum _{i=1}^n{S}_{{\alpha }_i}^{\tau }\hfill \\ =-1\hfill \end{array}$$

(19)

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

$$\tau =\frac{n\mathsf{\text{ln}}4}{\alpha }$$

(20)

$${\text{Φ}}_{\tau }=\frac{1}{n}$$

(21)

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 '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. ${\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 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, ${\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

$$\begin{array}{c}\frac{d{x}_1}{dt}={\beta }_1\cdot \theta \left(K_3,x_3\right)-{\alpha }_1{x}_1\\ \frac{d{x}_2}{dt}={\beta }_2\cdot \theta \left(x_1,K_1\right)-{\alpha }_2{x}_2\\ \frac{d{x}_3}{dt}={\beta }_3\cdot \theta \left(x_2,K_2\right)-{\alpha }_3{x}_3\end{array}$$

(22)

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

$$\frac{d{x}_1}{dt}={\beta }_1-{\alpha }_1{x}_1$$

(23)

Integrating Eq (23), we obtain

$$x_1\left(t\right)=\frac{{\beta }_1}{{\alpha }_1}+C_{11}{e}^{-{\alpha }_{1}t}$$

(24)

Similarly, in $I_{21}$, $I_{31}$ and $I_{12}$, $x_2$ follows

$$x_2\left(t\right)=\frac{{\beta }_2}{{\alpha }_2}+C_{21}{e}^{-{\alpha }_2\left(t-{\tau }_{11}\right)}$$

(25)

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

$$x_3\left(t\right)=\frac{{\beta }_3}{{\alpha }_3}+C_{31}{e}^{-{\alpha }_3\left(t-{\tau }_{11}-{\tau }_{21}\right)}$$

(26)

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

$$\frac{d{x}_1}{dt}=-{\alpha }_1{x}_1$$

(27)

Integrating Eq (27), we obtain

$$x_1\left(t\right)=C_{12}{e}^{-{\alpha }_1\left(t-{\tau }_{11}-{\tau }_{21}-{\tau }_{31}\right)}$$

(28)

Similarly, in $I_{22}$, $I_{32}$ and $I_{11}$, $x_2$ follows

$$x_2\left(t\right)=C_{22}{e}^{-{\alpha }_2\left(t-{\tau }_{11}-{\tau }_{21}-{\tau }_{31}-{\tau }_{12}\right)}$$

(29)

In $I_{32}$, $I_{11}$ and $I_{21}$, $x_3$ follows

$$x_3\left(t\right)=C_{32}{e}^{-{\alpha }_3\left(t-{\tau }_{11}-{\tau }_{21}-{\tau }_{31}-{\tau }_{12}-{\tau }_{22}\right)}$$

(30)

$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

$$x_i\left({\displaystyle \sum _{j=1}^i}{\tau }_{j1}\right)=K_i$$

(31)

By inserting Eqs (24)-(26), Eq (31) becomes

$$\frac{{\beta }_i}{{\alpha }_i}+C_{i1}{e}^{-{\alpha }_i{\tau }_{i1}}=K_i$$

(32)

By solving Eq (32) for ${\tau }_{i1}$, we obtain

$${\tau }_{i1}=-\frac{\mathsf{\text{ln}}{\gamma }_i}{{\alpha }_i}$$

(33)

where

$${\gamma }_i=\frac{K_i-{\beta }_i/{\alpha }_i}{C_{i1}}$$

(34)

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

$$x_i\left({\displaystyle \sum _{j=1}^3}{\tau }_{j1}+{\displaystyle \sum _{j=1}^i}{\tau }_{j2}\right)=K_i$$

(35)

By inserting Eqs (28)-(30), Eq (35) becomes

$$C_{i2}{e}^{-{\alpha }_i{\tau }_{i2}}=K_i$$

(36)

By solving Eq (36) for ${\tau }_{i2}$, we obtain

$${\tau }_{i2}=-\frac{\mathsf{\text{ln}}{\delta }_i}{{\alpha }_i}$$

(37)

where

$${\delta }_i=\frac{K_i}{C_{i2}}$$

(38)

By inserting Eqs (33) and (37) into Eqs (24)-(26), we get the peaks of $x_1$, $x_2$ and $x_3$:

$$x_1\left({\displaystyle \sum _{i=1}^3}{\tau }_{i1}\right)=\frac{{\beta }_1}{{\alpha }_1}+C_{11}{\gamma }_1{{\gamma }_2}^{{\alpha }_1/{\alpha }_2}{{\gamma }_3}^{{\alpha }_1/{\alpha }_3}$$

(39)

$$x_2\left({\displaystyle \sum _{i=1}^3}{\tau }_{i1}+{\tau }_{12}\right)=\frac{{\beta }_2}{{\alpha }_2}+C_{21}{\gamma }_2{{\gamma }_3}^{{\alpha }_2/{\alpha }_3}{{\delta }_1}^{{\alpha }_2/{\alpha }_1}$$

(40)

$$x_3\left({\displaystyle \sum _{i=1}^3}{\tau }_{i1}+{\tau }_{12}+{\tau }_{22}\right)=\frac{{\beta }_3}{{\alpha }_3}+C_{31}{\gamma }_3{{\delta }_1}^{{\alpha }_3/{\alpha }_1}{{\delta }_2}^{{\alpha }_3/{\alpha }_2}$$

(41)

By inserting Eqs (33) and (37) into Eqs (28)-(30), we get the troughs of $x_1$, $x_2$ and $x_3$:

$$x_1\left({\displaystyle \sum _{i=1}^3}{\tau }_{i1}+{\displaystyle \sum _{i=1}^3}{\tau }_{i2}\right)=K_1{{\delta }_2}^{{\alpha }_1/{\alpha }_2}{{\delta }_3}^{{\alpha }_1/{\alpha }_3}$$

(42)

$$x_2\left({\displaystyle \sum _{i=1}^3}{\tau }_{i1}+{\displaystyle \sum _{i=1}^3}{\tau }_{i2}+{\tau }_{11}\right)=K_2{{\delta }_3}^{{\alpha }_2/{\alpha }_3}{{\gamma }_1}^{{\alpha }_2/{\alpha }_1}$$

(43)

$$x_3\left({\displaystyle \sum _{i=1}^3}{\tau }_{i1}+{\displaystyle \sum _{i=1}^3}{\tau }_{i2}+{\tau }_{11}+{\tau }_{21}\right)=K_3{{\gamma }_1}^{{\alpha }_3/{\alpha }_1}{{\gamma }_2}^{{\alpha }_3/{\alpha }_2}$$

(44)

By summing up ${\tau }_{ij}$, we obtain the period:

$$\begin{array}{c}\tau ={\displaystyle \sum _{i=1}^3}\left({\tau }_{i1}+{\tau }_{i2}\right)\\ =-{\displaystyle \sum _{i=1}^3}\frac{\mathsf{\text{ln}}\left({\gamma }_i{\delta }_i\right)}{{\alpha }_i}\end{array}$$

(45)

The oscillation amplitude is calculated by subtracting the trough from the peak. From Eqs (24) and (39), the amplitude of $x_1$ is

$$\begin{array}{c}{\epsilon }_1=x_1\left({\displaystyle \sum _{i=1}^3}{\tau }_{i1}\right)-x_1\left(0\right)\\ =C_{11}\left({\gamma }_1{{\gamma }_2}^{{\alpha }_1/{\alpha }_2}{{\gamma }_3}^{{\alpha }_1/{\alpha }_3}-1\right)\end{array}$$

(46)

Similarly, from Eqs (25) and (40), the amplitude of $x_2$ is

$$\begin{array}{c}{\epsilon }_2=x_2\left({\displaystyle \sum _{i=1}^3}{\tau }_{i1}+{\tau }_{12}\right)-x_2\left({\tau }_{11}\right)\\ =C_{21}\left({\gamma }_2{{\gamma }_3}^{{\alpha }_2/{\alpha }_3}{{\delta }_1}^{{\alpha }_2/{\alpha }_1}-1\right)\end{array}$$

(47)

From Eqs (26) and (41), the amplitude of $x_3$ is

$$\begin{array}{c}{\epsilon }_3=x_3\left({\displaystyle \sum _{i=1}^3}{\tau }_{i1}+{\tau }_{12}+{\tau }_{22}\right)-x_3\left({\tau }_{11}+{\tau }_{21}\right)\\ =C_{31}\left({\gamma }_3{{\delta }_1}^{{\alpha }_3/{\alpha }_1}{{\delta }_2}^{{\alpha }_3/{\alpha }_2}-1\right)\end{array}$$

(48)

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

$$C_{11}=K_1{{\delta }_2}^{{\alpha }_1/{\alpha }_2}{{\delta }_3}^{{\alpha }_1/{\alpha }_3}-\frac{{\beta }_1}{{\alpha }_1}$$

(49)

$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),

$$C_{21}=K_2{{\delta }_3}^{{\alpha }_2/{\alpha }_3}{{\gamma }_1}^{{\alpha }_2/{\alpha }_1}-\frac{{\beta }_2}{{\alpha }_2}$$

(50)

$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),

$$C_{31}=K_3{{\gamma }_1}^{{\alpha }_3/{\alpha }_1}{{\gamma }_2}^{{\alpha }_3/{\alpha }_2}-\frac{{\beta }_3}{{\alpha }_3}$$

(51)

$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),

$$C_{12}=C_{11}{\gamma }_1{{\gamma }_2}^{{\alpha }_1/{\alpha }_2}{{\gamma }_3}^{{\alpha }_1/{\alpha }_3}+\frac{{\beta }_1}{{\alpha }_1}$$

(52)

$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),

$$C_{22}=C_{21}{\gamma }_2{{\gamma }_3}^{{\alpha }_2/{\alpha }_3}{{\delta }_1}^{{\alpha }_2/{\alpha }_1}+\frac{{\beta }_2}{{\alpha }_2}$$

(53)

$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),

$$C_{32}=C_{31}{\gamma }_3{{\delta }_1}^{{\alpha }_3/{\alpha }_1}{{\delta }_2}^{{\alpha }_3/{\alpha }_2}+\frac{{\beta }_3}{{\alpha }_3}$$

(54)

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

$$\begin{array}{c}{x}_1\left(0\right)=0\\ x_i\left({\displaystyle \sum _{j=1}^{i-1}}{\tau }_{j1}\right)=0\left(i>1\right)\end{array}$$

(55)

Using Eq (55) and Eqs (24)-(26), the integral constants $C_{i1}$ are

$$C_{i1}=-\frac{{\beta }_i}{{\alpha }_i}$$

(56)

The assumption that the peak is ${\beta }_i/{\alpha }_i$ gives

$$\begin{array}{c}{x}_1\left({\displaystyle \sum _{j=1}^n}{\tau }_{j1}\right)=\frac{{\beta }_i}{{\alpha }_i}\\ x_i\left({\displaystyle \sum _{j=1}^n}{\tau }_{j1}+{\displaystyle \sum _{j=1}^{i-1}}{\tau }_{j2}\right)=\frac{{\beta }_i}{{\alpha }_i}\left(i>1\right)\end{array}$$

(57)

Using Eqs (57) and Eqs (28)-(30), the integral constants $C_{i2}$ are

$$C_{i2}=\frac{{\beta }_i}{{\alpha }_i}$$

(58)

By inserting Eqs (56) and (58) into Eqs (34) and (38), respectively, the period $\tau$ described by Eq (3) (or Eq (45)) becomes

$$\tau =-{\displaystyle \sum _{i=1}^n}\frac{\mathsf{\text{ln}}\left[{\rho }_i\left(1-{\rho }_i\right)\right]}{{\alpha }_i}$$

(59)

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

$${\epsilon }_i=\frac{{\beta }_i}{{\alpha }_i}$$

(60)

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

$$S_{{\alpha }_i}^{\tau }=\frac{{\alpha }_i}{\tau }\frac{\partial \tau }{\partial {\alpha }_i}={\tau }^{-1}{{\alpha }_i}^{-1}\left[\mathsf{\text{ln}}\left[{\rho }_i\left(1-{\rho }_i\right)\right]-\frac{1-2{\rho }_i}{1-{\rho }_i}\right]$$

(61)

$$S_{{\beta }_i}^{\tau }=\frac{{\beta }_i}{\tau }\frac{\partial \tau }{\partial {\beta }_i}={\tau }^{-1}{{\alpha }_i}^{-1}\frac{1-2{\rho }_i}{1-{\rho }_i}$$

(62)

$$S_{K_i}^{\tau }=\frac{K_i}{\tau }\frac{\partial \tau }{\partial K_i}=-{\tau }^{-1}{{\alpha }_i}^{-1}\frac{1-2{\rho }_i}{1-{\rho }_i}$$

(63)

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

$$\begin{array}{c}{\text{Φ}}_{\tau }={\displaystyle \sum _{i=1}^n}\left[{\left(S_{{\alpha }_i}^{\tau }\right)}^2+{\left(S_{{\beta }_i}^{\tau }\right)}^2+{\left(S_{K_i}^{\tau }\right)}^2\right]\\ =\frac{{\displaystyle \sum _{i=1}^n}\left[{\left(\frac{\mathsf{\text{ln}}\left[{\rho }_i\left(1-{\rho }_i\right)\right]}{{\alpha }_i}-\frac{1-2{\rho }_i}{{\alpha }_i\left(1-{\rho }_i\right)}\right)}^2+2{\left(\frac{1-2{\rho }_i}{{\alpha }_i\left(1-{\rho }_i\right)}\right)}^2\right]}{{\left[{\displaystyle \sum _{i=1}^n}\frac{\mathsf{\text{ln}}\left[{\rho }_i\left(1-{\rho }_i\right)\right]}{{\alpha }_i}\right]}^2}\end{array}$$

(64)

The single-parameter sensitivities of amplitude are

$$S_{{\alpha }_j}^{{\epsilon }_i}=\frac{{\alpha }_j}{{\epsilon }_i}\frac{\partial {\epsilon }_i}{\partial {\alpha }_j}=\left\{\begin{array}{c}\hfill -1\begin{array}{}\end{array}\left(i=j\right)\hfill \\ \hfill 0\begin{array}{}\end{array}\left(i\ne j\right)\hfill \end{array}\right.$$

(65)

$$S_{{\beta }_j}^{{\epsilon }_i}=\frac{{\beta }_j}{{\epsilon }_i}\frac{\partial {\epsilon }_i}{\partial {\beta }_j}=\left\{\begin{array}{c}\hfill 1\begin{array}{}\end{array}\left(i=j\right)\hfill \\ \hfill 0\begin{array}{}\end{array}\left(i\ne j\right)\hfill \end{array}\right.$$

(66)

$$S_{K_j}^{{\epsilon }_i}=\frac{K_j}{{\epsilon }_i}\frac{\partial {\epsilon }_i}{\partial K_j}=0$$

(67)

Amplitude MPS is given by

$${\text{Φ}}_{{\epsilon }_i}={\displaystyle \sum _{j=1}^n}\left[{\left(S_{{\alpha }_j}^{{\epsilon }_i}\right)}^2+{\left(S_{{\beta }_j}^{{\epsilon }_i}\right)}^2+{\left(S_{K_j}^{{\epsilon }_i}\right)}^2\right]=2$$

(68)

Multiparameter sensitivity (MPS)

Generally a dynamic model for biochemical networks is described by ordinary differential equations:

$$\stackrel{·}{\mathbf{\text{x}}}=\mathbf{\text{F}}\left(t,\mathbf{\text{x}},\mathbf{\text{p}}\right)$$

(69)

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 single-parameter sensitivity of the output with respect to a change in the i th parameter is given by

$$S_{p_i}^q=\frac{p_i}{q}\frac{\partial q}{\partial p_i}=\frac{\partial \mathsf{\text{ln}}q}{\partial \mathsf{\text{ln}}{p}_i}$$

(70)

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

$${\text{Φ}}_q={\displaystyle \sum _{i=1}^m}{\left(S_{p_i}^q\right)}^2$$

(71)

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 (${\alpha }_i$, ${\beta }_i$ and $K_i$ in this study) [8]. Eq (70) can be rewritten as

$$S_{p_i}^q=\frac{p_i}{q\left(p\right)}\underset{\Delta p_i\to 0}{\mathsf{\text{lim}}}\frac{q\left(p^{\prime }\right)-q\left(p\right)}{\Delta p_i}=\underset{\Delta p_i\to 0}{\mathsf{\text{lim}}}\frac{\mathsf{\text{ln}}q\left(p^{\prime }\right)-\mathsf{\text{ln}}q\left(p\right)}{\mathsf{\text{ln}}\left(p_i+\Delta p_i\right)-\mathsf{\text{ln}}{p}_i}$$

(72)

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 single-parameter sensitivity, $\Delta p_i$ is not infinitesimal but a small value ($\Delta p_i=10^{-3}{p}_i$ 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].