 Research Article
 Open Access
 Published:
A novel interaction perturbation analysis reveals a comprehensive regulatory principle underlying various biochemical oscillators
BMC Systems Biology volume 11, Article number: 95 (2017)
Abstract
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 naturallyoccurring 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.
Background
Oscillations are commonly observed phenomena in biological systems and perform crucial functions in regulating physiological or cellular processes [1]. The beating of the heart, the breathing motion of the lungs, and the circadian rhythm of sleep and wakefulness can be regarded as oscillations to maintain physiological homeostasis [2]. Glucose metabolism, cyclic adenosine monophosphate (cAMP) generation, mitogenactivated protein kinase (MAPK) signaling, and cell cycle progression are wellknown cellular oscillations [3].
Oscillators appear to have different requirements for regulating the frequency and amplitude depending on their biological functions. Both frequency and amplitude of a circadian oscillator need to be regulated against fluctuations in order to maintain robust 24h rhythms [4,5,6]. In the heart beating, the frequency has to be increased according to the intensity of physical activities [7]. In neuronal firings, proper regulation of the frequency is essential for information transmission in the brain. However, both the heart beating and neuronal firings are seldom required to modulate the amplitude. On the other hand, in glycolytic and cAMP oscillators, the regulation of the amplitude is as important as the regulation of the frequency since the amplitude plays a significant role in the activities of glycolysis and the protein kinase A (PKA) signaling pathway [8, 9].
Thus, how do these oscillators meet such different requirements for regulating the frequency and amplitude? Novak et al. classified biochemical oscillators into three classes: class I oscillators (delayed negative feedback oscillators), class II oscillators (amplified negative feedback oscillators), and class III oscillators (incoherently amplified negative feedback oscillators) [10]. This classification was based solely on the network structure of the oscillator. However, interestingly, the different regulatory requirements seem to be reflected in this classification in view of the fact that (i) the circadian rhythm oscillator belongs to class I oscillators; (ii) the sinus node oscillator and neuronal oscillator belong to class II oscillators; and (iii) the glycolytic and cAMP oscillators fall into class III oscillators. Therefore, a particular type of network structure appears to serve a particular regulatory requirement better than other types, and this implies that there is an association between network structures and regulatory functions.
Such an association between them could also be inferred from the previous study by Tsai et al. in which it was revealed that an interlinked positive and negative feedback structure outperforms a simple negative feedback structure in tuning the frequency of an oscillator [11]. In addition, a positive feedback was revealed to promote the oscillation of a negative feedback oscillator [12]. However, the detailed relationship between various network structures and regulatory functions has only been partially explored till now. To investigate the relationship in detail from a systems perspective, we constructed all possible threenode oscillator models of maximum four links using ordinary differential equations (ODEs) to represent six conceptual network structures of biochemical oscillators, and then performed interaction perturbation to systematically analyze the regulatory pattern of the frequency and amplitude of each model.
So far, the parameter perturbation method has been used to study the properties of oscillators. However, this method might not be adequate to analyze the networklevel characteristics of oscillators. Because a parameter can represent various biological functions (e.g., the rate of synthesis or degradation of a molecule, the strength of binding between two molecules and the sensitivity of a reaction), perturbation of a parameter may not correspond to the variation of an interaction in the network. Moreover, the same molecular interaction can be represented in multiple ways (see Additional file 1: Notes). For instance, the interaction ‘X activates Y’ can be represented in several ways:
and
In the three equations, the parameter k _{1} represents different biological processes, and thus, the perturbation of k _{1} will yield various results. In particular, a single parameter can be involved in two interactions (eq. (2) represents two interactions, ‘X on Y’ and ‘Y on Y’), and more than two parameters can represent a single interaction (in eq. (3), k _{1} and K _{ m } are involved in the interaction ‘X on Y’). In these cases, the role of an interaction cannot be assessed independently through parameter perturbation analysis. A solution to such limitations could be to modulate the interactions rather than the parameters. For this purpose, we developed a novel perturbation strategy called the interaction perturbation method which directly modulates the strength of an interaction between two nodes in a network. By using this method, we found that a strong association exists between network structures and regulatory patterns of the frequency and amplitude of biochemical oscillators. In simple negative feedback oscillators, both the frequency and amplitude were found to be rarely modulated. In contrast, the frequency could be tuned in amplified negative feedback oscillators while the amplitude could be modulated in incoherently amplified negative feedback oscillators. Our analysis shows that different regulatory properties can emerge from different network structures of biochemical oscillators.
Results
Analyses of 3node biochemical oscillators
Because biochemical oscillator models are too diverse in their size and complexity to be investigated individually, we constructed all possible representative threenode oscillator models that consist of maximum four links. We began by determining the parameter sets and then conducted analyses of each model based on the interaction perturbation method. The procedures are provided in detail in the METHODS (Fig. 1 and Additional file 1: Figure S1).
Network structures of six 3node biochemical oscillator models
Each model includes three nodes (X, Y, and Z) and nine possible types of interactions (Lxx, Lxy, Lxz, Lyx, Lyy, Lyz, Lzx, Lzy, and Lzz) (Fig. 2). In all six models, the term X, Y, and Z denote the ‘activator’, ‘inhibitor’, and ‘mediator’, respectively, that is, X activates Y; Y inhibits X, and Z mediates the activation or inhibition.
The simple negative feedback oscillator (simple NFO) has a negative feedback loop only (Fig. 2a). For the simple NFO to be able to oscillate, at least three nodes have to be included in its negative feedback loop since the time delay required for the sustenance of the oscillation cannot be sufficiently provided with only two nodes. Adding a third node (denoted by Z in Fig. 2a) generates an appropriate time delay. In this structure, Lyx, Lzy, and Lxz form a negative feedback loop where X activates Y directly, and Y inhibits X through Z, which is consistent with the denoted function of X and Y: X as ‘activator’, and Y as ‘inhibitor’. In the activatoramplified negative feedback oscillator (activatoramplified NFO), X and Y form a twonode negative feedback loop (Lxy and Lyx), and Z amplifies X (Lxz and Lzx) (Fig. 2b). This amplification by Z plays an important role in maintaining oscillation [13]. In the inhibitoramplified negative feedback oscillator (inhibitoramplified NFO), X and Y form a twonode negative feedback loop (Lxy and Lyx), and Z amplifies Y (Lyz and Lzy) (Fig. 2c). Like the activatoramplified NFO, the amplification by Z is essential for the maintenance of oscillation [13]. Each incoherentlyamplified negative feedback oscillator (incoherentlyamplified NFO) has a negative feedback loop containing one incoherent link. We constructed three incoherentlyamplified NFOs: type 1 incoherentlyamplified NFO; type 2 incoherentlyamplified NFO; and type 3 incoherentlyamplified NFO (Fig. 2d, e, and f). Among all the types of incoherentlyamplified NFOs, Lyx, Lzy and Lxz form a negative feedback loop containing an incoherent link. An incoherent link makes the function of a node become inconsistent with its denoted function as an activator or inhibitor. For instance, in the type 1 incoherentlyamplified NFO, Y inhibits X via Z (Lzy and Lxz) and directly activates X (Lxy) simultaneously. Thus, Y is no longer an ‘inhibitor’. The incoherent link is Lxy in the type 1 and type 2 incoherentlyamplified NFO and Lyz in the type 3 incoherentlyamplified NFO.
According to the classification by Novak et al., the above six models can be classified into three classes: the simple NFO belongs to class I oscillators; the activatoramplified NFO and the inhibitoramplified NFO belong to class II oscillators; and all types of incoherentlyamplified NFOs belong to class III oscillators [10].
Analysis results of the six 3node oscillation models
To investigate the regulatory patterns of the threenode oscillators, we performed perturbations by weakening each interaction by 1%, 2%, 4%, and 8% (a weakening factor of 0.99, 0.98, 0.96, and 0.92, respectively) and observed the changes in the frequency and amplitude. The interaction perturbation was implemented by multiplying a weakening factor to the element of Jacobian matrix that is to be perturbed during one period of oscillation. Fig. 3 shows the results of the perturbations in the density plots (see the METHODS for details). In these plots, the regulatory characteristics of the frequency and amplitude are represented by the distribution patterns of the density. The concentrated density near the point (1, 1) indicates that the frequency and amplitude are robust to perturbations. The horizontal distribution of the density denotes that the change in the frequency is larger than the change in the amplitude, and the vertical distribution of the density denotes the opposite.
To represent the results quantitatively, we grouped the changes in the frequency and amplitude into three patterns: In pattern R, both the frequency and amplitude changed by less than 1%; in pattern F, either the frequency or the amplitude changed by more than 1% and the changes in the frequency were greater than the changes in the amplitude; in pattern A, either the frequency or the amplitude changed by more than 1% and the changes in the amplitude were greater than the changes in the frequency. Table 1 shows the distribution of the patterns in each 3node oscillator model (a full description of the distribution of the patterns generated by the perturbation of each interaction is provided in Additional file 1: Table S1).
The simple NFO showed the highest robustness to perturbations regardless of the types of perturbed interactions (Fig. 3a and Additional file 1: Figure S2). The rates of change in both the frequency and amplitude were less than 1% in 92.9% of the perturbation results (Table 1) and are depicted as the darkest density concentrated on the point (1, 1) (Fig. 3a).
In the activatoramplified NFO and inhibitoramplified NFO, the change in the frequency was larger than the change in the amplitude. The results are shown in Fig. 3b and c, in which the density is spread in a nearly horizontal direction. These changes in the frequency were caused by the perturbations of Lxx or Lxy (Additional file 1: Figure S3 and Additional file 1: Figure S4). In both oscillators, pattern F was observed in more than 30% of the perturbation results.
In contrast to the regulatory patterns observed in the activatoramplified NFO and the inhibitor amplified NFO, all the incoherentlyamplified NFOs showed moderate changes in the amplitude. In the type 1 incoherentlyamplified NFO, the amplitude was slightly more adjustable to perturbations than the frequency while in the type 2 incoherentlyamplified NFO, the amplitude was changed to a large extent (Fig. 3d and e). In the type 3 incoherentlyamplified NFO, the frequency and amplitude were changed to various extents (Fig. 3f). Both in the type 2 incoherentlyamplified NFO and the type 3 incoherentlyamplified NFO, pattern A was observed in more than 50% of the perturbation results.
For the six 3node oscillator models, the perturbation results obtained by weakening each interaction by 2%, 4%, or 8% had qualitatively the same regulatory patterns as those obtained by weakening each interaction by 1% (Additional file 1: Figures S2S7).
In summary, a distinct regulatory pattern was observed in each 3node oscillator. Class I oscillator (the simple NFO) is robust to perturbations while for class II oscillators (the activatoramplified NFO and the inhibitoramplified NFO), the frequency can be selectively regulated. In class III oscillators (types 1, 2, and 3 incoherentlyamplified NFOs), the amplitude can be regulated. Based on these observations, we deduced the regulatory principle that the differences in network structures give rise to different regulatory patterns of the frequency and amplitude.
Mathematically controlled comparisons between structurally related biochemical oscillators
Class I, class II and class III oscillators are structurally related to one another, and their structural differences arise from an added link. A class II oscillator can be formed by adding a link to the activator (X) or inhibitor (Y) of a class I oscillator. A class III oscillator can be formed by adding an incoherent link to a class I oscillator.
This prompted us to assume that the added links could bring about changes to the regulatory patterns of the oscillators. To examine this idea, we performed mathematically controlled comparisons between structurally related oscillators.
Construction of structurally related threenode models for mathematically controlled comparisons
We developed three additional threenode oscillator models which contain one additional link to the backbone of a simple NFO. A selfactivating positive feedback link was added to the activator (X) and inhibitor (Y) of the simple NFO to generate an activatoramplified NFO variant and an inhibitoramplified NFO variant, both of which belong to class II oscillators. Adding Lxy to the simple NFO generated a variant of the type 1 incoherentlyamplified NFO (hereinafter called the type 1 incoherentlyamplified NFO variant), which belongs to class III oscillators. Simulation of the simple NFO was performed with a representative parameter set suggested by Novak et al. [10]. The parameters of the newly generated oscillators (the activatoramplified NFO variant, the inhibitoramplified NFO variant and the type 1 incoherentlyamplified NFO variant) were determined with methods for mathematically controlled comparisons [14]. The full ODEs are provided in Additional file 1: Eq. A1, and the full parameters are provided in Additional file 1: Table S3.
Analysis results of the structurally related models
Perturbations on the oscillators were performed by weakening each interaction by 1% during one period of oscillation. A distinct regulatory pattern for each model could be identified despite the fact that the frequency and amplitude changed by less than 1% compared to the unperturbed cases in all four models (Fig. 4). In the activatoramplified NFO variant and the inhibitoramplified NFO variant, the frequency was more adjustable than the amplitude, whereas in the type 1 incoherentlyamplified NFO variant, the amplitude was more adjustable than the frequency. Overall, adding an amplifying link could enhance the ability to regulate the frequency of the oscillator whereas adding an incoherent link could enhance the ability to regulate the amplitude of the oscillator.
Analyses of naturallyoccurring biochemical oscillator models
To explore whether the regulatory principle suggested here could also apply to naturallyoccurring biochemical oscillator models, we performed analyses of nine wellknown biochemical oscillator models which were constructed based on experimental data. For each model, perturbations were conducted by weakening each interaction by 1% during oneperiod of oscillation. The subsequent changes in the frequency and amplitude are shown in Fig. 5. The ODE equations and the parameters are provided in the Additional file 1: Eq. A2.
The circadian rhythm model by Goldbeter [15], the circadian rhythm model by Leloup et al. [16] and the repressilator [17] are wellknown examples of class I oscillators (a class I oscillator consists of a negative feedback only). These oscillators showed the highest robustness to perturbations: both the frequency and amplitude rarely changed in response to the perturbations. The sinus node model by Yanagihara et al. [18], the neuronal excitation model by Hodgkin and Huxley [19] and the cell cycle model by Pomerening et al. [20] can be classified as class II oscillators (a class II oscillator includes a positive feedback). These class II oscillators worked better in adjusting the frequency than in adjusting the amplitude; the perturbations induced more changes to the frequency than to the amplitude. On the other hand, in the cAMP oscillator model [21] and the glycolysis models [22, 23] that belong to class III oscillators (a class III oscillator includes an incoherent link), the amplitude was more adjustable than the frequency: the amplitude changed more than the frequency.
In summary, the regulatory principle suggested here in the threenode oscillator models could also apply to naturallyoccurring biochemical oscillator models.
Discussion
Our analysis based on the interaction perturbation method revealed the regulatory principle that different network structures of biochemical oscillators give rise to different regulatory patterns of the frequency and amplitude; for class I oscillators, the frequency and amplitude are seldom regulated; for class II oscillators, the frequency is more adjustable than the amplitude; for class III oscillators, the amplitude is more adjustable than the frequency. The results of the mathematically controlled comparisons further demonstrated the reliability of this regulatory principle and its potential for application to naturallyoccurring biochemical oscillator models.
In systems biological studies, the parameter perturbation method has been widely used to investigate the relationship between network structures and their biological functions [24,25,26,27,28,29]. In addition, various mathematical methods have been developed to analyze the characteristics of oscillators. The sensitivity heat map and parameter sensitivity spectrum developed by Rand et al. have been utilized to provide a more integrated picture of the overall sensitivities of a system and to probe how the function of a network depends upon its structure and parameters [30]. Irene et al. proposed an optimizationbased approach to investigate what environmental conditions drive specific oscillatory network [31]. The state sensitivity decomposition method developed by Wilkins et al. is useful in analyzing the influence of parameter changes on period, amplitude and relative phase of oscillation [32]. In addition, robustness and dynamical characteristics between various oscillatory systems could be effectively compared by using bifurcation analysis, parameter sensitivity analysis, and stochastic simulation [33, 34]. However, since all these methods somehow focus on “parameters”, they can be classified as a parameter perturbation analysis in a broad sense. On the other hand, as far as the network topology is concerned, the biological meaning of a parameter does not always correspond to specific interaction. Hence, there is still difficulty in attributing the perturbation of a particular interaction in a regulatory network to the perturbation of a parameter in the corresponding mathematical model.
The interaction perturbation method proposed in this study has several advantages over the parameter perturbation method. First, the result of an interaction perturbation analysis can be properly interpreted in the context of a biological network since the perturbation directly modulates a link of the network structure. Second, this method can provide a more pertinent comparison between different network structures by allowing the focus of the comparison to be placed on the difference of the interaction in the network, not on the indirect difference of the underlying biological process. If the network structures have the same number of nodes, the comparison can be performed more effectively between them as they have a Jacobian matrix of the same size (threenode networks have a Jacobian matrix of 3 by 3). Third, this method can simplify analysis procedures. Previously suggested methods (e.g., optimizationbased method, state sensitivity decomposition method, etc) for the analysis of biochemical oscillators can provide meaningful insights into the nature of oscillators, but most of them require a certain level of expert knowledge on mathematics [31, 32]. In contrast, with a given parameter set, we just need to transform firstorder ODEs into secondorder ODEs and integrate the secondorder ODEs using a perturbed Jacobian matrix without going through any other complicated procedures such as selection of a bifurcation parameter, identification of a Hopf bifurcation point and numerical continuation [35].
In this study, we demonstrated that the regulatory patterns of the frequency and amplitude depend on the network structures of the biochemical oscillators. Notably, even for the same class of network structures, different regulatory patterns were observed. For instance, for the activatoramplified NFO, the amplitude was adjustable although the range was narrow whereas for the inhibitoramplified NFO, modulation of the amplitude was negligible. The regulatory range of the amplitude was wider for the type 2 incoherentlyamplified NFO than that for the type 1 or a type 3 incoherentlyamplified NFO. For the type 3 incoherentlyamplified NFO, both the amplitude and frequency could be regulated to a various extent. Thus, not only overall network topologies but the interlinkage of nodes appear to be involved in the formation of regulatory patterns of the frequency and amplitude.
It may also be noteworthy to mention that, in this study, the chosen parameter set and the kind of interactions were identified as a minor contributory factor that could affect the regulatory patterns of the frequency and amplitude, though not significantly. For class I oscillators, the frequency and amplitude were changed by less than 1% for most of the parameter sets (pattern R) except for a few parameter sets where the amplitude was changed by more than 1% (pattern A). For class II oscillators, the frequency was adjustable for more than one third of the parameter sets (pattern F) whereas, for the others, the frequency was not changed (patterns R and A). For class III oscillators, the amplitude was adjustable for a relatively greater part of the parameter sets (pattern A) whereas, for the others, the amplitude was not changed (patterns R and F).
For the activatoramplified NFO and inhibitoramplified NFO, by perturbation of Lxx or Lxy, the frequency was adjusted whereas the perturbations did not significantly change the frequency or amplitude. In the incoherently amplified NFOs, only some kinds of interactions seem to be involved in modulation of the amplitude.
Our analyses of naturallyoccurring biochemical oscillator models showed that the regulatory principle suggested here may have applications in naturallyoccurring biochemical oscillation models. A question then might arise as to what functional benefits can be derived from a particular network structure in naturallyoccurring biochemical oscillators? Both frequency and amplitude of a circadian oscillator need to be regulated against fluctuations in order to maintain robust 24h rhythms [36, 37]. This requirement can be satisfied by the network structure of a class I oscillator. When an incoherent feedforward structure is added to such a circadian oscillator, a stable oscillator with a different frequency can be generated and used to meet other biological needs [38]. Sinus nodal cells and neurons should be able to tune their frequency to transmit the information to neighboring cells appropriately [39] and cell cycles have to regulate the rate of their progression appropriately in response to environmental changes. To this end, the network structure of a class II oscillator might be a suitable one. In the glycolysis and cAMP models, regulation of the amplitude has greater importance since the amplitude of phosphofructokinase and cAMP, which are actively involved in metabolic and signaling processes respectively, has a significant role in cellular functions. Therefore, the network structure of a class III oscillator can be a better choice for these cases.
From an evolutionary point of view, the fact that class II and class III oscillators are frequently found in natural biological examples can be considered as important evidence implicating that they might have some performance advantages over class I oscillators. This study suggests that such advantages might be found from the ability to tune the frequency in class II oscillators and the ability to regulate the amplitude in class III oscillators.
Our study has the following limitations. Firstly, our study employed the Jacobian matrix which describes the interactions between state variables instead of employing the monodromy matrix (i.e., the state transition matrix over one period) or Floquet multipliers (i.e., eigenvalues of the monodromy matrix) which have been widely used to determine the oscillatory properties of a system with a limit cycle. Therefore, our scope of analysis was largely confined to examining influences of interaction perturbation on the oscillatory properties. Secondly, interaction perturbation was performed during one cycle of oscillation because longer duration of perturbation destabilized the oscillation in most cases. Thirdly, our analysis of structurally related models may not be sufficient to investigate the general characteristics of each structure in greater detail since it was performed under the predefined parameter combinations.
Conclusions
From the analyses based on the interaction perturbation method, we found a new regulatory principle that differences in network structures can give rise to different regulatory patterns of the frequency and amplitude. This finding could serve as a basis for further investigation into the underlying mechanism for the regulation of the frequency and amplitude in existing biochemical oscillators as well as for designing synthetic oscillators with a specific regulatory function.
Methods
Analysis procedures for 3node biochemical oscillators
We constructed six representative oscillator models and generated random parameter sets for each model ensuring its sustained oscillation under the parameter sets. Then, we converted the interactions in the model into corresponding elements of the Jacobian matrix and performed perturbations of the elements. After the perturbations, we measured resulting changes of the frequency and amplitude. The analysis workflows are described in Fig. 1 and Additional file 1: Figure S1.
Construction of six representative models for biochemical oscillators
Each threenode oscillator model was described in terms of three coupled ODEs with the combinatorial use of mass action laws and MichaelisMenten kinetics. The ODEs of the six 3node oscillator models are provided in the Additional file 1: Eq. A1. In every oscillator model, the oscillations were sustained under specific parameter sets. After an initial transient, integrations that started under different initial conditions quickly converged to a common trajectory with the same frequency and amplitude, namely, a limitcycle oscillation.
Random parameter generation for the six 3node biochemical oscillator models
Determining the parameter values constitutes an important process to create sustained oscillations. We chose to determine parameter values by extensive search of the parameter space because an analytical approach does not lend itself to dealing with a large number of parameters.
Random parameter sets were generated for each threenode oscillator model by selecting parameters from an exponential distribution within the range of 0.001 to 1000, using the Latin hypercube sampling method [40]. This range corresponds to biologically reasonable values typically used to model biological systems [41,42,43,44,45]. All parameters except for the Hill coefficients were randomly generated. For each parameter set, we verified whether the model produced a limitcycle oscillation [46]. Through this process, a total of 1000 parameter sets were generated for each model, and consequently, each threenode oscillator model yielded 1000 parametervalueassigned models.
Algebraic representation of interaction using Jacobian matrix
A network topology shows clearly whether an interaction is activating or inhibiting. However, when this network is represented by the system of ODEs, the function of interaction is not easily identifiable. To represent an interaction as an algebraic object, the Jacobian matrix will be a reasonable choice since an element of the Jacobian matrix corresponds to an interaction. Because an element of the Jacobian matrix cannot be obtained directly from firstorder ODEs, we differentiated the firstorder ODEs with respect to time to generate secondorder ODE systems (step 3 in Fig. 1). Thus, this system can be represented by the matrix product of the Jacobian matrix and firstorder ODEs. Here, the Jacobian matrix is defined as a _{ ij } = ∂f _{ i }/∂x _{ j }, where i and j denote the row and column indices of the Jacobian matrix, respectively. A nonzero value of a _{ ij }means that variable x _{ j } influences the evolution of variable x _{ i }; in other words, an interaction from node j to node i exists [47].
Firstorder ODEs were numerically integrated using various initial conditions until they converged to a limitcycle oscillation. In order to confirm that secondorder ODEs are good approximations to firstorder ODEs, we simulated secondorder ODEs and firstorder ODEs using a point in a limitcycle oscillation as an initial point. The results showed almost no differences (see Additional file 1: Table S2 for the mean differences between the two simulations).
Perturbation conditions
Because the aim of this study is to investigate the difference in the regulatory patterns of the frequency and amplitude between network structures, the following perturbation conditions for each parametervalueassigned model were set so that the regulatory patterns should not be influenced by perturbation conditions: the type of interaction to be perturbed; the perturbation strength; the perturbation duration; and the perturbation starting point.
Every type of interaction was perturbed one by one since the function of an interaction may not be distinguishable if more than two interactions were perturbed simultaneously. The strength of the perturbations was weakened by 1%, 2%, 4%, and 8%, which correspond to weakening factors of 0.99, 0.98, 0.96, and 0.92, respectively. These weakening degrees were selected because weakening perturbations by more than 8% often destabilized the limit cycle oscillation.
After determining the type of interactions and weakening factors, we multiplied a weakening factor by the corresponding element of the Jacobian matrix to construct the Jacobian matrix of the perturbation. For instance, when we perturbed Lyx with a weakening factor of 0.99, we multiplied 0.99 to the element at the second row and first column of the Jacobian matrix leaving all other remaining elements the same.
The perturbations were performed during one period of oscillation since the results could vary according to the oscillatory phases. We established 40 starting points of perturbations which were evenly distributed along the time cycle to prevent the trajectories from being influenced by the positions where the perturbations began.
Perturbation processes
A firstorder ODE was numerically integrated using various initial conditions until it reached a starting point of perturbation. We simulated a secondorder ODE using the perturbed Jacobian matrix during one period of oscillation. After that, the secondorder ODE was integrated using the unperturbed Jacobian matrix until the oscillation was stabilized.
Representation of perturbation results
We measured the frequency and amplitude when the oscillation was stabilized after completion of a perturbation while excluding the case of damped oscillations. The percentages of the parameter sets out of total parameter sets for 3node oscillators that did not show sustained limitcycle oscillations are provided in Additional file 1: Table S4. The change in the frequency and amplitude is presented as a ratio to the frequency and amplitude of the oscillation before the perturbation. For instance, the ratio (2, 0.5) means that the frequency doubled and the amplitude halved.
The perturbations of a threenode oscillator model generated around 1000,000 results. To make the results more intuitive and easily understood, we depicted them in density plots. We divided the whole frequencyamplitude domain into 10,000 equalsized subdomains, and calculated the number of results that belong to each subdomain and the percentage that it occupies in the total number of results, and then converted the calculated percentage into density and each colorcoded subdomain according to its density; in the plots, the density is increasing over a continuum starting from white followed by yellow, red and black.
Mathematically controlled comparisons
To determine the parameter sets for the activatoramplified NFO variant, the inhibitoramplified NFO variant, and the type 1 incoherentlyamplified NFO variant, we adopted a method for mathematically controlled comparisons proposed by Michael A. Savageau [14]. First, the values of the parameters for the unaltered processes of one system are assumed to be identical with those of the corresponding parameters of the other system. For instance, in the activatoramplified NFO variant, degradation of X, synthesis of Y, degradation of Y, synthesis of Z, and degradation of Z are unaltered processes in comparison to the simple NFO. So, parameters for those processes (k_{dx}, k_{1}, k_{2}, K_{m}, and k_{3}) in the activatoramplified NFO variant are equal to the corresponding parameters of the simple NFO. Second, parameters associated with altered processes are free to assume any values. In the activatoramplified NFO variant, parameters related with the synthesis of X are assumed to have any values. Third, the parameters of the altered processes are determined by imposing constraints on the external behavior of the system. For the above three oscillation models, the following two constraints were imposed to determine the free parameters: (i) integration of the ODE models with the specified parameter sets have to be able to generate a limitcycle oscillation; (ii) the frequency and amplitude of the oscillation have to be similar to those of the simple NFO.
After determination of the parameters, perturbations were performed in the same manner as in the conceptual threenode models. Then, the changes in the frequency and amplitude in the four oscillator models (the simple NFO, the activatoramplified NFO variant, the inhibitoramplified NFO variant, and the type 1 incoherentlyamplified NFO variant) were compared.
Abbreviations
 cAMP:

cyclic adenosine monophosphate
 MAPK:

mitogenactivated protein kinase
 NFO:

negative feedback oscillator
 ODE:

ordinary differential equation
 PKA:

protein kinase A
References
 1.
Hess B, Boiteux A. Oscillatory phenomena in biochemistry. Annu Rev Biochem. 1971;40:237–58.
 2.
Glass L. Synchronization and rhythmic processes in physiology. Nature. 2001;410:277–84. Nature Publishing Group
 3.
Maini P. Biochemical oscillations and cellular rhythms: the molecular basis of periodic and chaotic behaviour by Albert Goldbeter, Cambridge University press, 1996. ISBN 0 521 40307 3. Trends Biochem Sci. 1996;21:403.
 4.
Reppert SM, Weaver DR. Coordination of circadian timing in mammals. Nature. 2002;418:935–41.
 5.
Dunlap JC. Molecular bases for circadian clocks. Cell. 1999;96:271–90. Elsevier
 6.
Rand DA, Shulgin BV, Salazar D, Millar AJ. Design principles underlying circadian clocks. J R Soc Interface. 2004;1:119–30. The Royal Society
 7.
Hall JE, Guyton AC. The circulation. Guyton and hall physiology review. Philadelphia: Elsevier. 2011; p.41–70.
 8.
Tasken K, Aandahl EM. Localized effects of cAMP mediated by distinct routes of protein kinase a. Physiol Rev. 2004;84:137–67.
 9.
Pilkis SJ, Granner DK. Molecular physiology of the regulation of hepatic gluconeogenesis and glycolysis. Annu Rev Physiol. 1992;54:885–909.
 10.
Novák B, Tyson JJ. Design principles of biochemical oscillators. Nat Rev Mol Cell Biol. 2008;9:981–91. Nature Publishing Group
 11.
Tsai TYC, Choi YS, Ma W, Pomerening JR, Tang C, Ferrell JE. Robust, tunable biological oscillations from interlinked positive and negative feedback loops. Science. 2008;321:126–9.
 12.
Ananthasubramaniam B, Herzel H. Positive feedback promotes oscillations in negative feedback loops. PLoS One. 2014;9:e104761. Thattai M, editor
 13.
Ferrell JE, Tsai TYC, Yang Q. Modeling the cell cycle: why do certain circuits oscillate? Cell. 2011;144:874–85. Elsevier
 14.
Savageau MA. Design principles for elementary gene circuits: elements, methods, and examples. Chaos. 2001;11:142. American Institute of Physics AIP
 15.
Goldbeter A. A model for circadian oscillations in the drosophila period protein (PER). Proc R Soc B Biol Sci. 1995;261:319–24.
 16.
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 Rhythm. 1998;13:70–87. Sage Publications Sage CA: Thousand Oaks, CA
 17.
Elowitz MB, Leibler S. A synthetic oscillatory network of transcriptional regulators. Nature. 2000;403:335–8. Nature Publishing Group
 18.
Yanagihara K, Akinori N, Irisawa H. Reconstruction of sinoatrial node pacemaker potential based on the voltage clamp experiments. Jpn J Physiol. 1980;30:841–57.
 19.
Hodgkin AL, Huxley AF. A quantitative description of membrane current and its application to conduction and excitation in nerve. J Physiol. 1952;117:500–44.
 20.
Pomerening JR, Kim SY, Ferrell JE. Systemslevel dissection of the cellcycle oscillator: bypassing positive feedback produces damped oscillations. Cell. 2005;122:565–78. Elsevier
 21.
Martiel JL, Goldbeter A. A model based on receptor desensitization for cyclic AMP signaling in Dictyostelium cells. Biophys J. 1987;52:807–28. Elsevier
 22.
Sel'kov EE. Selfoscillations in Glycolysis. 1. A simple kinetic model. Eur J Biochem. 1968;4:79–86. Blackwell Publishing Ltd
 23.
Higgins J. A chemical mechanism for oscillation of glycolytic intermediates in yeast cells. Proc Natl Acad Sci. 1964;51:989–94.
 24.
Kim D, Rath O, Kolch W, Cho KH. A hidden oncogenic positive feedback loop caused by crosstalk between Wnt and ERK pathways. Oncogene. 2007;26:4571–9. Nature Publishing Group
 25.
Shin SY, Rath O, Zebisch A, Choo SM, Kolch W, Cho KH. Functional roles of multiple feedback loops in extracellular signalregulated kinase and Wnt signaling pathways that regulate epithelialMesenchymal transition. Cancer Res. 2010;70:6715–24.
 26.
Shin SY, Yang HW, Kim JR, Do Heo W, Cho KH. A hidden incoherent switch regulates RCAN1 in the calcineurin–NFAT signaling network. J Cell Sci. 2011;124:82–90. The Company of Biologists Ltd
 27.
Kwon YK, Cho KH. Coherent coupling of feedback loops: a design principle of cell signaling networks. Bioinformatics. 2008;24:1926–32. Oxford University Press
 28.
Shin SY, Choo SM, Kim D, Baek SJ, Wolkenhauer O, Cho KH. Switching feedback mechanisms realize the dual role of MCIP in the regulation of calcineurin activity. FEBS Lett. 2006;580:5965–73.
 29.
Lee HS, Hwang CY, Shin SY, Kwon KS, Cho KH. MLK3 is part of a feedback mechanism that regulates different cellular responses to reactive oxygen species. Sci Signal. 2014;7:ra52.
 30.
Rand DA. Mapping global sensitivity of cellular network dynamics: sensitivity heat maps and a global summation law. J R Soc Interface. 2008;5(Suppl 1):S59–69.
 31.
OteroMuras I, Banga JR. Design principles of biological oscillators through optimization: forward and reverse analysis. PLoS One. 2016;11:e0166867. Poyatos JF, editor.
 32.
Wilkins AK, Tidor B, White J, Barton PI. Sensitivity analysis for oscillating dynamical systems. SIAM J Sci Comput. 2009;31:2706–32.
 33.
CaicedoCasso A, Kang HW, Lim S, Hong CI. Robustness and period sensitivity analysis of minimal models for biochemical oscillators. Sci Rep. 2015;5:13161. Nature Publishing Group
 34.
CaicedoCasso A, Kang HW, Lim S, Hong CI. Corrigendum: robustness and period sensitivity analysis of minimal models for biochemical oscillators. Sci Rep. 2016;6:18504. Nature Publishing Group
 35.
Rinaldi S, Muratori S, Kuznetsov Y. Multiple attractors, catastrophes and chaos in seasonally perturbed predatorprey communities. Bull Math Biol. 1993;55:15–35.
 36.
Hastings JW, Sweeney BM. On the mechanism of temperature independence in a biological clock. Proc Natl Acad Sci. 1957;43:804–11.
 37.
Pittendrigh CS. On temperature independence in the clock system controlling emergence time in drosophila. Proc Natl Acad Sci. 1954;40:1018–29.
 38.
Martins BM, Das AK, Antunes L, Locke JC. Frequency doubling in the cyanobacterial circadian clock. Mol Syst Biol. 2016;12:896.
 39.
Izhikevich EM, Desai NS, Walcott EC, Hoppensteadt FC. Bursts as a unit of neural information: selective communication via resonance. Trends Neurosci. 2003;26:161–7.
 40.
Dubitzky W, Wolkenhauer O, Cho KH, Yokota H, editors. Latin hypercube sampling. Encyclopedia of systems biology. New York, NY: Springer; 2013. p. 1105.
 41.
Chickarmane V, Kholodenko BN, Sauro HM. Oscillatory dynamics arising from competitive inhibition and multisite phosphorylation. J Theor Biol. 2007;244:68–76.
 42.
Liu P, Kevrekidis IG, Shvartsman SY. Substratedependent control of ERK phosphorylation can lead to oscillations. Biophys J. 2011;101:2572–81. Elsevier
 43.
Markevich NI, Hoek JB, Kholodenko BN. Signaling switches and bistability arising from multisite phosphorylation in protein kinase cascades. J Cell Biol. 2004;164:353–9. Rockefeller University Press
 44.
Qiao L, Nachbar RB, Kevrekidis IG, Shvartsman SY. Bistability and oscillations in the HuangFerrell model of MAPK signaling. PLoS Comput Biol. Public Library of Science. 2007;3:e184.
 45.
Shankaran H, Ippolito DL, Chrisler WB, Resat H, Bollinger N, Opresko LK, et al. Rapid and sustained nuclear–cytoplasmic ERK oscillations induced by epidermal growth factor. Mol Syst Biol. EMBO Press. 2009;5:332.
 46.
Dhooge A, Govaerts W, Kuznetsov YA. MATCONT: a MATLAB package for numerical bifurcation analysis of ODEs. ACM transactions on mathematical software (TOMS). ACM. 2003;29:141–64.
 47.
Thomas R. Circular causality. IEE Proc Syst Biol. IET Digital Library. 2006;153:140–53.
Acknowledgements
We would like to thank HoSung Lee and JeHoon Song for their valuable discussions on the initial manuscript.
Funding
This work was supported by the National Research Foundation of Korea (NRF) grants funded by the Korea Government, the Ministry of Science, ICT & Future Planning (2017R1A2A1A17069642, 2015M3A9A7067220, and 2013M3A9A7046303). It was also supported by the KAIST Future Systems Healthcare Project from the Ministry of Science, ICT & Future Planning, the KUSTARKAIST Institute, Korea under the R&D program supervised by the KAIST, and the grant of the Korean Health Technology R&D Project, Ministry of Health & Welfare, Republic of Korea (HI13C2162).
Availability of data and materials
All model equations and parameters for the replication of the results are provided in the Additional file 1.
Author information
Affiliations
Contributions
KHC supervised this study; JHK and KHC designed the study; JHK and KHC conducted modeling, simulations, and analysis of the simulation results; JHK and KHC wrote the manuscript. All authors read and approved the final manuscript.
Corresponding author
Correspondence to KwangHyun Cho.
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Additional file
Additional file 1:
Details on methods and results of the interaction perturbation analysis for comprehensive assessment of the regulatory principle underlying various biochemical oscillators. (PDF 3894 kb)
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Kang, J., Cho, K. A novel interaction perturbation analysis reveals a comprehensive regulatory principle underlying various biochemical oscillators. BMC Syst Biol 11, 95 (2017) doi:10.1186/s1291801704727
Received
Accepted
Published
DOI
Keywords
 Biochemical oscillators
 Network structure
 Regulation of frequency and amplitude
 Perturbation analysis
 Systems biology