Abstract
Biological oscillators coordinate individual cellular components so that they function coherently and collectively. They are typically composed of multiple feedback loops, and period mismatch is unavoidable in biological implementations. We investigated the advantageous effect of this period mismatch in terms of a synchronization response to external stimuli. Specifically, we considered two fundamental models of genetic circuits: smooth and relaxation oscillators. Using phase reduction and Floquet multipliers, we numerically analysed their entrainability under different coupling strengths and period ratios. We found that a period mismatch induces better entrainment in both types of oscillator; the enhancement occurs in the vicinity of the bifurcation on their limit cycles. In the smooth oscillator, the optimal period ratio for the enhancement coincides with the experimentally observed ratio, which suggests biological exploitation of the period mismatch. Although the origin of multiple feedback loops is often explained as a passive mechanism to ensure robustness against perturbation, we study the active benefits of the period mismatch, which include increasing the efficiency of the genetic oscillators. Our findings show a qualitatively different perspective for both the inherent advantages of multiple loops and their essentiality.
1. Introduction
Biochemical oscillators exist primarily at the genetic level, such as in the interplay among mRNAs and proteins [1]. At least one negative feedback loop (NFL) is essential for limitcycle oscillation, and in the minimum configuration, one NFL with timedelay is sufficient to generate oscillation. However, many realistic oscillators possess extra components, such as positive feedback loops (PFL) and more than one NFL, which are apparently redundant. Since more regulatory components consume more resources (e.g. amino acids for synthesis and adenosine triphosphate for operation), this would seem energetically disadvantageous, and thus not preferable through evolution. In addition, the coexistence of positive and negative loops may make the resulting behaviour incoherent. Still, the circadian clock of even the simplest life form (e.g. the marine dinoflagellate Gonyaulax polyedra) is composed of multiple oscillators whose rhythms exhibit different periodic patterns [2,3]. Usually this observation is explained by the notion of biological robustness, in which there is a backup mechanism or a resistance to perturbation. However, such observations do not suggest the essentiality of multiple loops: why does not even a single species possess the minimum configuration? In this study, we numerically analyse the inherent advantages of multiple loops from the perspective of entrainment and the adaptation of the system to external input [4].
Biologists have investigated coupled oscillators at the molecular level in many different organisms. BellPedersen et al. [5] reviewed circadian clocks in several species from different clades. They suggested that multiple loops comprise pacemaker and slave oscillators in unicellular organisms, whereas in mammals and avians, a centralized pacemaker entrains downstream systems. Circadian systems have also been investigated theoretically. Wagner et al. studied the stability of an oscillator (called the Goodwin model [6]) against perturbation of the kinetic rate, and found that interlinked loops are more robust (by nearly one order of magnitude) than a nonlinked counterpart [7]. Other benefits of multiple feedback loops are found in [8,9] and their intuitive advantages are summarized in the review by Hastings [10].
When multiple loops are involved, synchronization becomes a crucial issue. In biological implementations, no loops can share an identical period in a strict sense. At first glance, such disorder seems to degrade the system performance. However, the results in nonlinear physics indicate that a certain amount of disorder actually enhances performance through mechanisms similar to those in stochastic resonance [11–18], Brownian transport [19–23], noiseinduced synchronization [24–26] and disorderinduced resonance [27].
Our focus is on the property called entrainability or the ability of oscillators to be synchronized with an external periodic signal. It is directly connected to the adaptability of the species and therefore its survival [28]. Not only does a Drosophila with a circadian clock survive better than one without [29], but also an Arabidopsis with better entrainment with the environment has increased photosynthesis and consequently grows faster [30]. Even cyanobacteria with better entrainability have an advantage over nonsynchronizable ones, although this advantage disappears in constant environments [31].
At the molecular level, entrainability has been attributed to several different mechanisms. Gonze et al. [32] studied a population of suprachiasmatic nuclei (SCN) neurons and showed that efficient global synchronization is achieved via damping. Liu et al. [33] experimentally confirmed that robust oscillatory behaviour can be achieved by tissuelevel communication between SCN cells, even if critical clock genes have been knocked out. Tsai et al. [34] reported that the presence of PFLs allowed for easier tuning of the period without changing the amplitude. Zhou et al. [35] showed that PFLs assist noiseinduced synchronization through their sensitivity to external signals.
Still, a fundamental question has not been answered: what is the advantage of harbouring multiple loops with mismatched periods? We investigate the effects of mismatched periods on limitcycle oscillations from the viewpoint of entrainment. Specifically, we consider the minimal fundamental model: an oscillator made of two loops, each of a different period, that are connected. By using phase reduction and Floquet multipliers, we analysed the response of this system to external periodic stimuli as a function of the coupling strength and the period ratio. Our main finding is the large effect of mismatched periods on entrainability: that is, better entrainment is achieved when the two oscillatory components have different natural periods. This enhancement effect is observed in two types of coupled oscillators: smooth and relaxation oscillators (figure 1). The enhanced entrainability seen in the present study is similar to that reported by Komin et al. [36]. However, our mechanism is different from theirs: the enhancement in our case is caused by the oscillators in the vicinity of the bifurcation on the limit cycle whereas their enhancement results from damped oscillation (oscillator death).
2. Methods
We studied the effects of mismatched periods on entrainment by using smooth and relaxation oscillators. Before explaining the details of these models, we will first explain period scaling and a coupling scheme.
2.1. The base model
An Ndimensional differential equation for a single oscillator is represented by 2.1where v and f(v) are Ndimensional column vectors defined by where ⊤ denotes a transposed vector (matrix). Equation (2.1) yields autonomous oscillation (a limitcycle oscillation) whose period is T. The period can be tuned without changing the amplitude of the oscillation by introducing a scaling parameter τ 2.2where the period of equation (2.2) is given by τ T. Other than their periods, solutions to equation (2.2) with different τ values will have the same properties [32,36]. We used equation (2.2) to study the effects of mismatched periods between two oscillatory components: 2.3and 2.4where the subscript i denotes the ith component and C_{i}(v_{i}, v_{j} ) represents a coupling term. To investigate the effects of mismatched periods, we set 2.5where R represents the ratio of the periods of the two oscillators.
2.2. Genetic oscillator models
We show two specific models of f(v), each for the smooth and relaxation oscillators.
2.2.1. Smooth oscillator
A smooth oscillator is the most basic structure for limitcycle oscillations [1], and many biochemical oscillators (e.g. the Goodwin model [6] and the repressilator [37]) belong to this class. We adopted the smooth oscillator presented in Novák & Tyson [1] (figure 1). A diagram and the simplified structure of the smooth oscillator are shown in figure 1a,b, where the transcription of X (mRNA) is inhibited by Z (protein). Protein Y plays a role in the timedelay of the loop X → Y → Z ⊣ X, where X → Y and X ⊣ Y denote activation and repression, respectively. In the presence of sufficient time delay in the NFL, the regulatory network exhibits limitcycle oscillations. The dynamics of the smooth oscillator are described by the following threedimensional differential equation on v = (x,y,z)^{⊤}, where x, y and z denote the concentrations of mRNA X, protein Y and protein Z, respectively: 2.6
Because the amplitude and the period of the limit cycle can be tuned separately by proper scaling, we employed the parameters V_{sx} = 5.44, k_{dx} = 8.99, k_{sz} = 4.49, k_{dz} = 4.49, k_{sy} = 18.0, V_{dy} = 2.72, K_{m} = 0.00303, K_{d} = 0.303, and p = 4.0 to achieve unit period and unit amplitude for y (we define the amplitude as the difference between the maximum and the minimum of the oscillation (see figure 2a)). These parameter settings are essentially identical to those presented in the electronic supplementary material of Novák & Tyson [1], except for the scale transform (V_{sx} = 0.2, k_{dx} = 0.1, k_{sz} = 0.05, k_{dz} = 0.05, k_{sy} = 0.2, V_{dy} = 0.1, K_{m} = 0.01, K_{d} = 1.0 and p = 4.0 in the original model).
We analysed two coupled smooth oscillators, as shown in figure 1c, where the subscript i of X_{i} denotes the ith oscillator component. We assumed that the protein translated from X_{1} positively regulates the transcription of X_{2}. The coupling term in equations (2.3) and (2.4) can be written in the following linear relation: 2.7where ε denotes the coupling strength.
2.2.2. Relaxation oscillator
The existence of PFLs makes oscillators excitatory, and such oscillators are known as relaxation oscillators [1]. Relaxation oscillators have cusps in their nullcline curves, and one example of this class is the FitzHugh–Nagumo oscillator [38,39].
We again employ a relaxationtype genetic oscillator presented by Novák & Tyson [1]. Figure 1d,e show specific and simplified models, respectively, of a relaxation oscillator in which transcription of X (mRNA) is inhibited by Y (protein), which serves as a NFL. Following this, Y (mRNA) is degraded by Z (enzyme) and also forms an enzyme–substrate complex ZY to inhibit the functionality of Z (). The interaction functions as the PFL in this network. In many cases, it is assumed that the concentration of enzyme Z equilibrates (ż = 0), yielding the following twodimensional differential equation in terms of , where x and y denote the concentrations of mRNA X and protein Y, respectively: 2.8
As in the smooth oscillator, we scaled the parameters so that the period and the amplitude y both became 1 (figure 2b). The parameters, V_{sx} = 0.815, K_{m} = 0.0268, k_{dx} = 3.04, k_{sy} = 60.8, k_{dy} = 3.04, V_{dy} = 16.3, K_{1} = 0.134, K_{d} = 0.268 and p = 4.0, are essentially equivalent to those presented in the electronic supplementary material of Novák & Tyson [1], except for the scale transform (V_{sx} = 0.05, K_{m} = 0.1, k_{dx} = 0.05, k_{sy} = 1.0, k_{dy} = 0.05, V_{dy} = 1.0, K_{1} = 0.5, K_{d} = 1.0, and p = 4.0 in the original model).
We analysed two coupled relaxation oscillators, as shown in figure 1f. The coupling term C_{i} was written as a linear relation: 2.9where ε denotes the coupling strength.
2.3. Analytical approaches
For an analytical analysis, we used phase reduction and Floquet multipliers. Even a system of only two coupled oscillators exhibits very complicated behaviours, such as n : m synchronization (n, m; positive integers) and chaos, depending on the model parameters (see §3). Because detailed bifurcation analysis of the coupled oscillator is outside the scope of this paper, we assumed 1 : 1 synchronization between two oscillatory components (tight coupling). Under this assumption, the coupled oscillator can be viewed as one oscillator, whose dynamics is represented by 2.10with 2.11
Incorporating an input signal in equation (2.10), we have 2.12where I(ωt) denotes an input signal with angular frequency ω and χ is a scalar representing the signal strength. We will use equation (2.12) in the calculations for entrainability and Floquet multipliers.
2.3.1. Phase reduction
A phasereduction approach uses a phase variable to reduce a multidimensional system into a onedimensional differential equation [40,41]. In the absence of external perturbations (equation (2.10)), a limitcycle oscillation can be described in terms of a phase variable ϕ ∈ [0,2π]: 2.13where Ω is the natural angular frequency of the limit cycle (Ω = 2π/T, where T is the period). Although the phase ϕ in equation (2.13) is only defined on the stable limitcycle trajectory, its definition can be expanded to the entire u space, where the equiphase surface is called the isochron ℐ(ϕ). An example of an isochron ℐ(ϕ) defined on a hypothetical limitcycle oscillator, shown by dashed lines, in figure 3a where the isochron is drawn at intervals of π/6.
In the presence of external signals (equation (2.12)), the time evolution of ϕ is described by 2.14where a dot (·) represents an inner product and U(ϕ) is the phase response curve (PRC) defined by 2.15Here, represents a point on limitcycle oscillation at the phase ϕ′. It is assumed that the perturbed trajectory is in the vicinity of the unperturbed limitcycle trajectory (i.e. χ is sufficiently small). The introduction of a slow variable ψ = ϕ− ωt in equation (2.14) yields 2.16where . Because ψ is a slow variable ( is very small), when χ is sufficiently small, the inner product term can be approximated by its average (separating the timescales of ψ and ): 17with 18(the integration of equation (2.18) is evaluated numerically in practical calculations). If stable solutions exist, this oscillator synchronizes to external signals. Θ(ψ) is a 2π periodic function, and its shape is typically given by sinusoidallike functions, as shown in figure 3b. The fixed points are the intersection points of Θ(ψ) and −ΔΩ, and only the point with is a stable solution (an empty circle in figure 3b). Therefore, the oscillator is synchronized to external signals if the following relation holds: 2.19with
The extent of synchronization in terms of χ (signal strength) and ω (signal angular frequency) is often described by a domain known as the Arnold tongue or the synchronization region [42], inside of which the oscillator can synchronize to an external periodic signal (the shaded domain in figure 3c). Although the border of the Arnold tongue is generally nonlinear as a function of χ, the upper and lower borders can be approximated linearly by equation (2.19) when χ is sufficiently small (these are shown with the dashed line in figure 3c). Thus, the range in which ω can be synchronized with an external signal can be quantified by , which represents the width of the Arnold tongue for sufficiently small χ. We define the entrainability W as a χindependent variable (see figure 3b): 2.20Writing equation (2.19) in terms of the period of the external signal, we have 2.21where T_{s} is the period of the external signal (T_{s} = 2π/ω) (note that equations (2.19) and (2.21) only hold for sufficiently small χ). The PRC equation (2.15) can be calculated by the following differential equation [41,43]: 2.22with the constraint 2.23where is a stable periodic solution of equation (2.10) (, where T is the period), is a Jacobian matrix of around , and we abbreviated (ϕ(t) is the phase as a function of t according to equation (2.13)). Because equation (2.22) is unstable, we numerically solved equation (2.22) backward in time.
2.3.2. Floquet multiplier
In the presence of external signals, the solution deviates from a stable orbit , where the deviation obeys the differential equation 2.24
Here, we assume that the deviation is sufficiently small. Because is a periodic function with period T, equation (2.24) has the Floquet solutions, generally represented by 2.25where are functions with period T, c_{i} is a coefficient determined by the initial conditions, and μ_{i} are the Floquet exponents. The Floquet multipliers ρ_{i} are related to the exponents via 2.26
Note that the Floquet exponents are a special case of the Lyapunov exponents for periodic systems (e.g. limit cycles) [44]. The Floquet multipliers can be calculated numerically from a matrix differential equation in terms of the matrix Q(t), after Klausmeier [44]: 2.27
Equation (2.27) is integrated numerically from t = 0 to t = T with an initial condition that Q(t = 0) is the identity matrix. The Floquet multipliers are obtained as the eigenvalues of Q(t = T), and hence Ndimensional limitcycle models yield N Floquet multipliers.
When the Floquet multipliers ρ_{i} cross the unit circle on the complex plane, the limit cycle bifurcates. Bifurcations of codimension 1 can occur in three patterns [45–48]:

— Saddlenode bifurcation. A single real Floquet multiplier crosses ρ = 1 (this includes transcritical or pitchfork bifurcations as special cases).

— Hopf bifurcation. Complexconjugate Floquet multipliers cross ().

— Perioddoubling bifurcation. A single real Floquet multiplier crosses ρ =−1.
In autonomous systems, one of the Floquet multipliers is ρ_{i} = 1, which corresponds to a pure phase shift. Wiesenfeld and McNamara showed that Floquet multipliers whose absolute values are close to 1 are responsible for signal amplification [47,48]. We calculated the Floquet multipliers up to the secondlargest absolute value (except for the inherent multiplier ρ = 1), which are denoted as ρ_{1} (the largest) and ρ_{2} (the second largest) in our models.
3. Results
We performed a bifurcation analysis and examined the relationship between entrainability and the Floquet multipliers. We set τ_{1} = 1 and τ_{2} = R in equations (2.3) and (2.4), where R represents the ratio of the periods of the two oscillators. When R ≠ 1, the two oscillatory components have different natural periods (without coupling) and hence the periods are mismatched. As the controllable parameters for the model calculations, we adopted R and ε, where ε denotes the coupling strength defined in equations (2.7) and (2.9).
3.1. Bifurcation diagram
Bifurcations of the limit cycle are classified into bifurcations in equilibria and those on the limit cycle. Here, bifurcation on the limit cycle was studied by numerically solving differential equations (2.6) and (2.8).
Specifically, we described the bifurcation diagrams as a function of R (period ratio) and ε (coupling strength) and calculated the local maxima of y_{1}. If the oscillator showed regular limitcycle oscillations, the maxima of y_{1} was represented by a single point. If the limit cycle underwent bifurcation, the maximal points were described by more than two points. We iterated by changing R from R_{min} to R_{max} (R_{min}< R_{max}) and then adopted the last values of the preceding R values as the initial values for the next R value. The incremental step for R was set at , where R_{min} and R_{max} are indicated as the lower and upper limits on figures 4 and 6. At each parameter value, we discarded as transient phases those during t ∈ [0,200], and we recorded the local maxima of the trajectories in during the intervals t ∈ [200, 300]. We carried out these same procedures on the coupling parameter ε (ε_{min}, ε_{max}, with the increment ).
In figure 4, the bifurcation diagram of y_{1} as a function of R in the smooth oscillator is shown for (figures 4a) ε = 1.5 and (figures 4b) ε = 2.5. Although the maxima of y_{1} slightly depend on R, there is no bifurcation. In figure 5, the bifurcation diagram as a function of ε is shown for (figure 5a) R = 1.4 and (figure 5b) R = 2.5. In both cases, bifurcations occurred at for R = 1.4 and for R = 2.5, and chaotic behaviour was observed when the coupling strength was small.
Figure 6 shows the bifurcation diagram as a function of R for the relaxation oscillator for (figure 6a) ε = 2.0 and (figure 6b) ε = 2.5. In the former case, a perioddoubling bifurcation occurred in the region 1.7 < R < 2.7. In the latter, a cusp appeared around R = 1.9, in the vicinity of a saddlenode bifurcation (see the section on entrainability and Floquet multipliers). In figure 7, the bifurcation diagram as a function of ε is shown for (figure 7a) R = 1.3 and (figure 7b) R = 1.9. A complex behaviour (chaos) appeared in a small region of ε. Especially in figure 7b, we found a branch of a maxima around 1.7 < ε < 2.3, indicating the occurrence of the perioddoubling bifurcation.
In comparison, bifurcations are less probable in a smooth oscillator than in a relaxation oscillator; bifurcations are only observed in the case of very weak coupling in smooth oscillators (figure 5). In contrast, the relaxation oscillator exhibited perioddoubling and chaotic oscillations in a wider range of coupling strengths (figures 6 and 7).
Next, we examine the time course of the coupled oscillators with mismatched periods (figures 8 and 9, respectively, for the smooth and relaxation oscillators). Without mismatched periods, the temporal trajectory of the oscillators would be identical to that of the case of a single oscillator (figure 2). Figure 8 shows the case for coupled smooth oscillators with (figure 8a) ε = 1.5, R = 1.5, and (figure 8b) ε = 2.5, R = 1.5; figure 9 shows the case for coupled relaxation oscillators with (figure 9a) ε = 2.0, R = 2.8, and (figure 9b) ε = 2.5, R = 1.9. In all the cases shown in figures 8 and 9, we see the 1 : 1 synchronization between the two oscillators. Without period mismatch, the amplitude of oscillation y is 1 (see §2). We see from figures 8 and 9 that the amplitudes of y_{1} and y_{2} depend on the period ratio R, where the amplitude of y_{1} is smaller than 1 except in figure 9b. From figures 4 and 6, which describe maxima of y_{1} as a function R, we see that y_{1} is larger than 1 in both oscillator models when the period mismatch is smaller, and the maxima of y_{1} start to decrease with increasing R. The decrease in the amplitude is larger for the smooth oscillator than for the relaxation oscillator.
3.2. Entrainability and Floquet multiplier
In using phase reduction and Floquet multipliers to calculate how the entrainability W (equation (2.20)) depends on R (equation (2.5)) and ε (equations (2.7) and (2.9)), we employed the following sinusoidal input signal: 3.1where the input signal only affects the concentration of mRNA in the first oscillator (x_{1} in equation (2.10)). We assumed 1 : 1 synchronization between the two oscillators (the regions where the maximum of y_{1} is represented by a single point in figures 4–7) so that the coupled oscillator could be viewed as a single case. Although we calculated the entrainability and the Floquet multipliers inside the 1 : 1 synchronized regions, the entrainability very close to the bifurcation points is not shown for some parameters, because equation (2.22) did not exhibit stable periodic solutions even after remaining there for a long time ( in figure 11a, in figure 12a and in figure 13d). In most cases, we calculated the Floquet multipliers up to the second most dominant ones. When the multipliers were given by complexconjugate values, we checked the three multipliers ρ_{1}, and ρ_{2}, where is the complex conjugate of ρ. When two real multipliers (ρ_{1} and ρ_{1}′) eventually emerged because of annihilation of the complexconjugate multipliers (ρ_{1} and ) by changing parameters, we checked the three multipliers ρ_{1}, ρ′_{1} and ρ_{2}.
By calculating the entrainability W and the Floquet multipliers ρ_{i}, we obtained stable limitcycle trajectories in a way similar to what we did in the bifurcation analysis. We iterated over R and adopted the values of the preceding R values as the initial values for the next R values. For each parameter value, we calculated the stable limitcycle trajectories for one period, and these trajectories were used for the calculation of the PRC functions (equation 2.15). These procedures were also carried out for the coupling parameter ε. The Floquet multipliers also require stable limitcycle trajectories, and we calculated them in the same way as described above.
3.2.1. Smooth oscillator
We first investigated the smooth oscillator as a function of R (figure 10). Specifically, figure 10a is a dualaxis plot of the entrainability W (solid line; left axis) and the period T (dashed line; right axis) for ε = 1.5. We can uniquely define the period T of the coupled oscillator because we are considering a 1 : 1 synchronization between the two oscillators. In figure 10a, the period T is not strongly affected by period mismatch and is around 1 in the whole R region, although the entrainability W depends on R. The entrainability is when the two oscillators have the same period (R = 1), whereas for . The first and the second dominant Floquet multipliers for ε = 1.5 are shown in figure 10b,c, where solid, dashed, and dotteddashed lines, respectively, denote , and (ρ_{i}: Floquet multiplier defined by equation (2.26)). In figure 10b, symmetric imaginary curves represent complexconjugate multipliers, and a branch of the real part around is owing to annihilation of the complexconjugate multipliers (two real multipliers emerge in place of the complexconjugate multipliers). We see that the ρ_{1} in figure 10b are imaginary values, and their absolute values are close to 1 for the R > 1.4 region. The coupled oscillator is close to the Hopf bifurcation points. The second dominant Floquet multiplier ρ_{2} is described in figure 10c, and we see that its magnitude is comparable to that of ρ_{1}. The coupled oscillator is near the saddlenode bifurcation points (the second Floquet multiplier ρ_{2} is a pure real value). We also carried out the same calculations with ε = 2.5 and found enhanced entrainability in the presence of a period mismatch (R > 1.5), as in figure 10c (data not shown). According to the analysis with Floquet multipliers, this enhancement can also be accounted for by the bifurcation.
Figure 11 shows the entrainability and Floquet multipliers as a function of ε in the smooth oscillator while keeping R constant. Figure 11a is a dualaxis plot of the entrainability W (left axis) and the period T (right axis) for R = 1.4, where the period T does not strongly depend on the coupling strength ε. In contrast, the entrainability W achieves a maximum, W ≃ 15 at ε ≃ 2.5, which is more than seventimes larger than the nomismatch case (W ≃ 2). Figure 11b shows the first dominant Floquet multiplier ρ_{1} as a function of ε, where the notation is the same as in figure 10. ρ_{1} is a complex multiplier for ε <5.5, and its absolute value (dotteddashed line) approaches 1 as ε decreases, showing that the oscillator is in the vicinity of the Hopf bifurcation. We see consistent results in figure 5 where the oscillator is quasiperiodic (chaotic) in the ε < 1.3 region. At ε ≃ 5.5, the complexconjugate multipliers annihilate, and two real multipliers emerge instead. Figure 11c shows the second dominant Floquet multiplier ρ_{2} as a pure positive real number, responsible for the saddlenode bifurcation points (this multiplier is close to ρ_{i} = 1). We see that ρ_{2} achieves its maximal value around ε ≃ 2.5, where the entrainability W also exhibits a maximum. Because in figure 11b monotonically decreases as a function of ε for ε < 5.5, ρ_{1} does not explain the qualitative behaviour of W, whereas ρ_{2} does. This result shows that being close to the saddlenode bifurcation is important for achieving entrainment. We investigated the entrainability with a different parameter for the period mismatch R = 2.5, and observed that W ≃ 26 was maximized at an intermediate coupling strength ε ≃ 5 (data not shown). The Floquet multiplier analysis also indicated that enhancement in the R = 2.5 case is caused by the oscillators being in the vicinity of the saddlenode bifurcation on the limit cycle.
3.2.2. Relaxation oscillator
We next apply the same analysis to the relaxation oscillator: we calculated the entrainability W, the period T and the Floquet multipliers ρ_{i}, as functions of R (figure 12) and ε (figure 13). Figure 12a shows a dualaxis plot of the entrainability W and the period T, for ε = 2.0 (see also figure 10). The coupled oscillator is not 1 : 1 synchronized in the region 1.7 < R < 2.7, and these unsynchronized regions were not plotted. Since the entrainability W grows rapidly in the vicinity of the bifurcation points (R ≃ 1.7) in figure 12a, we calculated the first ρ_{1} and second ρ_{2} dominant Floquet multipliers in figure 12b,c, where ρ_{1} is a complexconjugate multiplier around R = 1.7 and is a pure real multiplier very near the bifurcation points. We see that ρ_{1} of ε = 2.0 exhibits a sudden decrease towards −1 near R ≃ 1.7 and 2.7, indicating that the coupled oscillator is close to the perioddoubling bifurcation points. The bifurcation diagram of figure 6 shows a branch, i.e. the perioddoubling bifurcation. ρ_{1} is a real multiplier for R < 3.3, and its value approaches −1 as R decreases toward R = 2.7 in the vicinity of the perioddoubling bifurcation. At R ≃ 3.3, complexconjugate multipliers emerge instead of the annihilation of two real multipliers. Figure 12c describes the second dominant Floquet multiplier ρ_{2} of ε = 2.0, and both the imaginary and real parts of ρ_{2} vanish for a range of R values. Because the rest of the Floquet multipliers for ε = 2.0 vanish (not shown here), all the dominant Floquet multipliers of ε = 2.0 are related to the perioddoubling bifurcation points.
Figure 12d shows the entrainability W and the period T for ε = 2.5. We also see the rapid increase in W around R ≃ 1.9, and the period T increases up to T = 1.3 for 1 ≤ R < 1.9 and then decreases for R > 1.9. Figure 12e,f show the first ρ_{1} and the second ρ_{2} dominant Floquet multipliers for ε = 2.5, respectively. ρ_{1} is a pure real multiplier with a minimum around R ≃ 2.2, which does not coincide with the position of the maximum of W. Although the magnitude of ρ_{2} is smaller than ρ_{1}, it gains a peak at R ≃ 1.9 identical to the peak position of W. Thus, the system is better entrained in the vicinity of the saddlenode bifurcation, and there exists an optimal coupling strength for the entrainability.
We next show the entrainability W, the period T and the Floquet multipliers ρ_{i} as a function of ε for the relaxation oscillator. Figure 13a is a dualaxis plot of the entrainability W and the period T for R = 1.3. We see very little dependence of the entrainability or the period on ε. Figure 13b,c show the first ρ_{1} and the second ρ_{2} dominant Floquet multipliers, respectively; ρ_{1} approaches 1 with decreasing ε, indicating that the oscillator approaches the saddlenode bifurcation points, whereas ρ_{2} vanishes for a range of values of ε. Figure 13d shows the entrainability W and the period T as a function of ε for R = 1.9, where the period T does not strongly depend on the coupling strength ε. In contrast, the entrainability W exhibits nonmonotonic behaviour and achieves the local maximum around ε = 2.5. The first dominant Floquet multiplier ρ_{1} in figure 13e shows that the oscillator approaches the perioddoubling bifurcation when approaching ε = 2.3 from above. Although the first dominant Floquet multiplier shows the perioddoubling bifurcation, its magnitude does not qualitatively explain the existence of the entrainability peak around ε = 2.5. In contrast, the second dominant Floquet multiplier ρ_{2} in figure 13f shows a positivevalued peak around ε = 2.5, and the position agrees with that of the entrainability. In figure 13d, enhancement is only observed inside a narrow region where ε < 2.8.
4. Discussion
Mismatched periods enhance entrainability for both smooth and relaxation oscillators, and the enhancement is related to the dominant Floquet multipliers. The behaviour as a function of the model parameters (R and ε) can be explained qualitatively by a pure real Floquet multiplier. It has previously been shown that in the vicinity of the bifurcation points, the oscillator is more sensitive to an external signal [47,48]. Entrainability is thus enhanced because the periodmismatch drives the coupled oscillator to the vicinity of saddlenode bifurcation. The pure real multiplier has a strong impact on the entrainment property (figures 10, 11, 12d–f and 13d–f), as do the Floquet multipliers related to other bifurcation types (figure 12a–c). In figure 12a, we see that even when neither the first nor the second dominant Floquet multiplier is a positive real multiplier, enhancement is induced by the perioddoubling bifurcation points.
From our analysis above, we saw that bifurcation is less likely to occur in a smooth oscillator than in a relaxation oscillator. Nevertheless, the enhancement of entrainability (an acute peak for a relaxation oscillator and an obtuse peak for a smooth oscillator) is seen in a wider region for a smooth oscillator than for a relaxation oscillator (figures 10a and 12d for R dependence, and figures 11a and 13d for ε dependence). Although a smooth oscillator may be better entrained, Figures 8 and 9 show that the amplitude of the oscillation is more strongly affected in the smooth oscillator. When increasing the period mismatch R, the amplitude of y_{1} increases at the outset and then starts to decrease; the decrease starts sooner for the smooth oscillator than it does for the relaxation oscillator (see figures 4 and 6). Since limitcycle oscillations with smaller amplitude are more sensitive to external stimuli, the reduction of the amplitude that is induced by the period mismatch is also a cause of the improved entrainability. However, the main contributing factor to the enhancement is the bifurcation on the limit cycle. This is because the behaviour of the entrainability W agreed with that of the Floquet multiplier, which does not depend on the scaling of the values.
As mentioned in §1, many circadian oscillators consist of multioscillatory networks, each of which has different components. As an example, early works on Gonyaulax reported periodic bioluminescence (24 to 27 h: average 25 h) and aggregation (15 to 23 h: average 20 h) patterns [2,3]. These two oscillators are synchronized in the presence of white light, whereas they are decoupled in the presence of dim red light. The ratio of the two periods is about 1.25 (see fig. 2c in Roenneberg [2]), which coincides with our results of R ≃ 1.4 for enhancement in the smooth oscillator (see figure 10a).
We have shown that the oscillators are better synchronized to an external signal when they are close to bifurcation points (mainly at the saddlenode but also partially at the perioddoubling bifurcations). Regarding the perioddoubling bifurcation in realistic systems, we note that experimental results showed the existence of a circabidian rhythm, with a 48 h period (double the circadian period). It has been reported that such circabidian rhythms can occasionally emerge under controlled conditions [49]. The bidian rhythm may result from a perioddoubling bifurcation of a circadian rhythm that is close to a bifurcation point.
Regarding stochasticity, gene expression is subject to two types of noise (stochasticity): one from the discrete molecular species (intrinsic noise) and the other from environmental variability (extrinsic noise) [50–53]. We may take into account stochastic effects in order to make our findings biochemically realistic. For example, a collection of identical oscillators may exhibit damped oscillation as the averaging effect of stochasticity [54]. Nevertheless, we employed a deterministic model (i.e. ordinary differential equations) in our analysis for two reasons. First, detailed stochastic modelling, such as a continuoustime discrete Markov chain, is not tractable: the only solution is an exhaustive, timeconsuming calculation using the Monte Carlo method without any theoretical guarantees. Second, without analytical calculations, the theoretical arguments for the enhancement would be difficult. The enhancement we detected at the bifurcation on the limit cycle (the saddle node bifurcation) is unlikely to be found by numerical simulations, but can be deduced from Floquet multipliers. Although our present analyses do not incorporate stochastic aspects, we consider that our results help in the understanding of the stochastic dynamics in modelling mismatched periods. We will consider such stochasticity in a future study.
Another topic yet to be investigated is when three or more oscillators have mismatched periods. When considering more than two oscillators, besides the oscillatory model, the topology of the coupling is important. We can make an educated guess at the consequences of a period mismatch in such a multioscillator model: if period inconsistency drives the coupled system to the vicinity of bifurcation on the limit cycle, the period mismatch will enhances entrainability. As long as the overall connection of the multiple oscillators can be simplified as ‘two mismatched oscillators’, the entrainability of such a component is expected to be enhanced.
In summary, we have used phase reduction and Floquet multipliers to show that entrainability is enhanced by mismatched periods in both smooth and relaxation oscillators. We focused on two identical oscillators that were symmetrically coupled to each other by using a deterministic approach. In order to make our results more useful for real biological systems, however, it would be important to consider situations of three or more oscillators, asymmetric coupling, or structurally heterogeneous cases with stochasticity. These extensions are left to our future studies.
Acknowledgements
This work was supported by the Global COE program ‘Deciphering Biosphere from Genome Big Bang’ from the Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan (Y.H. and M.A.); a GrantinAid for Young Scientists B (no. 23700263) from MEXT, Japan (Y.H.); and a GrantinAid for Scientific Research on Innovative Areas ‘Biosynthetic Machinery’ (no. 11001359) from MEXT, Japan (M.A.).
 Received December 12, 2012.
 Accepted January 14, 2013.
 © 2013 The Author(s) Published by the Royal Society. All rights reserved.