A novel interaction perturbation analysis reveals a comprehensive regulatory principle underlying various biochemical oscillators

Background Biochemical oscillations play an important role in maintaining physiological and cellular homeostasis in biological systems. The frequency and amplitude of oscillations are regulated to properly adapt to environments by numerous interactions within biomolecular networks. Despite the advances in our understanding of biochemical oscillators, the relationship between the network structure of an oscillator and its regulatory function still remains unclear. To investigate such a relationship in a systematic way, we have developed a novel analysis method called interaction perturbation analysis that enables direct modulation of the strength of every interaction and evaluates its consequence on the regulatory function. We have applied this new method to the analysis of three representative types of oscillators. Results The results of interaction perturbation analysis showed different regulatory features according to the network structure of the oscillator: (1) both frequency and amplitude were seldom modulated in simple negative feedback oscillators; (2) frequency could be tuned in amplified negative feedback oscillators; (3) amplitude could be modulated in the incoherently amplified negative feedback oscillators. A further analysis of naturally-occurring biochemical oscillator models supported such different regulatory features according to their network structures. Conclusions Our results provide a clear evidence that different network structures have different regulatory features in modulating the oscillation frequency and amplitude. Our findings may help to elucidate the fundamental regulatory roles of network structures in biochemical oscillations. Electronic supplementary material The online version of this article (10.1186/s12918-017-0472-7) contains supplementary material, which is available to authorized users.


Additional Information
.
Step-by-step procedures for the interaction perturbation method to analyze an activator-amplified NFO.
Additional Figure S2. Density plots of the simple NFO.
Additional Figure S3. Density plots of the activator-amplified NFO.
Additional Figure S4. Density plots of the inhibitor-amplified NFO.
Additional Figure S5. Density plots of the type 1 incoherently-amplified NFO.
Additional Figure S6. Density plots of the type 2 incoherently-amplified NFO.
Additional Figure S7. Density plots of the type 3 incoherently-amplified NFO.

II. Additional Tables
Additional Table S1. Regulatory patterns of the frequency and amplitude according to the types of interactions Additional Table S2. Relative mean differences between 1 st -order ODE and 2 nd -order ODE in the results of the simulations during ten periods of oscillations.
Additional Table S3. Parameter sets for the 3-node biochemical oscillator models for mathematically controlled comparison Additional Table S4. Percentage of parameter sets for 3-node biochemical oscillator models that did not sustain limit-cycle oscillation after interaction perturbation from total parameter sets.

III. Additional Equations
Additional Equation A1. ODEs of the 3-node biochemical oscillators Additional Equation A2. ODEs and parameters for naturally occurring biochemical oscillator models.

IV. Additional Notes
Additional Note A1. Advantages of interaction perturbation over parameter perturbation in this study Additional Note A2. Difference between 1 st -order ODE and 2 nd -order ODE in respect of responses to a perturbation of the system

I. Additional Figures
Additional Figure A1. Step-by-step procedures for the interaction perturbation method to analyze an activator-amplified NFO. The steps in Figure A1 are an example of the steps in Figure 1.
Step 1. An activator-amplified NFO is chosen for interaction perturbation; Step 2. A parameter set is defined so that a limit cycle oscillation can be generated; Step 3. The first-order ODEs are transformed into the equivalent second-order ODEs by differentiation. The second-order ODEs are represented by the matrix product of the Jacobian matrix and first-order ODEs. Here, the Jacobian matrix shows its complete algebraic form. The element of the Jacobian matrix at the intersection of the i th row and j th column denotes the interaction from node j to node i. For example, the expression, -kdx×x, corresponds to the first row and the second column element of the Jacobian matrix. Therefore, it denotes the interaction from Y to X); Step 4. The interaction from Y to X is weakened by 2% during one period of oscillation. The weakening factor of 0.98 is applied to the expression, -kdx×x, for this perturbation; Step 5-1. The first-order ODEs are numerically integrated until a limit cycle oscillation is generated.
The time course of the variable X is plotted. The point P1 ( indicates the initial values for this integration; Step 5-2. A point is taken along a limit cycle (i.e., P2), and then the second-order ODEs are integrated using the perturbed Jacobian matrix during one period of oscillation. The point P2 is ( 1500 1500 1500 1500 1500 1500  0.1111, 0,1777, 0.0506, 0.0009, 0.0266). Although the second-order ODEs are integrated using the parameter sets which are the same as those for the integration of the first-order ODEs, the different trajectories are obtained because of the different subsets of the initial conditions (i.e., in the case of the first order ODEs: X, Y, Z; in the case of the second order ODEs: X, Y, Z, dXdt, dYdt, dZdt); Step 6.
The changes in the frequency and amplitude produced by the interaction perturbation are calculated.
Here, the frequency is increased and the amplitude is decreased.

II. Additional Tables
Additional Table A1. Regulatory patterns of the frequency and amplitude according to the types of interactions. a Both the frequency and amplitude changed by less than 1%.
b The frequency or amplitude changed by more than 1%, and the changes in the frequency are greater than the changes in the amplitude.
c The frequency or amplitude changed by more than 1%, and the changes in the amplitude are greater than the changes in the frequency.
Additional Table A2. Relative mean differences between 1 st -order ODE and 2 nd -order ODE in the results of the simulations during ten periods of oscillations. Each relative mean difference was calculated by dividing the absolute difference by the maximum value of oscillation of each oscillator.
For six representative 3-node oscillator models (*), each relative mean difference was calculated for each parameter set and represented by mean ± standard deviation.   Table A4. Percentage of parameter sets for the 3-node biochemical oscillator models that did not sustain limit-cycle oscillation after interaction perturbation from total parameter sets.

M M M P C P C PC PC PC PC B B B B I
(2) Circadian rhythm model by Goldbeter.  13.3 (1 )  ( 1 )

IV. Additional Notes Additional Note A1. Advantages of interaction perturbation over parameter perturbation in this study
The purpose of this study is to investigate the different regulatory functions of networks arising from the differences in their structure. This purpose can be achieved by examining how a network responds to a perturbation of each link within the network. Then, a question is raised as to whether a parameter perturbation in an ODE model can be regarded as a link perturbation.
The number and the function of the parameters are determined by the structure of the ODE model.
There are three parameters (i.e. ksx, Km, and kdx) in the equation (1.1) and two parameters (i.e. ksx and kdx) in the equation (1.2) although those two equations represent the same network structure which consists of two links (i.e. 'X on X' and 'Y on X'). In the equation (1.1), kdx represents 'X on X'; ksx and Km represent 'Y on X'. Herein, the perturbation of kdx clearly means the perturbation of the link 'X on X'.
However, in the link 'Y on X', the perturbation of either ksx or Km may not exactly correspond to the perturbation of the link. In the equation (1.2), kdx is involved in two links (i.e. 'X on X' and 'Y on X') simultaneously. In this structure, it is impossible to perturb a link independently of other links by parameter perturbation.
In summary, more than two parameters may represent a single link and a single parameter may be involved in more than two links. Therefore, this problem may undermine the reliability of the results of parameter perturbations.
Direct modulation of specific interaction can be a solution to this problem. All the ODE models which represent 'Y inhibits X' will consist of two interactions (i.e. X on X and Y on X). Herein 'Y inhibits X' represents the negative interaction from Y to X. What will be the results when we directly modulate the interaction from Y to X? In the equation (1.1), dX dt will increase due to the increased first In the equation (1.2), dX dt will also increase due to the decreased second term ( dx k Y X   ). Regardless of the structure of ODE models, the results of the direct perturbation of an interaction reflect the link structure inside the network. Therefore, perturbation of interactions will be methodologically suitable for analyses of the different functions of networks arising from their structure. In this study, interactions between molecules were represented algebraically by using Jacobian matrix and the functional characteristics of the networks were investigated by direct perturbation of interactions.

Additional Note A2. Difference between the first-order ODE and the second-order ODE in respect of responses to a perturbation of the system
In the first-order ODE, applying transient perturbation does not change the frequency and amplitude of the limit cycle oscillator. However, the transient perturbation in the second-order ODE can change them. This is because, for the second-order ODE, the first derivatives of the state variables as well as the state variables themselves change their values during the integration process.
Suppose that the system is defined by coupled ODEs with two state variables, X and Y, for convenience of explanation. If it is a first-order ODE system, we need two initial conditions (i.e., X(t=0) and Y(t=0)) and parameter values for numerical simulation. In this case, the simulation using identical parameter values and appropriate initial conditions would generally lead the system to converge to a same fixed point (or limit cycle). In contrast, for the simulation of an equivalent secondorder ODE system, the number of required initial conditions becomes double compared to that for the first-order ODE system (i.e., X(t=0), Y(t=0), ( 0), ( 0) dX dY t t dt dt   ). In this case, the simulation using identical parameter values might not converge to a same fixed point (or limit cycle) even if the initial conditions for X and Y (i.e., X(t=0), Y(t=0)) are identical. This is because the additionally required initial conditions (i.e., ( 0), ( 0) dX dY t t dt dt   ) are different. In Figure R1 shown below, we presented the simulation results that show a difference between a first-order ODE system and the equivalent second-order ODE system: different limit cycles are generated by integrations of the equivalent secondorder ODE system with the same parameter set and same initial values of the state variables, but with different initial time derivatives of the state variables.
These results indicate that every time when the interaction perturbation is performed, the second-order ODE system may reach a different initial condition (i.e., different initial values and different time derivatives), consequently leading to a different limit cycle.

Figure R1.
Comparison between a first-order ODE system and the equivalent second-order ODE system of a two-node activator amplified NFO. (A) Simulation results of the first-order ODE system.
Integrations from two different initial values (i.e., P1, P2) lead to the same limit cycle. (B) Simulation results of the equivalent second-order ODE system. Integrations from the same initial values (e.g. P1) can result in different limit cycles due to the different initial time derivatives of the state variables (e.g. trajectory 1 and trajectory 3: their limit cycles are different). Detailed information on the simulation is provided in Table R1.