Phase Responses of Oscillating Components in a Signaling Pathway

Signal transduction pathways control various events in mammalian cells such as growth, proliferation, differentiation, apoptosis, or migration in response to environmental stimuli. Because of their importance, the activity of signaling pathways is controlled by multiple modes of positive and negative feedback regulation. Although negative feedback regulation primarily functions to stabilize a system, it also becomes a source of emerging oscillations. For example, the oscillatory behavior of mitogen-activated protein kinase (MAPK) activity has been theoretically proposed earlier and experimentally verified recently. However, the physiological function of such oscillatory behavior in biological systems remains unclear. To understand the functional aspects of this behavior, one should analyze the oscillation dynamics from a mathematical point of view. In this study, we applied the phase reduction method to two simple, structurally similar phosphorylation-dephosphorylation cycle models with negative feedback loops (Models A and B) and a MAPK cascade model, whose dynamics all show oscillation. We found that all three models we tested have a Type II phase response. In addition, we found that when a pair of each models A and B coupled through a weak diffusion interaction, they could synchronize with a zero phase difference. A pair of MAPK cascade models also showed synchronous oscillation, however, when a time delay was introduced into the coupling, it showed an asynchronous response. These results imply that structurally similar or even identical biological oscillators can produce differentiated dynamics in response to external perturbations when the cellular environment is altered. Synchronous or asynchronous oscillation may add strength to or dampen the efficiency of signal propagation, depending on subcellular distances and cell density. Phase response analysis allows prediction of dynamics changes in oscillations in multiple cellular environments.

1. Introduction

Oscillatory dynamics are widely distributed in nature (Strogatz, 2003). In biological systems, circadian rhythm, heart rhythm, locomotion, and electrical activity of the brain are well-known oscillation generators (Winfree, 2001; Khalsa et al., 2003). Although oscillatory behavior is a product of negative feedback regulation, the question of how the oscillatory information is processed in biological systems is still unresolved. Mammalian cells respond to extra-cellular signals and transfer this information to the nucleus to express/repress genes necessary for adaptation to a new environment or differentiation state. Signal transduction pathways play important roles to control expression of the correct genes and with the precise timing to satisfy cellular needs. Therefore, signaling pathways are spatio-temporally controlled by many positive and negative feedback loops through transcriptional and post-transcriptional modification. As a result, several types of oscillatory behaviors in the components of signaling pathways can be observed when negative feedback regulation is introduced into a system. For example, Shankaran et al. have shown persistent periodic shuffling of fluorescent-labeled extra-cellular signal-regulated kinase (ERK, a subset of the family of mitogen-activated protein kinases, MAPK) between cytosol and nucleus in epidermal growth factor (EGF) stimulated cells at the single cell level. Intriguingly, these periodic cycles among neighboring cells were asynchronous (Shankaran et al., 2009). ERK is one of the deterministic kinases that control transcription when translocated into the nucleus, therefore this nuclear shuffling process is highly regulated. The work of Shankaran et al. was the first experimental demonstration of the oscillatory behavior of ERK, although this was predicted earlier on theoretical grounds (Kholodenko, 2006). Given this asynchronous oscillation, one would think that it would be difficult to identify ERK dynamics in a population of cells, where the signal would be averaged. In addition, since periodic activation of ERK has been difficult to demonstrate experimentally, the cellular conditions leading to oscillatory ERK activation are likely quite narrow and restricted. In general, when oscillators interact with each other through a strong coupling, they tend to synchronize. Therefore, the asynchronous oscillation observed by Shankaran et al. suggests that the coupling strength of ERK in neighboring cells is weak, at least under the experimental conditions used in these studies. The question then arises, what kind of conditions allows a pair of cells to achieve asynchronous oscillation? Is the weak coupling enough to cause asynchronous oscillation? We have investigated these questions by applying the phase reduction method to two models of phosphorylation-dephosphorylation cycles and in a MAPK cascade, all of which exhibit negative feedback regulation.

The mechanism that causes the emergence of oscillations has been energetically studied (for example, Guckenheimer and Holmes, 1983; Kuramoto, 2003). Theoretically, oscillators are classified according to their bifurcation types, such as saddle-node and Hopf bifurcations. Oscillation by the saddle-node bifurcation emerges when a half-stable cycle splits into a pair of stable and unstable limit cycles, but oscillation by the Hopf bifurcation emerges when a stable spiral fixed point changes to a unstable spiral fixed point surrounded by a stable limit cycle (Hale and Koçak, 1991; Strogatz, 2001). Although phase space structures can be partly derived from such bifurcation types, the dynamic properties of the system have to be evaluated by other methods. To investigate the underlying oscillatory mechanism, a framework termed the phase reduction method has been developed in mathematics and non-linear physics (Hansel et al., 1995; Hoppensteadt and Izhikevich, 1997; Kuramoto, 2003). By using this method, an oscillator in a high dimensional space can be described by only one variable, phase, and its dynamics are packed into a phase response (or phase sensitivity) function. The phase response function has been derived analytically, not only from mathematical models but also from experimental biological data (Reyes and Fetz, 1993; Khalsa et al., 2003; Lahav et al., 2004; Stricker et al., 2008). This method facilitates the classification of structurally related but dynamically differentiated biochemical oscillators. The phase response functions have been classified into two types, commonly referred to as Type I and Type II (Hansel et al., 1995; Hoppensteadt and Izhikevich, 1997; Kuramoto, 2003). A Type I phase response function generally attains a positive value during an oscillation period, whereas a Type II phase response function possesses a significantly large region of negative values. A small perturbation to an oscillator advances its phase when it is in a phase that generates a positive phase response, but retards its phase when in a negative phase response. It is known that Type I and Type II phase response functions correspond to saddle-node and Hopf bifurcations, respectively (Rinzel and Ermentrout, 1998).

In this study we have used three models and the phase reduction method to investigate the type of phase response function and phase difference of two weakly coupled oscillators in the steady state.

2. Results

2.1. Single Oscillator: A Phosphorylation-Dephosphorylation Cycle

Several modes of negative feedback regulation have been identified in signal transduction pathways, and these are potential candidates for emerging oscillatory phenomena. Many fundamental negative feedback models that cause the emergence of oscillations have been proposed (Kholodenko, 2006; Novák and Tyson, 2008). Here, we adopt the simple phosphorylation-dephosphorylation cycle models proposed by Kholodenko (2006). While he has considered all the possible topologies of feedback regulation to phosphorylation and dephosphorylation steps in the cycle, we use two, Models A and B (the latter of which corresponds to Model C in the original paper) in our study (Figures 1A,F, Section 4.2). In these models, negative feedback is realized by inhibiting kinase (Kin) production (or its activity) in Model A and enhancing phosphatase (Phos) production (or its activity) in Model B. First, we explore the parameters exhibiting oscillation by varying Phos and Kin for Models A and B, respectively (Figure 1). The parameter regions that can induce oscillations, resultant oscillation periods, and frequencies are shown in Figures 1B,G. The long-dashed lines indicate the parameter values that we have adopted in the following analyses. Figures 1C,H show the periodic orbits of the models. The oscillation periods in the two models are clearly very different from each other.

FIGURE 1

Figure 1. Phosphorylation-dephosphorylation cycle models. (A) Schematic of Model A, in which the phosphorylated M, Mp, decreases the production of kinase or inhibits its activity. Here, Phos takes a constant value, 200. (B) Period (blue line) and frequency (green line) of oscillation as a function of Phos. (C) Oscillatory behavior of the model network. Blue and red lines correspond to Mp and Kin, respectively. (D) Periodic orbit during one period of the model. The horizontal axis represents the phase. Colors are the same as in (C). (E) Phase response functions of the model, which are obtained by solving the adjoint equation given in Equation 4 (Ermentrout, 1996). Colors are the same as in (C). (F) Schematic of Model B, in which the phosphorylated M, Mp, increases the production of phosphatase or activates its activity. Here, Kin takes a constant value, 150. (G–J) Same as in (B–E), but here they are obtained using model B. Blue and red lines in (H–J) correspond to Mp and Phos, respectively.

We next applied a phase reduction method to these models and calculated the phase response functions of Models A and B. The details of phase reduction method are found in Section 4.1, in which the state vector X(t) is given by (Mp(t), Kin(t)) and (Mp(t), Phos(t)), where Mp represents an activated and phosphorylated form of M, for Models A and B, respectively. Figures 1D,I show the periodic orbits for one oscillation period in Models A and B, respectively, in which the peaks of the Mp (blue lines) are located at the origin of the phase. The phase response functions of Mp of both models are similar to each other and have significantly large regions of both positive and negative values (Figures 1E,J for Models A and B, respectively), which means that they can be categorized as Type II oscillators. This result implies that these regulatory networks have similar phase responses of Mp to a small perturbation regardless of the difference in their biological feedback targets.

2.2. Coupled Oscillators: Interacting Phosphorylation-Dephosphorylation Cycles

Next, we investigated the behavior of the above models in the presence of a weak interaction in the steady state by presuming that two cells are located next to each other. The phase reduction method allows calculation of the fixed phase difference of a weakly coupled pair of identical oscillators in the steady state (for details, see Section 4.1). Here, we consider the case of two identical oscillators interacting through a weak diffusion coupling (see Figures 2A,C). In addition to the model equations for Models A and B (Section 4.2), we adopt the following diffusion couplings. For Model A, the interaction function in Equations (7) and (8) is given by C((Mp_i(t), Kin_i(t)), (Mp_j(t), Kin_j(t))) = (−gMp(Mp_i(t) − Mp_j(t)), −gKin(Kin_i(t) − Kin_j(t))). Similarly, for Model B, it is given by C((Mp_i(t), Phos_i(t)), (Mp_j(t), Phos_j(t))) = (−gMp(Mp_i(t) − Mp_j(t)), −gPhos(Phos_i(t) − Phos_j(t))). We assume that the coefficients of interactions, gMp, gKin, and gPhos, are sufficiently small for each oscillator to remain in the basin of the periodic orbit. Under weak coupling conditions, a coupled dynamical system can be reduced to a system governed only by phase difference. The reduced dynamical system is given by Equations (17) and (18). Once we obtain the phase response function, we can calculate the gamma function (Γ ⁻(Φ), defined by Equations (17) and (18)), thereby delineating the dynamics of the phase difference. The gamma functions for Models A and B are shown in Figures 2B,D, respectively. As shown in the figures, Φ = 0 and 0.5 satisfy Γ ⁻(Φ) = 0 for both models, and only Φ = 0 satisfies the stability conditions (Equations (19) and (20)). Therefore, the oscillators in both coupled dynamical systems are asymptotically synchronized in the steady state. Next, we evaluated the shape of the gamma function by varying the ratio gKin/gMp (gPhos/gMp) between 0 and 1 for Model A (Model B) and find that Φ = 0 is the only stable solution in this parameter region. Therefore, these results suggest that in a many body system consisting of each model, their oscillation cycles can synchronize with an almost zero phase difference in a noisy environment.

FIGURE 2

Figure 2. (A) Schematic of a pair of the identical oscillators of Model A. Here, we assume weak interactions between each Mp and Kin. (B) Gamma function of the coupled system shown in (A), in which we adopt gMp = 1 and gKin = 1. (C) Same as (A), but Model B was used by assuming weak interactions between each Mp and Kin. (D) Gamma function of the coupled system shown in (C), in which we adopted gMp = 1 and gPhos = 1.

2.3. MAPK Cascade Model

Next, we considered the asynchronous oscillations experimentally observed by Shankaran et al. (2009). They reported that EGF induces oscillations in the nuclear localization (an indication of activation) of ERK in living cells and that the oscillations are asynchronous between neighboring cells (Figure 1 in Shankaran et al., 2009). Here, we adopt the MAPK cascade model that was originally developed by Huang and Ferrell (1996) and later modified by Qiao et al. (2007). The Huang-Ferrell model has been widely used for analyzing the dynamic behavior of the MAPK system (for example, Ferrell and Machleder, 1998). The original model and parameters produce an ultrasensitive MAPK activity with a Hill coefficient of 4-5, but later, Qiao et al. found that the model can also produce oscillations with other parameter sets in a wide range of the parameter space. Here, we used the Huang-Ferrell model modified by Qiao et al. to analyze the oscillatory behavior of MAPK activity (Figure 3A, Section 4.3).

FIGURE 3

Figure 3. (A) Schematic of the Huang-Ferrell model. Here, MAPKKK* indicates the activation form of MAPKKK. (B) Oscillatory behaviors of MAPKPP and MAPKPP in the steady state of the model. (C,D) Phase response functions of MAPKPP (C) and MAPKPP (D), respectively. (E) Gamma function without a time delay (black line) and one with a time delay τ = 0.25 (red line). (F) Stable (blue) and unstable (red) fixed phase differences as a function of time delay τ.

The time evolutions of MAPKKPP and MAPKPP in an oscillatory state are shown in Figure 3B. Here, we assume that MAPKPP is an activated ERK as observed by Shankaran et al. We applied the phase reduction method to this model and evaluated the stable fixed phase difference by calculating the phase response function of the model. The phase response functions of MAPKKPP and MAPKPP are shown in Figures 3C,D, respectively, in which the origin of the phase corresponds to the peak of MAPKPP. As shown in the Figure 3D, MAPKPP showed a Type II phase response, the same as those of Models A and B described in the previous section.

Next, we considered the phase difference of a pair of identical MAPK oscillators. Because the entire signal transduction network in the cells used by Shankaran et al. remains unknown, we substituted a MAPK cascade as an EGF-induced signaling pathway in the observed neighboring cells. In addition, we assumed that the interaction between a pair of cells can be effectively modeled as a function of the difference between the MAPKPPs. Then we adopted a simple diffusion type interaction function: – gMAPKPP(MAPKPP_i(t) – MAPKPP_j(t)). Using this function, we calculated the gamma function as shown in Figure 3E (black line). The value Φ = 0 is the only stable solution and the two MAPK oscillators are synchronized asymptotically. Thus, employment of only a diffusion interaction between the two MAPK models (or cells) was insufficient to reproduce the asynchronous oscillation of MAPKPP that has been detected in living cells.

To achieve an asynchronous oscillation, the sign of the coefficient can be changed. When doing so, we obtained a gamma function whose shape is the reflection of the one shown as a black line in Figure 3E and Φ = 0.5 results in a stable solution. Because we considered an effective interaction in our model, a negative diffusion coefficient can be considered. Another possible scenario is to incorporate a time delay within the interaction, – gMAPKPP(MAPKPP_i(t) – MAPKPP_j(t − τ)), where τ > 0. This scenario is reasonable because spatially isolated cells need to communicate with each other, but the signal from one cell will be delivered with a time delay to another cell when the two cells remain apart. In agreement with this theoretical assumption, oscillatory ERK activity was only observed in cells grown at low density but not in cells at confluency (Shankaran et al., 2009). When we incorporated a time delay into the interaction, the gamma function could be very simply calculated (Equations (30) and (31)). Considering τ = 0.25 as a delay time, we calculated the gamma function in Equation (31), as illustrated in Figure 3E (red line). As a result, Φ = 0.5 transforms into a stable solution. Therefore, ERK signals potentially asynchronously oscillate in the delayed system. In Figure 3F, we show how stable and unstable fixed phase solutions alter by varying the time delay, τ, between 0 and 1. A time delay 0.24 ≤ τ ≤ 0.63 results in Φ = 0.5, which is the only stable fixed phase difference. Interestingly, a wide range of time delays could induce asynchronous oscillations of MAPK.

3. Discussion

In this study, we first evaluated the phase response function, a characteristic property related to oscillation, for two simple phosphorylation-dephosphorylation cycle models A and B, in which negative feedback regulation either inhibits the kinase activity or enhances the phosphatase activity. The two models have relatively different oscillatory periods. However, their phase response functions corresponding to the reaction product Mp are very similar and both have a Type II phase response, i.e., their phase can be retarded or advanced depending on the timing of an external stimulus. This result suggests that even if cells use different negative feedback regulatory mechanisms (e.g., kinase inhibition or phosphatase activation), they can produce similar, although not exactly the same, oscillatory dynamics. It also indicates that similar oscillatory signals can be generated using different biological components and resources, such as kinases and phosphatases, produced by regulated transcription or activated by post-transcriptional modification in a cell context manner. Furthermore, we investigated the fixed phase differences of a pair of identical oscillators interacting through weak diffusion coupling. Here we found that the oscillations originating from each model can synchronize even if their interaction is weak. This result implies that when two almost identical oscillations (i.e., the biochemical reaction dynamics generated from Model A or B) can interact in spatially separated subcellular locations (such as the cytoplasm and nucleus), their oscillatory signals can still be synchronized. These results suggest that oscillation dynamics can be robustly maintained within the cell.

We next analyzed the MAPK cascade and found that the phase response function corresponding to MAPKPP is also Type II. When two identical MAPK oscillators are coupled through diffusion, these signals are synchronized. Again, in addition to the above two models, this result confirmed that oscillatory signals can be robust within the cell. However, the presumed interaction of MAPKs mediated by diffusion in our analysis did not satisfactorily reproduce the asynchronous ERK oscillation dynamics reported by Shankaran et al. (2009). We reasoned that this inconsistency could be the result of a time delay in the coupling and, under this condition, asynchronous oscillation was observed with a wide range of parameters in time delays. Our results suggest that whether a pair of cells oscillates synchronously or asynchronously may depend on the distance between the cells. In the experimental setting of Shankaran et al., oscillatory ERK activity was only observed in cells grown at low density but not in confluent cells. Therefore, we presume that compensatory mechanisms would mask the oscillatory dynamics of ERK in cells at high density. In reality, thousands of MAPK molecules exist in a cell, therefore molecular regulation of MAPK in a cell in a tightly packed cell population may be quite different from that in a sparse cell population, and therefore the former circumstance would interfere with experimental visualization of oscillatory dynamics of MAPK.

In this study, we only dealt with identical pathways under weak interaction to examine the effect of two units coupling for an oscillatory response. However, our current study suggests that a pair of slightly different oscillators under weak interaction conditions can result in similar synchronous behavior: therefore the method should be applicable for an evaluation of non-identical oscillatory units as well. The phase reduction method thus can be applied to pairs of widely different oscillator modules, as are frequently observed in biological networks.

4. Materials and Methods

4.1. Phase Reduction Method

The phase reduction method has been widely applied to coupled oscillator systems. Here, we briefly describe the derivation of the phase model. First, we define the phase. We consider a system whose dynamics is described in the vector form as follows:

$\frac{d\mathbf{\text{X}}\left(t\right)}{dt}=\mathbf{\text{F}}\left(\mathbf{\text{X}}\left(t\right)\right),$ (1)

where X(t) represents the dynamical variable, and F(X(t)) is a vector function that determines the dynamics of the system. Hereafter, we restrict the phase space of the system to a region in its basin of attraction in which a stable limit cycle, χ(t), exists. We assume that the period of the solution is T, i.e., χ(t + T) = χ(t). According to the phase reduction method, such a system can be reduced to a system consisting of the phase degree of freedom only. We can define a phase Φ(t) ∈ [0, 1) along the periodic orbit χ(t) with constant time derivative as follows:

$\frac{d\varphi \left(t\right)}{dt}=1,$ (2)

in which we can define the origin arbitrarily. Because the periodic orbit can be parameterized by Φ(t), we describe it by χ(Φ(t)). Although the phase is defined along the periodic orbit, it can be extended into the neighborhood of the periodic orbit, which is known as the asymptotic phase. Because the asymptotic phase can be described as a function of the state X(t), we can derive its dynamical equation as follows:

$\begin{array}{ll}\frac{d\varphi \left(\mathbf{\text{X}}\left(t\right)\right)}{dt}& =\frac{\partial \varphi \left(\mathbf{\text{X}}\right)}{\partial X}\cdot \frac{d\mathbf{\text{X}}\left(t\right)}{dt}\\ & =\frac{\partial \varphi \left(\mathbf{\text{X}}\right)}{\partial \mathbf{\text{X}}}\cdot \mathbf{\text{F}}\left(\mathbf{\text{X}}\left(t\right)\right),\end{array}$ (3)

where · denotes the inner product. From Equations (2) and (3), Φ(X(t)) must satisfy the following expression

$\frac{\partial \varphi \left(\mathbf{\text{X}}\right)}{\partial \mathbf{\text{X}}}\cdot \mathbf{\text{F}}\left(\mathbf{\text{X}}\right)=1.$ (4)

Next, we derive the phase response of the system. We assume that the system Equation (1) receives a weak perturbation:

$\frac{d\mathbf{\text{X}}\left(t\right)}{dt}=\mathbf{\text{F}}\left(\mathbf{\text{X}}\left(t\right)\right)+\mathbf{\text{p}}\left(t\right).$ (5)

Similar to Equation (3), the asymptotic phase obeys the following dynamical equation:

$\begin{array}{ll}\frac{d\varphi \left(\mathbf{\text{X}}\left(t\right)\right)}{dt}& =\frac{\partial \varphi \left(\mathbf{\text{X}}\right)}{\partial \mathbf{\text{X}}}\cdot \frac{d\mathbf{\text{X}}\left(t\right)}{dt}\\ & =\frac{\partial \varphi \left(\mathbf{\text{X}}\right)}{\partial \mathbf{\text{X}}}\cdot \left(\mathbf{\text{F}}\left(\mathbf{\text{X}}\left(t\right)\right)+\mathbf{\text{p}}\left(t\right)\right)\\ & =1+\frac{\partial \varphi \left(\mathbf{\text{X}}\right)}{\partial \mathbf{\text{X}}}\cdot \mathbf{\text{p}}\left(t\right)\\ & =1+\mathbf{\text{Z}}\left(\varphi \right)\cdot \mathbf{\text{p}}\left(t\right),\end{array}$ (6)

in which Z(Φ) is called as the “phase response function” or “phase sensitivity.” The phase response function describes the phase response (phase shift) of the system to small perturbations.

Now, we consider a coupled oscillator system consisting of two identical oscillatory units,

$\begin{array}{lll}& \frac{d{\mathbf{\text{X}}}_1\left(t\right)}{dt}=\mathbf{\text{F}}\left({\mathbf{\text{X}}}_1\left(t\right)\right)+\mathbf{\text{C}}\left({\mathbf{\text{X}}}_1\left(t\right),{\mathbf{\text{X}}}_2\left(t\right)\right),& \end{array}$ (7)

$\begin{array}{lll}& \frac{d{\mathbf{\text{X}}}_2\left(t\right)}{dt}=\mathbf{\text{F}}\left({\mathbf{\text{X}}}_2\left(t\right)\right)+\mathbf{\text{C}}\left({\mathbf{\text{X}}}_2\left(t\right),{\mathbf{\text{X}}}_1\left(t\right)\right).& \end{array}$ (8)

Here, X _i(t) represents the dynamical variables corresponding to the i-th unit, and C is the interaction function. When the interaction is weak (“weak” means that the units do not leave the basin of attraction of the periodic orbit), the coupled system can be described by the asymptotic phase as follows:

$\begin{array}{lll}\frac{d{\varphi }_1\left(t\right)}{dt}& =1+\mathbf{\text{Z}}\left({\varphi }_1\left(t\right)\right)\cdot \mathbf{\text{C}}\left({\varphi }_1\left(t\right),{\varphi }_2\left(t\right)\right),& \end{array}$ (9)

$\begin{array}{lll}\frac{d{\varphi }_2\left(t\right)}{dt}& =1+\mathbf{\text{Z}}\left({\varphi }_2\left(t\right)\right)\cdot \mathbf{\text{C}}\left({\varphi }_2\left(t\right),{\varphi }_1\left(t\right)\right).& \end{array}$ (10)

By replacing the phase Φ_1,2(t) with the phase ψ_1,2(t) defined by Φ_1,2(t) = t + ψ_1,2(t), Equations (9) and (10) are transformed into the following equations:

$\begin{array}{lll}& \frac{d{\psi }_1\left(t\right)}{dt}=\Gamma \left({\psi }_1\left(t\right)-{\psi }_2\left(t\right)\right),& \end{array}$ (11)

$\begin{array}{lll}& \frac{d{\psi }_2\left(t\right)}{dt}=\Gamma \left({\psi }_2\left(t\right)-{\psi }_1\left(t\right)\right),& \end{array}$ (12)

where

$\Gamma \left({\psi }_i\left(t\right)-{\psi }_j\left(t\right)\right)={\int }_0^{1}d\theta \mathbf{\text{Z}}\left(\theta +{\psi }_i\left(t\right)\right)\cdot \mathbf{\text{C}}\left(\theta +{\psi }_i\left(t\right),\theta +{\psi }_j\left(t\right)\right),$ (13)

in which we used the fact that ψ_1,2(t) vary slowly, the weak interaction assumption, and the average during one period. By substituting ψ_1,2(t) = −t + Φ_1,2(t) into Equations (11–13), we obtain the phase description of the coupled oscillator system as follows:

$\begin{array}{lll}\frac{d{\varphi }_1\left(t\right)}{dt}& =1+\Gamma \left({\varphi }_1\left(t\right),{\varphi }_2\left(t\right)\right),& \end{array}$ (14)

$\begin{array}{lll}\frac{d{\varphi }_2\left(t\right)}{dt}& =1+\Gamma \left({\varphi }_2\left(t\right),{\varphi }_1\left(t\right)\right),& \end{array}$ (15)

$\begin{array}{lll}\Gamma \left({\varphi }_i\left(t\right),{\varphi }_j\left(t\right)\right)& ={\int }_0^{1}d\theta \mathbf{\text{Z}}\left(\theta +{\varphi }_i\left(t\right)\right)\cdot \mathbf{\text{C}}\left(\theta +{\varphi }_i\left(t\right),\theta +{\varphi }_j\left(t\right)\right).& \end{array}$ (16)

By subtracting Equation (15) from Equation (14), the ordinary differential equation (ODE) is obtained that governs the evolution of the phase difference Φ(t) = Φ₁(t) − Φ₂(t),

$\begin{array}{lll}\frac{d\varphi \left(t\right)}{dt}& ={\Gamma }^{-}\left(\varphi \left(t\right)\right)\equiv \Gamma \left(\varphi \left(t\right)\right)-\Gamma \left(-\varphi \left(t\right)\right),& \end{array}$ (17)

$\begin{array}{lll}\Gamma \left(\varphi \left(t\right)\right)& ={\int }_0^{1}d\theta \mathbf{\text{Z}}\left(\theta +\varphi \left(t\right)\right)\cdot \mathbf{\text{C}}\left(\theta +\varphi \left(t\right),\theta \right).& \end{array}$ (18)

The fixed phase difference Φ is stable if it satisfies the following conditions:

$\begin{array}{lll}{\Gamma }^{-}\left(\varphi \right)& =0,& \end{array}$ (19)

$\begin{array}{lll}\frac{d{\Gamma }^{-}\left(\varphi \right)}{d\varphi }& <0,& \end{array}$ (20)

and it is unstable if it satisfies the following conditions:

$\begin{array}{lll}{\Gamma }^{-}\left(\varphi \right)& =0,& \end{array}$ (21)

$\begin{array}{lll}\frac{d{\Gamma }^{-}\left(\varphi \right)}{d\varphi }& >0.& \end{array}$ (22)

Finally, we derive a phase model corresponding to a coupled oscillator system with a delayed interaction. In this case, the system is expressed as follows:

$\begin{array}{lll}& \frac{d{\mathbf{\text{X}}}_1\left(t\right)}{dt}=\mathbf{\text{F}}\left({\mathbf{\text{X}}}_1\left(t\right)\right)+\mathbf{\text{C}}\left({\mathbf{\text{X}}}_1\left(t\right),{\mathbf{\text{X}}}_2\left(t-\tau \right)\right),& \end{array}$ (23)

$\begin{array}{lll}& \frac{d{\mathbf{\text{X}}}_2\left(t\right)}{dt}=\mathbf{\text{F}}\left({\mathbf{\text{X}}}_2\left(t\right)\right)+\mathbf{\text{C}}\left({\mathbf{\text{X}}}_2\left(t\right),{\mathbf{\text{X}}}_1\left(t-\tau \right)\right).& \end{array}$ (24)

Then, Equations (9) and (10) are changed into the following equations:

$\begin{array}{lll}& \frac{d{\varphi }_1\left(t\right)}{dt}=1+\mathbf{\text{Z}}\left({\varphi }_1\left(t\right)\right)\cdot \mathbf{\text{C}}\left({\varphi }_1\left(t\right),{\varphi }_2\left(t-\tau \right)\right),& \end{array}$ (25)

$\begin{array}{lll}& \frac{d{\varphi }_2\left(t\right)}{dt}=1+\mathbf{\text{Z}}\left({\varphi }_2\left(t\right)\right)\cdot \mathbf{\text{C}}\left({\varphi }_2\left(t\right),{\varphi }_1\left(t-\tau \right)\right).& \end{array}$ (26)

In a first order approximation, we can assume that Φ_1,2(t − τ) = Φ_1,2(t) − τ. Then, the phase description of the system is given by the following equations:

$\begin{array}{lll}\frac{d{\varphi }_1\left(t\right)}{dt}& =1+\Gamma \left({\varphi }_1\left(t\right),{\varphi }_2\left(t\right)\right),& \end{array}$ (27)

$\begin{array}{lll}\frac{d{\varphi }_2\left(t\right)}{dt}& =1+\Gamma \left({\varphi }_2\left(t\right),{\varphi }_1\left(t\right)\right),& \end{array}$ (28)

$\begin{array}{llll}\Gamma \left({\varphi }_i\left(t\right),{\varphi }_j\left(t\right)\right)& ={\int }_0^{1}d\theta \mathbf{\text{Z}}\left(\theta +{\varphi }_i\left(t\right)\right)\cdot \mathbf{\text{C}}(\theta +{\varphi }_i\left(t\right),\theta & & \\ & +{\varphi }_j\left(t\right)-\tau ).& \end{array}$ (29)

Therefore, the ODE of the phase difference attains the following form:

$\begin{array}{lll}\frac{d\varphi \left(t\right)}{dt}& ={\Gamma }^{-}\left(\varphi \left(t\right)\right)\equiv \Gamma \left(\varphi \left(t\right)\right)-\Gamma \left(-\varphi \left(t\right)\right),& \end{array}$ (30)

$\begin{array}{lll}\Gamma \left(\varphi \left(t\right)\right)& ={\int }_0^{1}d\theta \mathbf{\text{Z}}\left(\theta +\varphi \left(t\right)\right)\cdot \mathbf{\text{C}}\left(\theta +\varphi \left(t\right),\theta -\tau \right),& \end{array}$ (31)

in which the delay is no longer explicitly included in the equation.

4.2. Simple Phosphorylation Model

Model A

$\begin{array}{lll}\frac{dMp}{dt}& =v_{kin}-v_{phos}& \end{array}$ (32)

$\begin{array}{lll}& =\frac{k_{kin}^{cat}KinM}{K_{m1}+M}\cdot \frac{1+AMp∕K_a}{1+Mp∕K_a}-\frac{k_{phos}^{cat}PhosMp}{K_{m2}+Mp},& \end{array}$ (33)

$\begin{array}{lll}\frac{dKin}{dt}& =v_{kin}^{synth}-v_{kin}^{deg}& \end{array}$ (34)

$\begin{array}{lll}& =V_{kin}^0\frac{1+Mp∕K_I}{1+I\cdot Mp∕K_I}-k_{kin}^{deg}Kin,& \end{array}$ (35)

$\begin{array}{lll}M& =M^{tot}-Mp.& \end{array}$ (36)

Model B

$\begin{array}{lll}\frac{dMp}{dt}& =v_{kin}-v_{phos}& \end{array}$ (37)

$\begin{array}{lll}& =\frac{k_{kin}^{cat}KinM}{K_{m1}+M}\cdot \frac{1+AMp∕K_a}{1+Mp∕K_a}-\frac{k_{phos}^{cat}PhosMp}{K_{m2}+Mp},& \end{array}$ (38)

$\begin{array}{lll}\frac{dPhos}{dt}& =v_{phos}^{synth}-v_{phos}^{deg}& \end{array}$ (39)

$\begin{array}{lll}& =V_{phos}^0\frac{1+A_{dp}Mp∕K_d}{1+Mp∕K_d}-k_{phos}^{deg}Phos,& \end{array}$ (40)

$\begin{array}{lll}M& =M^{tot}-Mp.& \end{array}$ (41)

The parameters of Models A and B are given below.

4.3. Huang-Ferrell Model

The ODEs of the model are as follows:

$\begin{array}{lll}\frac{d{v}_0}{dt}& =a_0{c}_0{c}_1-\left(d_0+k_0\right)v_0,& \end{array}$ (42)

$\begin{array}{llll}\frac{d{v}_1}{dt}& =k_0{v}_0-a_1{v}_1{c}_2+d_1{v}_2-a_2{v}_1{c}_3+\left(k_2+d_2\right)v_3& & \\ & -a_4{v}_4{v}_1+\left(k_4+d_4\right)v_6,& \end{array}$ (43)

$\begin{array}{lll}\frac{d{v}_2}{dt}& =a_1{v}_1{c}_2-\left(d_1+k_1\right)v_2,& \end{array}$ (44)

$\begin{array}{lll}\frac{d{v}_3}{dt}& =a_2{c}_3{v}_1-\left(d_2+k_2\right)v_3,& \end{array}$ (45)

$\begin{array}{lll}\frac{d{v}_4}{dt}& =k_2{v}_3-a_3{v}_4{c}_4+d_3{v}_5-a_4{v}_4{v}_1+d_4{v}_6+k_5{v}_8,& \end{array}$ (46)

$\begin{array}{lll}\frac{d{v}_5}{dt}& =a_3{v}_4{c}_4-\left(d_3+k_3\right)v_5,& \end{array}$ (47)

$\begin{array}{lll}\frac{d{v}_6}{dt}& =a_4{v}_4{v}_1-\left(d_4+k_4\right)v_6,& \end{array}$ (48)

$\begin{array}{llll}\frac{d{v}_7}{dt}& =k_4{v}_6-a_5{v}_7{c}_4+d_5{v}_8-a_6{v}_7{c}_5+\left(d_6+k_6\right)v_9& & \\ & -a_8{v}_{10}{v}_7+\left(d_8+k_8\right)v_{12},& \end{array}$ (49)

$\begin{array}{lll}\frac{d{v}_8}{dt}& =a_5{v}_7{c}_4-\left(d_5+k_5\right)v_8,& \end{array}$ (50)

$\begin{array}{lll}\frac{d{v}_9}{dt}& =a_6{c}_5{v}_7-\left(d_6+k_6\right)v_9,& \end{array}$ (51)

$\begin{array}{lll}\frac{d{v}_{10}}{dt}& =k_6{v}_9-a_7{v}_{10}{c}_6+d_7{v}_{11}-a_8{v}_{10}{v}_7+d_8{v}_{12}+k_9{v}_{14},& \end{array}$ (52)

$\begin{array}{lll}\frac{d{v}_{11}}{dt}& =a_7{v}_{10}{c}_6-\left(d_7+k_7\right)v_{11},& \end{array}$ (53)

$\begin{array}{lll}\frac{d{v}_{12}}{dt}& =a_8{v}_{10}{v}_7-\left(d_8+k_8\right)v_{12},& \end{array}$ (54)

$\begin{array}{lll}\frac{d{v}_{13}}{dt}& =k_8{v}_{12}-a_9{v}_{13}{c}_6+d_9{v}_{14},& \end{array}$ (55)

$\begin{array}{lll}\frac{d{v}_{14}}{dt}& =a_9{v}_{13}{c}_6-\left(d_9+k_9\right)v_{14}.& \end{array}$ (56)

The algebraic equations and total quantities of the molecules are as follows:

$\begin{array}{lll}{c}_0& =\text{KK}{\text{K}}_{\text{tot}}-\left(v_0+v_1+v_2+v_3+v_6\right),& \end{array}$ (57)

$\begin{array}{lll}{c}_1& =\text{E}{\text{1}}_{\text{tot}}-v_0,& \end{array}$ (58)

$\begin{array}{lll}{c}_2& =\text{E}{\text{2}}_{\text{tot}}-v_5,& \end{array}$ (59)

$\begin{array}{llll}{c}_3& =\text{K}{\text{K}}_{\text{tot}}-(v_3+v_4+v_5+v_6+v_7& & \\ & +v_8+v_9+v_{12}),& \end{array}$ (60)

$\begin{array}{lll}{c}_4& ={\text{KKP}}^{\prime }\text{as}{\text{e}}_{\text{tot}}-\left(v_5+v_8\right),& \end{array}$ (61)

$\begin{array}{lll}{c}_5& ={\text{K}}_{\text{tot}}-\left(v_9+v_{10}+v_{11}+v_{12}+v_{13}+v_{14}\right),& \end{array}$ (62)

$\begin{array}{lll}{c}_6& ={\text{KP}}^{\prime }\text{as}{\text{e}}_{\text{tot}}-\left(v_{11}+v_{14}\right),& \end{array}$ (63)

$\begin{array}{lll}\text{KK}{\text{K}}_{\text{tot}}& =3.7112E-003,& \end{array}$ (64)

$\begin{array}{lll}\text{E}{\text{1}}_{\text{tot}}& =1.0000E-007,& \end{array}$ (65)

$\begin{array}{lll}\text{E}{\text{2}}_{\text{tot}}& =1.2115E-004,& \end{array}$ (66)

$\begin{array}{lll}\text{K}{\text{K}}_{\text{tot}}& =4.0658,& \end{array}$ (67)

$\begin{array}{lll}{\text{KKP}}^{\prime }\text{as}{\text{e}}_{\text{tot}}& =9.4579E-005,& \end{array}$ (68)

$\begin{array}{lll}{\text{K}}_{\text{tot}}& =2.0,& \end{array}$ (69)

$\begin{array}{lll}{\text{KP}}^{\prime }\text{as}{\text{e}}_{\text{tot}}& =7.5372E-002.& \end{array}$ (70)

The kinetic parameters are given below:

Molecular names, variables, and initial conditions.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

This work was supported by JSPS KAKENHI Grant (Number 22700321) to Masaki Nomura and in part by Cell Innovation Project (MEXT, Japan) to Mariko Okada-Hatakeyama.

References

Ermentrout, B. (1996). Type I x. Neural Comput. 8, 979–1001.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Guckenheimer, J., and Holmes, P. (1983). Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields. New York: Springer.

Hale, J. K., and Koçak, H. (1991). Dynamics and Bifurcations. New York: Springer-Verlag.

Hansel, D., Mato, G., and Meunier, C. (1995). Synchrony in excitatory neural networks. Neural Comput. 7, 307–337.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Hoppensteadt, F. C., and Izhikevich, E. M. (1997). Weakly Connected Neural Networks. New York: Springer-Verlag.

Huang, C.-Y. F., and Ferrell, J. E. Jr. (1996). Ultrasensitivity in the mitogen-activated protein kinase cascade. Proc. Natl. Acad. Sci. U.S.A. 93, 10078–10083.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Ferrell, J. E. Jr., and Machleder, E. M. (1998). The biochemical basis of an all-or-none cell fate switch in Xenopus oocytes. Science 280, 895–898.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Khalsa, S. B. S., Jewett, M. E., Cajochen, C., and Czeisler, C. A. (2003). A phase response curve to single bright light pulses in human subjects. J. Physiol. 549, 945–952.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Kholodenko, B. N. (2006). Cell-signalling dynamics in time and space. Nat. Rev. Mol. Cell Biol. 7, 165–176.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Kuramoto, Y. (2003). Chemical Oscillations, Waves, and Turbulence. New York: Dover Publications.

Lahav, G., Rosenfeld, N., Sigal, A., Geva-Zatorsky, N., Levine, A. J., Elowitz, M. B., et al. (2004). Dynamics of the p53-Mdm2 feedback loop in individual cells. Nat. Genet. 36, 147–150.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Novák, B., and Tyson, J. J. (2008). Design principles of biochemical oscillators. Nat. Rev. Mol. Cell Biol. 9, 981–991.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Qiao, L., Nachbar, R. B., Kevrekidis, I. G., and Shvartsman, S. Y. (2007). Bistability and oscillations in the Huang-Ferrell model of MAPK signaling. PLoS Comput. Biol. 3:e184. doi:10.1371/journal.pcbi.0030184

CrossRef Full Text

Reyes, A. D., and Fetz, E. E. (1993). Effects of transient depolarizing potentials on the firing rate of cat neocortical neurons. J. Neurophysiol. 69, 1673–1683.

Pubmed Abstract | Pubmed Full Text

Rinzel, J., and Ermentrout, B. (1998). “Analysis of neural excitability and oscillations,” in Methods in Neural Modeling, eds C. Koch and I. Segev (Cambridge, MA: The MIT Press), 251–291.

Shankaran, H., Ippolito, D. L., Chrisler, W. B., Resat, H., Bollinger, N., Opresko, L. K., et al. (2009). Rapid and sustained nuclear-cytoplasmic ERK oscillations induced by epidermal growth factor. Mol. Syst. Biol. 5, 332.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Stricker, J., Cookson, S., Bennett, M. R., Mather, W. H., Tsimring, L. S., and Hasty, J. (2008). A fast, robust and tunable synthetic gene oscillator. Nature 456, 516–519.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Strogatz, S. H. (2001). Nonlinear Dynamics and Chaos. Boulder, CO: Westview Press.

Strogatz, S. H. (2003). Sync: The Emerging Science of Spontaneous Order. New York: Hyperion Books.

Winfree, A. T. (2001). The Geometry of Biological Time. New York: Springer.