Nonlinear phase-amplitude reduction of delay-induced oscillations
Kiyoshi Kotani, Yutaro Ogawa, Sho Shirasaka, Akihiko Akao, Yasuhiko Jimbo, Hiroya Nakao
aa r X i v : . [ n li n . AO ] J u l Nonlinear phase-amplitude reduction of delay-induced oscillations
Kiyoshi Kotani a Research Center for Advanced Science and Technology,The University of Tokyo, 4–6–1 Komaba, Meguro-ku, Tokyo 153-8904, Japanand JST, PRESTO, 4–1–8 Honcho, Kawaguchi-shi, Saitama 332-0012, Japan
Yutaro Ogawa
Research Center for Advanced Science and Technology,The University of Tokyo, 4–6–1 Komaba, Meguro-ku, Tokyo 153-8904, Japan
Sho Shirasaka
Graduate School of Information Science and Technology,Osaka University, 1-5 Yamadaoka, Suita, Osaka 565-0871, Japanand Research Center for Advanced Science and Technology,The University of Tokyo, 4–6–1 Komaba, Meguro-ku, Tokyo 153-8904, Japan
Akihiko Akao
Department of Mathematics, University of Pittsburgh, Pittsburgh, Pennsylvania 15213, USAand Graduate School of Engineering, The University of Tokyo,7–3–1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan
Yasuhiko Jimbo
Graduate School of Engineering, The University of Tokyo,7–3–1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan
Hiroya Nakao
Department of Systems and Control Engineering,Tokyo Institute of Technology, Tokyo 152-8522, Japan (Dated: July 23, 2020) a Electric adress: [email protected] bstract Spontaneous oscillations induced by time delays are observed in many real-world systems. Phase reductiontheory for limit-cycle oscillators described by delay-differential equations (DDEs) has been developed toanalyze their synchronization properties, but it is applicable only when the perturbation applied to theoscillator is sufficiently weak. In this study, we formulate a nonlinear phase-amplitude reduction theory forlimit-cycle oscillators described by DDEs on the basis of the Floquet theorem, which is applicable when theoscillator is subjected to perturbations of moderate intensity. We propose a numerical method to evaluatethe leading Floquet eigenvalues, eigenfunctions, and adjoint eigenfunctions necessary for the reduction andderive a set of low-dimensional nonlinear phase-amplitude equations approximately describing the oscillatordynamics. By analyzing an analytically tractable oscillator model with a cubic nonlinearity, we show that theasymptotic phase of the oscillator state in an infinite-dimensional state space can be approximately evaluatedand non-trivial bistability of the oscillation amplitude caused by moderately strong periodic perturbationscan be predicted on the basis of the derived phase-amplitude equations. We further analyze a model of gene-regulatory oscillator and illustrate that the reduced equations can elucidate the mechanism of its complexdynamics under non-weak perturbations, which may be relevant to real physiological phenomena such ascircadian rhythm sleep disorders.
Keywords: time-delayed systems | Floquet theorem | limit-cycle oscillations | synchronization | adjoint method | dimensionality reduction | multistability . INTRODUCTION Time-delayed feedback can break continuous time-translational symmetry and lead to oscillatorybehavior in many physical, biological, social, and engineered systems [1–8]. In biology, for example,ultradian oscillations in the hypothalamic-pituitary-adrenal (HPA) system are induced by time-delayed synthesis of hormones in the adrenal cortex [9]. Also, somite segmentation in zebrafishis regulated by oscillatory dynamics induced by time delays in the synthesis of proteins [10], andmammalian circadian rhythm is generated by feedback regulations of clock genes in suprachiasmaticnucleus (SCN) [11]. Such oscillatory dynamics in systems with time delays can be described asstable limit-cycle orbits of delay-differential equations (DDEs).In many of such systems, each oscillatory unit, or oscillator, is not isolated but perturbedby external forcing or by mutual coupling with other oscillators, and the state of each oscillatormay deviate from the unperturbed limit cycle of an isolated oscillator when the perturbation isnot sufficiently weak. Therefore, it is important to understand how perturbations of moderateintensity can modulate the period, amplitude, and other properties of delay-induced oscillations.For example, in the case of zebrafish somite segmentation, it is known that strong couplingsbetween cells are necessary for the spatio-temporal oscillatory dynamics of her1 (zebrafish hairy-related gene1 ) expression [10]. In the case of circadian clock genes, the oscillatory period in the free-running condition is known to be slightly different from 24 h, but they are entrained by the periodicexternal day-and-night lights through retinal ganglion cells [12, 13]. Strong light stimulation canfurther induce large modulation in the activities of the clock genes [14]. Since irregular dynamics ofcircadian rhythms manifest themselves as diseases such as sleep disorders [13, 15–17], understandingof the dynamics of circadian clock genes under strong perturbations may facilitate therapies forsleep disorders.The phase reduction theory is a standard mathematical framework for characterizing responseproperties of weakly perturbed limit-cycle oscillators and analyzing their synchronization dynamicsvia dimensionality reduction [18–24]. Recently, the phase reduction theory has been extendedalso to DDEs exhibiting stable limit-cycle oscillations, which requires non-trivial mathematicalgeneralization because DDEs are infinite-dimensional dynamical systems [25, 26]. However, thephase reduction has a strong limitation in that it is applicable only when the oscillator stateremains sufficiently close to the unperturbed limit cycle. Specifically, when non-weak perturbationsare applied or relaxation time of the system state to the limit cycle is not sufficiently small, theamplitude degrees of freedom may no longer be enslaved by the phase, leading to the breakdown3f the lowest-order phase-only description. In such cases, the nonlinear interaction of the phaseand amplitude may lead to non-trivial dynamics that cannot be captured by the phase reduction.To overcome this difficulty, several mathematical frameworks have been proposed for oscillatorysystems described by ordinary differential equations (ODEs), such as higher-order phase resettingcurves [27], extended phase equations [28], and higher-order phase-amplitude equations [29]. Still,for oscillatory dynamics of DDEs away from the limit cycle, much remains unknown because of theirinfinite-dimensional nature. Thus, a general framework for dimensionality reduction of limit-cycleoscillators described by DDEs that can analyze the effect of moderately strong perturbations isneeded. Such a framework would shed light on oscillatory dynamical systems in which nonlinearity,time delay, and strong perturbations coexist.In this study, our interest lies in the situation where the phase and amplitude of DDEs interactsignificantly in a nonlinear manner. We develop a nonlinear phase-amplitude reduction theoryfor DDEs, which gives a general mathematical framework for reducing DDEs describing limit-cycleoscillators to low-dimensional ordinary differential equations on the basis of the Floquet theory [30–32]. We also propose a practical numerical method, which we call the extended adjoint method, toevaluate the Floquet eigenvalues, eigenfunctions, and their adjoint functions, which are necessaryfor the reduction. By using biorthogonality of the Floquet eigenfunctions and their adjoints, weproject the oscillator state onto an eigenspace spanned by a few slowly-decaying Floquet eigenfunc-tions and derive a set of phase-amplitude equations which takes into account nonlinear interactionsbetween the slowly-decaying Floquet eigenmodes. In contrast to the standard lowest-order phasereduction, the amplitude component associated with the second Floquet eigenfunction is included,which can play important roles when the relaxation of the system is slow or when the system isstrongly perturbed.We confirm the validity of the theory using an analytically-tractable DDE with a cubic non-linearity by showing that the reduced phase-amplitude equations accurately predict the amplitudeof the phase-locked oscillations under a periodic force, which exhibits non-trivial bistable responseinduced by the non-weak amplitude effects. We then apply the theory to a model of a gene-regulatory oscillator under moderately strong forcing and analyze its synchronization dynamics.We show that the reduced phase-amplitude equations can also predict nontrivial bistable dynamicsof the system, which is analogous to a circadian disorder called advanced sleep-phase syndrome(ASPS) [15–17]. 4
I. THEORY
In this section, we derive a set of reduced nonlinear phase-amplitude equations for limit-cycleoscillators described by DDEs on the basis of the Floquet theory and propose a practical numericalmethod to calculate the Floquet eigenvalues, eigenfunctions and their adjoints that are necessary forthe reduction. We also derive approximate phase-amplitude equations for the oscillators subjectedto periodic external forcing.
A. DDEs with a stable limit-cycle solution
We consider general delay-differential equations (DDEs) that have a stable limit-cycle solution.Mathematical analysis of such DDEs, for example, analyzing the synchronization properties whenthey are periodically perturbed, is not easy because they describe nonlinear infinite-dimensionaldynamical systems on Banach spaces. Our aim is to derive simpler tractable equations by reducingthem to low-dimensional ODEs while preserving their essential quantitative properties and to an-alyze synchronization dynamics of nonlinear oscillators described by such DDEs under moderatelystrong external perturbations. In previous studies [25, 26], phase reduction methods for stablelimit-cycle solutions of DDEs have been developed, which are applicable when the perturbationsgiven to the system is sufficiently weak. In this study, we develop a nonlinear phase-amplitudereduction theory for DDEs.We consider a DDE for X ( t ) ∈ R N , represented as a column vector, with a maximum delaytime τ >
0. To construct a solution of the DDE, we have to take into account the history of X ( t ) from t − τ to t . Thus, we introduce its history-function representation, X ( t ) ( σ ) ≡ X ( t + σ )( − τ ≤ σ ≤
0) [30–32]. Here, X ( t ) ( · ) ∈ C and C = C ([ − τ, → R N ) is a Banach space of(column) vector-valued continuous functions mapping [ − τ,
0] into R N , which is equipped with anorm || x || C = sup θ ∈ [ − τ, || x ( θ ) || , where || · || is the usual Euclidean norm on R N . This historyfunction X ( t ) represents the state of the dynamical system described by the DDE at time t , wherethe state space of the system is given by the infinite-dimensional Banach space C .Using the above notation, a DDE can generally be written as ddt X ( t ) ( σ ) = ddσ X ( t ) ( σ ) ( − τ ≤ σ < , N ( X ( t ) ( · )) + G (cid:0) X ( t ) ( · ) , t (cid:1) ( σ = 0) . (1)Here, the vector-valued functional N : C → R N represents the system dynamics and G : C × R → N denotes external perturbation applied to the system that depends on the system state X ( t ) .Both functionals are assumed to be sufficiently smooth. This DDE can describe not only systemswith discrete delays but also systems with distributed delays [33]; see Ref. [34] for the relationbetween nonlinear functionals and their kernel representations, which are widely used for systemswith distributed delays described by integro-differential equations.We consider a situation in which the DDE (1) without the external perturbation ( G = 0) hasa stable limit-cycle solution X ( t ) whose period is T , i.e., X ( t + T ) = X ( t ), and represent it as ahistory function X ( t )0 ( · ) ∈ C satisfying X ( t + T )0 = X ( t )0 , where X ( t )0 ( σ ) ≡ X ( t + σ ) ( − τ ≤ σ ≤ . (2)In what follows, we also denote the limit cycle as X ( φ )0 , where we use the phase φ (0 ≤ φ < T ) inplace of the time t to parametrize system state on the limit cycle. The phase φ increases from 0to T , where the origin φ = 0 can be chosen as a specific system state on the limit cycle. Whenthe system state evolves along the limit cycle without perturbation, the phase φ increases with aconstant frequency 1, i.e., φ = t (mod T ). Similarly, we also denote T -periodic history functions,such as the Floquet eigenfunctions, using the phase φ instead of t when necessary.The definition of the phase can further be extended to the basin of attraction of the limitcycle by assigning the same phase value φ to the set of system states { X ( t ) } that asymptoticallyconverge to the same system state as X ( φ )0 when the system evolves without perturbation [25, 26],i.e., lim t →∞ k X ( t ) − X ( φ + t )0 k C = 0, yielding the notion of asymptotic phase Φ( X ( t ) ) : C → [0 , T )that maps a system state X ( t ) in the basin to a phase value. The asymptotic phase Φ satisfies ddt Φ( X ( t ) ) = 1 (3)when the system state evolves in the basin of the limit cycle without perturbation. The iso-surfaces of Φ, called isochrons , are not simply hyperplanes in general. For ordinary differentialequations, the asymptotic phase has been used as a canonical representation of rhythms of stableoscillatory dynamics [18–21, 35] and provides in-depth insights into strongly-perturbed oscillatorydynamics [24, 27–29, 36]. Recently, it has also been defined for DDEs and other non-conventionaloscillatory systems [25, 26].We assume that the relaxation dynamics of the system state to the limit cycle can be decomposedinto a few slow modes and remaining faster modes, which are well separated in time scale from eachother. In this case, a rectangular coordinate frame moving along the periodic orbit, which was usedin Refs. [35, 37, 38], is not useful for reducing the dynamics to low-dimensional ODEs, because fast6nd slow components interact already at the lowest order in this coordinate frame. It is also noteasy to proceed with the asymptotic phase and associated amplitudes, because they are generallygiven by highly nonlinear functionals of the system state X ( t ). We therefore use a coordinate framedefined by the Floquet eigenfunctions to decompose the system state as discussed in Ref. [23] forODEs. The space spanned by the Floquet eigenfunctions with non-vanishing relaxation rates istangent to the isochron at each point on the limit cycle. For this purpose, we need to calculate theFloquet eigenvalues, eigenfunctions, and their adjoints of DDEs. B. Floquet theory for DDEs
We first describe the Floquet theory for the DDE (1) without the perturbation term, i.e., G = 0. We denote small deviation of X ( t ) from X ( t ) as Y ( t ) = X ( t ) − X ( t ), and introduce itshistory-function representation Y ( t ) ( · ) ∈ C with Y ( t ) ( σ ) ≡ Y ( t + σ ) ( − τ ≤ σ ≤
0) as Y ( t ) ( σ ) = X ( t ) ( σ ) − X ( t )0 ( σ ) ( − τ ≤ σ ≤ . (4)The linearized variational equation for Y ( t ) is given by ddt Y ( t ) ( σ ) = L ( t ) (cid:16) Y ( t ) (cid:17) ( σ ) ( − τ ≤ σ ≤ , (5)where L ( t ) ( Y ( t ) ) is a history representation of a linear functional defined by L ( t ) (cid:16) Y ( t ) (cid:17) ( σ ) = ddσ Y ( t ) ( σ ) ( − τ ≤ σ < , ˆ − τ dσ ′ ¯Ω ( t ) ( σ ′ ) Y ( t ) ( σ ′ ) ( σ = 0) . (6)Here, ¯Ω ( t ) ( σ ) ≡ δ N ( X ( t ) ( · )) δX ( t ) ( σ ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) X ( t ) = X ( t )0 (7)is a functional differentiation of N with respect to X ( t ) evaluated at the system state X ( t ) = X ( t )0 on the limit cycle. Note that Eq. (5) gives a periodically driven linear system because X ( t )0 is T -periodic. In what follows, we expand N in a functional Taylor series in Y ( t ) as N ( X ( t ) ( · )) = N ( X ( t )0 ( · )) + L ( t ) (cid:16) Y ( t ) (cid:17) (0) + F nl (cid:16) Y ( t ) ( · ) (cid:17) , (8)where L ( t ) (cid:0) Y ( t ) (cid:1) (0) represents a linear functional of Y ( t ) defined in Eq. (6) with σ = 0 and F nl ( Y ( t ) ( · )) represents the remaining nonlinear functional of Y ( t ) , respectively, and both of thesefunctionals are evaluated at X ( t ) = X ( t )0 . 7s an example, let us consider a simple DDE, ddt X ( t ) = N ( X ( t ) , X ( t − τ )) , (9)which is equivalent to ddt X ( t ) ( σ ) = ddσ X ( t ) ( σ ) ( − τ ≤ σ < , N ( X ( t ) (0) , X ( t ) ( − τ )) ( σ = 0) , (10)in the history-function representation. By using the chain rule for functional differentiation andrepresenting the terms in N as X ( t ) (0) = ˆ − τ dσ ′ δ ( σ ′ ) X ( t ) ( σ ′ ) (11)and X ( t ) ( − τ ) = ˆ − τ dσ ′ δ ( σ ′ + τ ) X ( t ) ( σ ′ ) (12)using Dirac’s delta function δ ( · ), we obtain ¯Ω ( t ) ( σ ) = N ( t ) δ ( σ ) + N ( t ) δ ( σ + τ ) , (13)where N j ( t ) ≡ ∂ x j N ( x , x ) ( j = 1 ,
2) is evaluated at ( x , x ) = ( X ( t )0 (0) , X ( t )0 ( − τ )). The linearizeddynamics for the deviation Y ( t ) can then be written as L ( t ) (cid:16) Y ( t ) (cid:17) ( σ ) = ddσ Y ( t ) ( σ ) ( − τ ≤ σ < , N ( t ) Y ( t ) (0) + N ( t ) Y ( t ) ( − τ ) ( σ = 0) . (14)Let us introduce a time-periodic linear operator ˆ L of period T , which acts on a complexifiedBanach space ( C ) C [39, Sec. III.7] as (cid:16) ˆ LY ( t ) (cid:17) ( σ ) ≡ − ddt Y ( t ) ( σ ) + L ( t ) (cid:16) Y ( t ) (cid:17) ( σ ) ( − τ ≤ σ ≤ , (15)and rewrite Eq. (5) as ( ˆ LY ( t ) )( σ ) = 0. Because Y ( t ) obeys a periodically driven linear system, bythe Floquet theorem for linear DDEs [30–32], the spectrum of ˆ L is at most countable and (cid:16) ˆ Lq ( t ) i (cid:17) ( σ ) = λ i q ( t ) i ( σ ) ( − τ ≤ σ ≤
0) (16)is satisfied, where λ i ∈ C is the i -th Floquet eigenvalue and q ( t ) i ∈ ( C ) C is the corresponding T -periodic Floquet eigenfunction ( i = 0 , , , ... ). Here, the largest eigenvalue, which is 0 and simple8y the Floquet theorem, is denoted as λ = 0 and the other eigenvalues are arranged in descendingorder of the real part.We also introduce adjoint eigenfunctions with respect to a bilinear form appropriate forDDEs [40]. Following Refs. [30–32], we define a bilinear form of two functions, A ∈ ( C ) C and B ∈ ( C ) ∗ C , as h B ( t ) , A ( t ) ; t i ≡ h B ( t ) (0) , A ( t ) (0) i − ˆ − τ dσ ˆ σ dξ h B ( t ) ( ξ − σ ) , ¯Ω ( t + ξ − σ ) ( σ ) A ( t ) ( ξ ) i . (17)Here, ( C ) ∗ C = C ([0 , τ ] → C N ∗ ) is the dual space of ( C ) C with respect to the bilinear form,consisting of (row) vector-valued functions that map the interval [0 , τ ] to C N ∗ , and [ · , · ] denotesthe Hermitian scalar product of V ∈ C N ∗ and U ∈ C N defined as [ V, U ] = P Nk =1 V k U k where V k and U k are vector components of V and U , respectively. An adjoint operator ˆ L ∗ of ˆ L with respectto this bilinear form can then be derived as (cid:16) ˆ L ∗ Y ( t ) ∗ (cid:17) ( s ) = ddt Y ( t ) ∗ ( s ) + L ( t ) ∗ (cid:16) Y ( t ) ∗ (cid:17) ( s ) (0 ≤ s ≤ τ ) , (18)where L ( t ) ∗ (cid:16) Y ( t ) ∗ (cid:17) ( s ) = − dds Y ( t ) ∗ ( s ) (0 < s ≤ τ ) , ˆ τ ds ′ Y ( t ) ∗ ( s ′ ) ¯Ω ( t + s ′ ) ( − s ′ ) ( s = 0) . (19)Here, Y ∗ ( t ) ∈ C N ∗ is a row vector of N complex components and Y ( t ) ∗ ( s ) ≡ Y ∗ ( t + s ) (0 ≤ s ≤ τ ) ∈ ( C ) ∗ C is its history-function representation.The adjoint eigenfunction q ( t ) ∗ i ∈ ( C ) ∗ C of q ( t ) i , which is also T -periodic, satisfies (cid:16) ˆ L ∗ q ( t ) ∗ i (cid:17) ( s ) = ¯ λ i q ( t ) ∗ i ( s ) (0 ≤ s ≤ τ ) , (20)where ¯ λ i is the complex conjugate of λ i . If λ i = λ j , q ( t ) i is orthogonal to q ( t ) ∗ j with respect tothe bilinear form Eq. (17), and hence they can be normalized to satisfy the biorthogonal relation h q ( t ) ∗ i , q ( t ) j ; t i = δ i,j . The zero eigenfunction of the linear operator ˆ L can be chosen as q ( t )0 ( σ ) = dX /dt | t + σ ( − τ ≤ σ ≤ t at X ( t ) = X ( t )0 on the periodic orbit [25]. Note that this definition specifies the normalization of q ( t )0 .For the other eigenfunctions q ( t ) i ( i = 1 , , ... ), we normalize them such that max ≤ t ≤ T (cid:16) q ( t ) i (0) (cid:17) =1. We use this convention for the normalization throughout this study. We note that the zeroeigenfunction of the linear operator ˆ L corresponds to the tangential component along the limitcycle, namely, the phase direction. Moreover, the zero eigenfunction q ( t ) ∗ of the adjoint operator ˆ L ∗ phase sensitivity function of the limit cycle, which characterizes linear response propertyof the oscillator phase to weak perturbations [25, 26]. Similarly, the other eigenfunctions q ( t ) ∗ i ( i = 1 , , ... ) characterize linear response properties of the amplitudes and called isostable responsecurves for the case of ODEs [41].The adjoint eigenfunctions can numerically be obtained by an extension of the adjoint methodfor DDEs, which was previously used to calculate the adjoint zero eigenfunction of ˆ L [25, 26]. Thatis, we numerically integrate the linearized and its adjoint equations while subtracting unnecessaryfunctional components by using the biorthogonality between the eigenfunctions and adjoint eigen-functions. The main difference from the adjoint method for q ( t ) ∗ developed in the previous studiesis that we calculate the adjoint eigenfunctions also for λ i ( i ≥ i − i -th component in order tocalculate the i -th eigenfunction precisely. For i ≥
1, we also need to renormalize the solutions ofthe equations by a factor e λ i t determined by the Floquet exponent in order to obtain the correcteigenfunctions.To numerically calculate the i -th eigenfunction q ( t ) i , we integrate the linearized equation ddt Y ( t ) ( σ ) = L ( t ) (cid:16) Y ( t ) (cid:17) ( σ ) ( − τ ≤ σ ≤
0) (21)forward in time. During the calculation, we subtract the 0-th to ( i − Y ( t ) , which are unnecessary but arises due to numerical errors. The Floquet eigenvalue λ i is numerically evaluated from the exponential decay rate of Y ( t ) . Then the eigenfunction q ( t ) i ( σ ) isobtained by compensating the exponential decay of Y ( t ) ( σ ) as q ( t ) i ( σ ) = e − λ i t Y ( t ) ( σ ) ( − τ ≤ σ ≤ i -th adjoint eigenfunction q ( t ) ∗ i is calculated by numerically integrating the adjoint linear equation ddt Y ( t ) ∗ ( s ) = − L ( t ) ∗ (cid:16) Y ( t ) ∗ (cid:17) ( s ) (0 ≤ s ≤ τ ) (22)backward in time while subtracting unnecessary eigencomponents and then compensating the nu-merical result by a factor e − λ i t . We call this procedure the extended adjoint method in this study. C. Nonlinear phase-amplitude equations
Our aim is to derive a set of low-dimensional dynamical equations from the original DDE byprojecting the system state onto a moving coordinate frame spanned by a small number of Floquet10igenfunctions. That is, we decompose the deviation of the system state X ( t ) from that on thelimit cycle X ( t )0 by using the eigenfunctions associated with the leading M eigenvalues other than0, which are assumed to be real and simple for the sake of simplicity [43], as X ( t ) ( σ ) ≃ X ( φ )0 ( σ ) + M X i =1 ρ i ( t ) q ( φ ) i ( σ ) , ( − τ ≤ σ ≤ , (23)where X ( φ )0 is a system state on the limit cycle parametrized by the phase φ ∈ [0 , T ), q ( φ ) i ( i =1 , ..., M ) is the Floquet eigenfunction associated with λ i and denoted as a function of φ rather than t , and { ρ i ( t ) } are real expansion coefficients representing amplitudes of the Floquet eigenmodes.The symbol ≃ indicates that we approximate X ( t ) ( σ ) by its projection on the space spanned bythe M eigenfunctions { q ( φ )1 , ..., q ( φ ) M } . We here use the term “amplitude” in a generalized sense,allowing it to take both positive and negative values; it is the component of the system statealong the Floquet eigenfunction corresponding to the direction transversal to the limit cycle andrepresents the deviation of the system state from the limit cycle. Here, the phase value φ fora given state X ( t ) is determined in such a way that the state difference X ( t ) − X ( φ )0 does nothave a tangential functional component q ( φ )0 along the limit cycle. Thus, we assume the followingorthogonality condition: D q ( φ ) ∗ , X ( t ) − X ( φ )0 ; φ E = 0 , (24)namely, the difference X ( t ) − X ( φ )0 is on the hyperplane that is tangent to the isochron on the limitcycle at X ( φ )0 . Note that the phase defined in this way is different from the asymptotic phase.Because we use a linear coordinate frame spanned by the Floquet eigenfunctions { q ( φ ) i } ( i =1 , ..., M ), nonlinear interactions between different eigenmodes generally arise. Specifically, whenthe eigenvalue λ with the largest non-zero part is close to 0, the perturbed system state does notgo back to the limit cycle quickly, and hence nonlinear interactions between the phase eigenmodeand the slowest-decaying amplitude eigenmode should be taken into account for better descriptionof the system.For ordinary differential equations, such coupled nonlinear phase-amplitude equations have beenderived by transforming the original equations around the limit cycle in several contexts [37, 46].Such transformation methods have also been developed for DDEs in Refs. [47, 48], though thetreatments of oscillatory dynamics in these studies are rather abstract. Quantitative analysis ofsynchronization dynamics of DDEs using the coordinate transform proposed therein have not beenvery fruitful despite their potential advantages, mainly due to the lack of practical methods fornumerically evaluating the Floquet eigenfunctions.11e hereafter restrict ourselves to the case in which λ takes a negative real value near zeroand Re { λ } ≪ λ for simplicity. To derive the phase and amplitude equations, we retain only theslowest two modes associated with λ and λ and approximate X ( t ) ( σ ) as X ( t ) ( σ ) ≃ X ( φ )0 ( σ ) + R ( t ) q ( φ )1 ( σ ) , (25)where R ( t ) = ρ ( t ) is the amplitude of the eigenmode corresponding to λ . The symbol ≃ hereindicates that we further approximate X ( t ) ( σ ) by its projection on a one-dimensional space spannedby q ( φ )1 . We substitute this expression into Eq. (1) and then project both sides of Eq. (1) onto theeigenfunctions q ( φ )0 and q ( φ )1 , respectively, by using biorthogonality of the eigenfunctions and derivethe equations for the phase φ and the amplitude R .As explained in Appendix A, the phase equation can be derived as dφdt = 1 + q ( φ ) ∗ (0) · ( F nl ( φ, R ) + G ( φ, R, t ))1 + R D q ( φ ) ∗ , L ( φ ) (cid:16) q ( φ )1 (cid:17) ; φ E , (26)or, by rewriting the right-hand side, dφdt =1 + q ( φ ) ∗ (0) · ( F nl ( φ, R ) + G ( φ, R, t )) − R D q ( φ ) ∗ , L ( φ ) (cid:16) q ( φ )1 (cid:17) ; φ E R D q ( φ ) ∗ , L ( φ ) (cid:16) q ( φ )1 (cid:17) ; φ E q ( φ ) ∗ (0) · ( F nl ( φ, R ) + G ( φ, R, t )) , (27)and the amplitude equation can similarly be derived as dRdt = λ R + q ( φ ) ∗ (0) · ( F nl ( φ, R ) + G ( φ, R, t )) − R hD q ( φ ) ∗ , L ( φ ) (cid:16) q ( φ )1 (cid:17) ; φ E − λ i R D q ( φ ) ∗ , L ( φ ) (cid:16) q ( φ )1 (cid:17) ; φ E q ( φ ) ∗ (0) · ( F nl ( φ, R ) + G ( φ, R, t )) , (28)where the nonlinear functional N in Eq. (8) is approximated by an ordinary function of φ and R , F nl ( φ, R ) ≡ F nl (cid:16) Rq ( φ )1 ( · ) (cid:17) = N (cid:16) X ( φ )0 ( · ) + Rq ( φ )1 ( · ) (cid:17) − N (cid:16) X ( φ )0 ( · ) (cid:17) − L ( φ ) (cid:16) Rq ( φ )1 (cid:17) (0) , (29)and the external perturbation is also approximated as G ( φ, R, t ) ≡ G (cid:16) X ( φ )0 ( · ) + Rq ( φ )1 ( · ) , t (cid:17) . (30)In Eq. (27) and Eq. (28), both the second and third terms on the right-hand side depend on F nl and G . Note that F nl ( φ, R ) includes only terms of O ( R ) or higher, because the constant termsand linear terms in R have already been subtracted in Eq. (29).12hus, by projecting the DDE onto the first two eigenfunctions, a set of two-dimensional coupledordinary differential equations for the phase φ and amplitude R is obtained. In order to considerthe higher-order effects of the amplitude deviations, we have not expanded the third-order termsin Eq. (27) and Eq. (28) in a series of R and hence the dynamics of φ and R are nonlinearlycoupled at the second and higher orders in R . This nonlinearity can be a source of intriguingoscillatory dynamics [27–29, 36]. We also note that the lowest-order phase-amplitude equations(see Refs. [23, 29, 41] for the case of ODEs) dφdt =1 + q ( φ ) ∗ (0) · G ( φ, R, t ) dRdt = λ R + q ( φ ) ∗ (0) · G ( φ, R, t ) (31)are obtained at the lowest-order approximation in R , where F nl ( φ, R ) is O ( R ) and does not appearat the lowest order.Finally, before proceeding, we note that there are also other formulations of phase or phase-amplitude reduced equations for analyzing higher-order effects of perturbations on limit cyclesdescribed by ODEs, such as non-pairwise phase interactions [23], higher-order phase reduction [49],nonlinear phase coupling function [50], and higher-order approximations of coupling functions [41],which can capture more detailed aspects of synchronization than the lowest-order phase equation. D. Averaged phase-amplitude equations
When the perturbation applied to the oscillator is a periodic external force whose frequencyis close to the natural frequency of the oscillator, we may further derive simpler, approximatephase-amplitude equations by averaging out the fast oscillatory component as follows.We assume that the perturbation G is purely external (i.e. independent of the system state andperiodic in t with period T ′ = T /r (frequency r ), i.e., G ( t + T /r ) = G ( t ) . (32)We also assume that the detuning between the natural frequency of the oscillator and the periodicforce is small and denote it as ∆ ω = 1 − r .We introduce a slow phase variable ψ ≡ φ − rt . The equations for ψ and R can be written as dψdt = ∆ ω + q ( ψ + rt ) ∗ (0) · ( F nl ( ψ + rt, R ) + G ( t ))1 + R D q ( ψ + rt ) ∗ , L ( ψ + rt ) (cid:16) q ( ψ + rt )1 (cid:17) ; ψ + rt E , (33)13nd dRdt = λ R + q ( ψ + rt ) ∗ (0) · ( F nl ( ψ + rt, R ) + G ( ψ + rt, R, t )) − R hD q ( ψ + rt ) ∗ , L ( ψ + rt ) (cid:16) q ( ψ + rt )1 (cid:17) ; ψ + rt E − λ i R D q ( ψ + rt ) ∗ , L ( ψ + rt ) (cid:16) q ( ψ + rt )1 (cid:17) ; ψ + rt E q ( ψ + rt ) ∗ (0) · ( F nl ( ψ + rt, R ) + G ( ψ + rt, R, t )) . (34)We also expand the nonlinear term F nl in Taylor series of R up to R N as F nl ( ψ + rt, R ) = N X ℓ =2 R ℓ F nl ,ℓ ( ψ + rt ) + O ( R N +1 ) , (35)where { F nl ,ℓ } ( ℓ = 2 , , ... ) are expansion coefficients. Note that the series for F nl starts from O ( R ).Considering that ψ evolves only slowly while rt rapidly increases, we approximate the termswith ψ + rt in Eqs. (33) and (34) by their one-period average, for example, as q ( ψ + rt ) ∗ (0) · F nl , ( ψ + rt ) ≈ T ′ ˆ T ′ q ( ψ + rs ) ∗ (0) · F nl , ( ψ + rs ) ds = 1 T ˆ T q ( θ ) ∗ (0) · F nl , ( θ ) dθ = a (36)and q ( ψ + rt ) ∗ (0) · G ( t ) ≈ T ′ ˆ T ′ q ( ψ + rs ) ∗ (0) · G ( s ) ds = 1 T ˆ T q ( θ ) ∗ (0) · G (cid:18) θ − ψr (cid:19) dθ = g ( ψ ) , (37)where ψ is kept constant during the integration. Expanding F nl ( ψ + rt, R ) up to O ( R ) andaveraging the coefficients, we obtain approximate equations for ψ and R as dψdt = ∆ ω + 11 + Ra (cid:0) a R + a R + g ( ψ ) (cid:1) , (38)and dRdt = λ R + b R + b R − R ( b − λ )1 + Ra (cid:0) a R + a R + g ( ψ ) (cid:1) + g ( ψ ) , (39)where the equations for the individual coefficients are given in Appendix B. We check the validityof the above averaging approximation numerically in the next section. E. Evaluation of the phase and amplitude
Numerically, the values of the phase φ and amplitude R can be evaluated from the system state X ( t ) by the following two-step procedure. First, we evaluate the phase of the state X ( t ) by choosing14he phase value φ so that it satisfies the orthogonality condition Eq. (24). Numerically, we find thevalue ˆ φ that minimizes the mean squared error, (cid:12)(cid:12)(cid:12)(cid:12)(cid:28) q ( ˆ φ ) ∗ , X ( t ) − X ( ˆ φ ) ; ˆ φ (cid:29)(cid:12)(cid:12)(cid:12)(cid:12) . (40)There exists a neighborhood of the periodic orbit where the phase and amplitude componentsdefined by using the Floquet eigenfunctions are uniquely determined [48, Lemma1]. However, ingeneral, there can exist multiple values of ˆ φ that satisfy Eq. (24) in the range 0 ≤ ˆ φ < T . Tochoose the appropriate value from them, for each candidate of ˆ φ , we evaluate the corresponding q component as ˆ R = (cid:28) q ( ˆ φ ) ∗ , X ( t ) − X ( ˆ φ ) ; ˆ φ (cid:29) (41)and adopt the value of ˆ φ that has the smallest (cid:12)(cid:12)(cid:12) ˆ R (cid:12)(cid:12)(cid:12) as the estimate of φ , and the smallest ˆ R as theestimate of R . F. Approximate evaluation of the asymptotic phase
The phase φ defined by the Floquet eigenfunction, which we use in the present study forthe phase-amplitude description, is different from the asymptotic phase Φ; the isosurface of Φ isgenerally curved and tangent to the isophase plane of φ at each point on the limit cycle. Sincethe asymptotic phase Φ provides useful information on the nonlinear dynamical properties of theoscillator, it is convenient if we can approximate Φ using φ and R . In this subsection, we proposea method to approximately evaluate the asymptotic phase of an unperturbed oscillator from φ and R defined by the Floquet eigenfunctions, which is valid when R is sufficiently small.When the perturbation is absent ( G = 0), Eq. (26) for φ can be written as dφdt = 1 + d ( φ, R ) (42)where d ( φ, R ) = q ( φ ) ∗ (0) · F nl ( φ, R )1 + R D q ( φ ) ∗ , L ( φ ) (cid:16) q ( φ )1 (cid:17) ; φ E . (43)The asymptotic phase Φ of the system state X ( t ) at time t with phase φ and amplitude R canapproximately be obtained by integrating d ( φ ( s ) , R ( s )) until the system state goes back sufficientlyclose to the limit cycle as Φ( X ( t ) ) = φ + ˆ ∞ t d ( φ ( s ) , R ( s )) ds. (44)15hen R is sufficiently small, we may ignore the higher-order terms in R in the equations for φ and R and assume that φ increases constantly with frequency 1 and R decays exponentially withrate λ as φ = φ + t − t , R ( t ) = R exp ( λ ( t − t )) , (45)at the lowest-order approximation. The asymptotic phase Φ of the system state X ( t ) can then beapproximately evaluated asˆΦ = φ + ˆ ∞ t d ( φ + s − t , R exp ( λ ( s − t ))) ds. (46)In Sec. III E and Sec. IV, we use the above method to estimate the asymptotic phase Φ of theoscillator and compare it with direct numerical results. III. ANALYTICALLY TRACTABLE MODEL
To demonstrate the validity of the proposed framework, we first consider a limit-cycle oscillatordescribed by a scalar DDE with a cubic nonlinearity, for which approximate expressions of theFloquet eigenfunctions and their adjoints can be analytically derived, and analyze the effect of aperiodic force on the dynamics.
A. Model
The model is represented as dx ( t ) dt = − x (cid:16) t − π (cid:17) + ǫx ( t ) (cid:20) − x ( t ) − x (cid:16) t − π (cid:17) (cid:21) + G ( t ) , (47)where x ( t ) ∈ R , ǫ = 0 .
05 is a small constant, and the external periodic force is described by G ( t ) = G sin (cid:18) πT rt (cid:19) , (48)where G is the intensity of the periodic force and r is the ratio of the natural frequency 2 π/T ofthe limit cycle to that of the external force. It is assumed that r is sufficiently close to 1.When G = 0, this DDE has a limit cycle of period T = 2 π given by x ( t ) = sin t , or x ( t )0 ( σ ) = sin( t + σ ) ( − τ ≤ σ ≤
0) (49)in the history-function representation, and its rate of attraction to the limit cycle is determined by ǫ . When ǫ is small, the relaxation time of the system state to the limit cycle is considerably largeas compared to the oscillation period as shown in Figs. 1(a) and (b).16e denote the small deviation of the system state from the limit cycle as y ( t ) ( σ ) = x ( t ) ( σ ) − x ( t )0 ( σ ) ( − τ ≤ σ ≤ L of this system is given by Eq. (6) with ¯Ω ( t ) ( σ ) = δ ( σ ) (cid:20) ǫ (1 − x ( t ) − x (cid:16) t − π (cid:17) (cid:21) − δ (cid:16) σ + π (cid:17) h ǫx ( t ) x (cid:16) t − π (cid:17)i , (50)where δ is Dirac’s delta function. By retaining the first two leading eigenvalues, the nonlinearphase-amplitude equations can be derived as Eqs. (27) and (28). B. Approximate analytical expressions of the eigenvalues and eigenfunctions
We first derive approximate Floquet eigenvalues, eigenfunctions, and adjoint eigenfunctions ofthe model Eq. (47) without the external force ( G = 0) analytically. In what follows, we considerthe case in which the relaxation of the system state to the limit cycle is slow and assume that λ is small and O ( ǫ ). First, the zero eigenfunction of ˆ L is given exactly as q ( t )0 ( σ ) = cos( t + σ ) ( − τ ≤ σ ≤
0) (51)and the adjoint eigenfunction is q ( t ) ∗ ( s ) = 8 ǫπ + 4 cos( t + s ) (0 ≤ s ≤ τ ) . (52)To find the exponent λ with the second largest real part, we introduce an ansatz q ( t )1 ( σ ) = Ce λ σ (sin( t + σ ) + l cos( t + σ )) ( − τ ≤ σ ≤
0) (53)where l is a constant and plug this into Eqs. (5) and (6). We then obtain the approximate eigenvalueand the associated eigenfunction up to O ( ǫ ) as λ = − ǫπ + 4 (54)and q ( t )1 ( σ ) = 2 √ π e − ǫσπ (cid:16) sin( t + σ ) − π t + σ ) (cid:17) ( − τ ≤ σ ≤ , (55)respectively. Similarly, for the corresponding adjoint eigenfunction, we approximately obtain q ( t ) ∗ ( s ) = C ∗ e ǫsπ (cid:16) sin( t + s ) + π t + s ) (cid:17) (0 ≤ s ≤ τ ) , (56)where the constant C ∗ is determined from the normalization condition h q ( t ) ∗ , q ( t )1 ; t i = 1.17 . Numerical evaluation of the eigenvalues and eigenfunctions To confirm the validity of the approximate analytical results for the Floquet eigenvalues, eigen-functions, and adjoint eigenfunctions obtained in the previous subsection, we numerically evaluatethese quantities by the extended adjoint method and compare with the approximate analyticalresults.First, as in the conventional adjoint method for DDEs [25, 26], we compute q ( t )0 ( σ ) ( − τ ≤ σ ≤ dX /dt | t + σ , and then q ( t ) ∗ ( σ ) ( − τ ≤ σ ≤
0) by backwardly integrating the adjointlinear equation. The adjoint eigenfunction q ( t ) ∗ is normalized such that h q ( t ) ∗ , q ( t )0 ; t i = 1. Next,we obtain the eigenfunction q ( t )1 with the largest negative eigenvalue ( λ < λ > λ i for i =2 , · · · , M ) [43]. As an initial function, we take an arbitrary function Y ( t =0)ini at t = 0 [51], subtractthe q ( t =0)0 component from this initial function as Y ( t =0) ( σ ) = Y ( t =0)ini ( σ ) − h q (0)0 ∗ , Y (0)ini ; 0 i q (0)0 ( σ )( − τ ≤ σ ≤ Y ( t =0)ini onto q (0)0 , and numericallyintegrate the linear equation (21) for Y ( t ) from this initial condition to t = T as explained before.Similarly, in order to compute the eigenfunction q ( t ) ∗ , we initialize Y ( t =0) ∗ ( s ) (0 ≤ s ≤ τ )appropriately and numerically integrate Eq. (22) backward in time, subtracting the q ( t ) ∗ componentat every period, and compensate the exponential decay in the numerical solution. The adjointeigenfunction q ( t ) ∗ is normalized so that h q ( t ) ∗ , q ( t )1 ; t i = 1.Figure 1(c) shows the exponential decay of the peak heights of Y ( t = nT ) (0), from which weobtain the Floquet eigenvalue λ . Figure 1(d) shows the time course of e − λ t Y ( t ) (0) that is usedfor numerical computation of eigenfunction q ( φ )1 . Figures 1(e) and (f) show the obtained pair ofFloquet eigenfunctions, where q ( φ )0 (0) and q ( φ ) ∗ (0) are plotted with respect to φ in Fig. 1(e), and q ( φ )1 (0) and q ( φ ) ∗ (0) are plotted with respect to φ in Fig. 1(f). We can confirm a good agreementbetween the numerical results and approximate analytical results for the eigenfunctions. Thenumerical value of the largest negative exponent λ is approximately evaluated as − . − ǫ/ ( π + 4) = − . . Phase-amplitude equations We now derive a set of nonlinear phase-amplitude equations from Eq. (47) with the periodicsinusoidal force. The nonlinear term F nl ( φ, R ) in Eq. (29) is explicitly given by F nl ( φ, R ) = ǫRq ( φ )1 (0) (cid:26) − (cid:16) x ( φ )0 (0) + Rq ( φ )1 (0) (cid:17) + (cid:16) x ( φ )0 (0) (cid:17) − (cid:16) x ( φ )0 (cid:16) − π (cid:17) + Rq ( φ )1 (cid:16) − π (cid:17)(cid:17) + (cid:16) x ( φ )0 (cid:16) − π (cid:17)(cid:17) (cid:27) + ǫx ( φ )0 (0) (cid:18) − (cid:16) Rq ( φ )1 (0) (cid:17) − (cid:16) Rq ( φ )1 (cid:16) − π (cid:17)(cid:17) (cid:19) (57)and the reduced equations (27) and (28) for φ and R can be derived using this equation.Expanding the nonlinear term F nl and applying the averaging procedure, the approximateequations for the phase difference ψ = φ − rt and R are given in the form of Eqs. (38) and (39)with g ( ψ ) = G T ˆ T q ( φ ) ∗ (0) · sin (cid:18) π ( φ − ψ ) T (cid:19) dφ = G (cid:18) g sin 2 πψT + g cos 2 πψT (cid:19) (58)and g ( ψ ) = ˆ T q ( φ ) ∗ (0) · sin (cid:18) π ( φ − ψ ) T (cid:19) dφ = G (cid:18) g sin 2 πψT + g cos 2 πψT (cid:19) . (59)Using numerically evaluated eigenvalues and eigenfunctions, the coefficients in Eqs. (38) and (39)can be calculated as λ = − . a = 1 . a = 0 . a = 0 . b = 1 . b = − . b = 0 . g = − . g = 0, g = − . g = 0 . ψ and the amplitude R are obtained as dψdt = ∆ ω + 11 + 1 . R (cid:18) . R + 0 . R − . G sin 2 πψT (cid:19) ,dRdt = λ R − . R + 0 . R − . G sin 2 πψT + 0 . G cos 2 πψT − R (1 . − λ )1 + 1 . R (cid:18) . R + 0 . R − . G sin 2 πψT (cid:19) . (60)Thus, we have approximately reduced an infinite-dimensional dynamical system described by aDDE to a set of ODEs for the phase and amplitude. E. Approximate evaluation of the asymptotic phase
In this subsection, we verify the validity of the approximate expression for the asymptotic phasederived in Sec. II F by evolving the present model from initial conditions far from the limit cycle.19rom the reduced phase-amplitude equations and Eq. (46), the asymptotic phase Φ for the presentmodel can be approximately evaluated from the phase φ and amplitude R asˆΦ = φ + (0 . . φ + 4 . . φ + 2 . R (61)up to O ( R ). For a given system state x ( t ) , the phase φ and amplitude R can be evaluated asexplained in Sec. II E, and the approximate asymptotic phase ˆΦ can then be obtained by Eq. (61).We also directly evaluate the asymptotic phase Φ for several initial conditions by numericallyintegrating the system and measuring the time necessary for the system state to converge sufficientlyclose to the limit cycle for comparison.As the first example, we try to estimate Φ when the initial function is on the φ - R plane, thatis, x ( t =0) ( σ ) = x ( φ )0 ( σ ) + Rq ( φ )1 ( σ ) ( − τ ≤ σ ≤ − φ for given initial valuesof φ and R obtained by direct numerical integration of the DDE, and Fig. 2 (b) shows analyticalresults of ˆΦ − φ obtained from Eq. (61). Figure 2 (c) shows the absolute difference between Φ andits analytical estimation ˆΦ. We can confirm a good agreement between the approximate analyticalcurve and direct numerical results for the whole range of φ when | R | is not too large.As the second example, we consider initial functions that are not on the φ - R plane. We set theinitial functions as x ( t =0) ( σ ) = sin σ + p sin( σ/
2) ( − τ ≤ σ ≤
0) with varying values of p [52], andevaluated their asymptotic phase Φ by direct numerical integration of the DDE. Figure 2 (d) showsthe phase φ , the asymptotic phase Φ estimated by Eq. (61), and the asymptotic phase Φ obtainedby direct numerical integration. We can confirm that the approximate analytical estimate of theasymptotic phase given by Eq. (61) gives reasonable agreement with the direct numerical resultseven though the system state is considerably far from the φ - R plane. F. Effect of a periodic force on the amplitude
In this subsection, we consider the effect of a periodic external force of moderate intensity withsmall frequency detuning. In particular, we focus on the average effect of the periodic force on theamplitude R in the phase-locked state, which cannot be analyzed without the amplitude equation.Since g = 0 in Eq. (58), the ψ -nullcline on which ˙ ψ = 0 is obtained from the averagedequation (38) as ψ = T π arcsin (cid:20) − Ra g G (cid:18) ∆ ω + 11 + Ra (cid:0) a R + a R (cid:1)(cid:19)(cid:21) (62)when the argument of arcsin is in the range [ − , R = F s ( R, G , ∆ ω ) = 0,20here F s represents the right-hand side of Eq. (39). The effect of the intensity of the periodic force G and the detuning ∆ ω on the stationary amplitude R of the oscillation in the steady state canbe evaluated from the partial derivatives of F s ( R, G , ∆ ω ) by the implicit function theorem.Figure 3 shows the predicted amplitude of the oscillation. The dependence of the amplitude on G at r = 1 is plotted in Fig. 3(a), where the stationary amplitude obtained from the averagedphase-amplitude equations (38) and (39) are compared with the linear approximation of the station-ary amplitude with a slope ∂R/∂G | R =0 ,G =0 , ∆ ω =0 = 18 .
2. Similarly, Fig. 3 (b) shows the depen-dence of the amplitude on r at G = 0 .
1, where the result of the phase-amplitude equations are com-pared with linear approximation of the amplitude with a slope ∂R/∂r | R =0 . ,G =0 . , ∆ ω =0 = 9 . G. Bistable response of delay-induced oscillation to a periodic force
In this subsection, we demonstrate that the present model can exhibit a nontrivial bistableresponse to a periodic force by a bifurcation analysis of Eq. (38) and Eq. (39). Such a phenomenonresults from higher-order amplitude effects and cannot be described by the phase-only equation northe lowest-order phase-amplitude equations. Using XPP-AUTO [53], we numerically find stationarysolutions in the range
R > − . / (1 + Ra ) exists (note that a = 1 . G and r , we observe quantitatively different behaviors of the systemstate as shown in Fig. 4.Figures 4 (a) and (b) show the stable and unstable fixed points on the ( R , r )-plane at twodifferent values of G . The system is always monostable when G = 0 .
02, while a bistable regionwhere R can take two stable fixed points is found around r = 1 .
052 when G = 0 .
1. Thus, itis expected that DDE (47) with G ( x, t ) = 0 . . t ) shows bistable dynamics. Figure 4 (c)shows the nullclines and stable fixed points on the ψ - R plane at r = 1 .
052 and G = 0 .
1. Thetwo crosses show the stable fixed points at ( − . , . − . , − . − ,
0) and ( − . , G ( x, t ) = 0 . . t ). We can clearly observethe bistable dynamics of the oscillator caused by moderately strong periodic forcing.21 V. GENE-REGULATORY OSCILLATOR
In this section, as a more complex, biologically-motivated example of DDEs, we investigate amodel of gene regulation [3] under a periodic sinusoidal force given by dx ( t ) dt = αC [ C + x ( t − τ )] − γx ( t ) R + x ( t ) − βx ( t ) + G ( t ) , (63)where x ( t ) ∈ R is the state variable representing protein concentration and α, β, γ, C , R , and thedelay time τ are real parameters. The first term of the right-hand side represents protein synthesiswith time delay for transcription and translation, while the second and the third terms representdegradation and dilution of the protein, respectively. Following the previous research [3], we set β = 0 . C = 10, and τ = 1. The external periodic force is G ( t ) = G sin (cid:0) πT rt (cid:1) with intensity G and frequency mismatch r . We set the rate constant of synthesis as α = 100, degradation as γ = 100, and Michaelis constant of degradation as R = 10 so that the system exhibits a slowconvergence to a limit cycle orbit and the effect of the amplitude dynamics can be clearly observed.This system has a stable limit cycle with a period T = 2 .
46, which can be obtained onlynumerically. Figures 5 shows the system state x ( t ) converging toward the limit-cycle attractor;Fig. 5(a) plots the time course of x ( t ) as a function of t and Fig. 5(b) shows the system trajectoryprojected on the ( x, dx/dt )-plane. The time constant of the relaxation to the limit cycle is muchlarger than the period of the oscillation as can be seen from the figures.For this model, the T -periodic linear operator ˆ L is given by Eq. (15) with ¯Ω ( t ) ( σ ) = δ ( σ ) (cid:26) − β − γR ( R + x ( t )) (cid:27) + δ ( σ + τ ) (cid:26) − αC ( C + x ( t − τ )) (cid:27) . (64)Figures 5 (c) and 5(d) show the first two eigenfunctions and adjoint eigenfunctions of ˆ L obtainedby the extended adjoint method [54], respectively. The second largest Floquet exponent is λ = − . a = 0 . a = − . × − , a = 1 . × − , g = − . g = − . × − , b = − . b = − . × − , b = 1 . × − , g = − . g = 0 . G = 0). We take the initial condition as a constant function, x ( t =0) ( σ ) = p , and evaluate the asymptotic phase by Eq. (46) and by direct numerical integrationof the DDE. Figure 6 (a) shows the phase φ , the asymptotic phase Φ estimated by using Eq. (46),and the asymptotic phase Φ evaluated by direct numerical integration of the DDE. It can be seen22hat the approximate analytical result reproduces the result of direct numerical measurement ofthe asymptotic phase.We next consider how the gene-regulatory oscillator behaves when it is subjected to a periodicexternal force. We conduct bifurcation analysis for different values of G and r in the same way asthat for Eq. (47) using XPP-AUTO. When the external periodic force is weak ( G = 0 .
05) and thefrequency mismatch is small enough, the system is synchronized to the periodic force with a singlestable amplitude as shown in Fig. 6(b), namely, the amplitude response is monostable. When weapply a stronger force, G = 0 .
4, the region of synchronization becomes wider. The amplitude ofsynchronized oscillations is positive when the frequency mismatch is small, whereas the amplitudeis negative when the mismatch is large. Moreover, there exists a bistable region around r = 0 . R can take either of two stable values, similar to the previous simplermodel with a cubic nonlinearity described by Eq. (47).Figure 6 (d) shows two time courses of DDE (63) with G = 0 . r = 0 . r . Both types of oscillations arestable. In the video in the Supplemental Material [55], the slow convergence of the system stateto either of these two oscillatory states are visualized by projecting the system state onto the( x, dx/dt )-plane. It is noteworthy that the frequency mismatch required for this bistable dynamicsis very small (less than 1 %) in this model. V. SUMMARY
In this study, we have developed a general mathematical framework for reducing delay-inducedlimit-cycle oscillators described by DDEs into a set of nonlinear phase-amplitude equations on thebasis of the Floquet theory. By projecting the original equation onto the reduced phase spacespanned by the first two Floquet eigenfunctions, we derived a set of nonlinear phase-amplitudeequations. We proposed an extended adjoint method for DDEs to numerically calculate the Flo-quet eigenfunctions and their adjoint eigenfunctions. We also developed a method to estimate theasymptotic phase of the system states in a neighborhood of the limit cycle from the phase and am-plitude defined by the Floquet eigenfunctions. The validity of the framework has been confirmedby analyzing two models of delay-induced oscillations. In the present framework, the derivation ofthe reduced equations requires only the calculation of the first two Floquet and adjoint eigenfunc-tions. Therefore, the reduction is practically manageable even though the dynamical system to be23educed is an infinite-dimensional DDE.Despite the simplicity, the resulting reduced equations convey richer information than simplylinearizing the system state around the periodic orbit. To illustrate this, we first studied ananalytically tractable DDE with a cubic nonlinearity. We derived an approximate expression ofthe nonlinear asymptotic phase in terms of the phase and amplitude and verified its validity usingdirect numerical integration of the original system. Moreover, we revealed nontrivial bistablesynchronization of the system with a periodic external forcing, where the amplitude can take twodifferent stable values depending on the initial condition, which cannot be analyzed within theconventional phase-only or the lowest-order phase-amplitude equations. We also analyzed a modelof gene-regulatory oscillator and showed that the reduced phase-amplitude equations also enabledus to capture the nontrivial bistable synchronization with a non-weak periodic force.The result for the gene-regulatory oscillator provides analytical insights into how the weakattraction of the limit cycle and nonlinear interactions between the phase and amplitude can alterthe synchronization dynamics of gene regulatory systems for circadian oscillations. For example,it is known that, in the case of ASPS, out-of-phase (phase-advanced) synchronized oscillation withthe day-and-night lights is stabilized in a similar manner to that is shown in Fig. 6(d) of thesecond model. It has also been reported that the free-running period of circadian oscillation inASPS patients is shorter than 24 h [16], and the temporal therapy (phase advance chronotherapy)can alter the out-of-phase synchronization into in-phase synchronization [17]. Our theoreticalresults imply that weak attraction of the limit cycle and nonlinear interactions between the phaseand amplitude could induce small-amplitude oscillations and bistability of the out-of phase andin-phase synchronized states. If this is the case, the rate of attraction of the system state to thelimit cycle, the Floquet exponent λ in our study, could be used as another effective index tounderstand circadian rhythm disorders in addition to conventional indices like the free-runningperiods and amplitudes of oscillation [11, 14–16]. Thus, the phase-amplitude analysis of delay-induced oscillations developed in this study can shed new light on the complex biological rhythms.There are many other examples of natural and artificial systems that exhibit complex oscillationsdue to the effect of time delay [1–8]. For example, breathing of chronic heart failure patients isa typical example of such natural systems [1]. The present study would provide further insightsinto nontrivial breathing dynamics. An example of artificial systems is the Mackey-Glass electricalcircuit [56] that can be modeled by a DDE, for which the present theory is readily applicableto analyze the synchronization dynamics. The present framework for reducing such time-delayedsystems to a set of nonlinear phase-amplitude equations can be useful as a general analytical24ethod to elucidate the origin of complex synchronization properties under the effect of non-weakperturbations or fluctuations. Further investigation on the nonlinear phase-amplitude equationswould provide us with more insight into the synchronization dynamics in time-delayed systems. VI. ACKNOWLEDGMENTS
We thank Yuki Shimono for helpful advice on computation of the eigenfunctions and G. BardErmentrout for fruitful discussion. This study is supported in part by JSPS KAKENHI (18H04122),The Asahi Glass Foundation, and JST PRESTO (JPMJPR14E2) to KK, and JSPS KAKENHI(JP16K13847, JP17H03279, 18K03471, and JP18H03287) and JST CREST (JPMJCR1913) to HN.
Appendix A: Derivation of the nonlinear phase-amplitude equations
In this section, details of the derivation of the phase-amplitude equations are presented. We firstdefine a phase φ ∈ [0 , T ) along the unperturbed limit cycle of Eq. (1), and represent the T -periodiceigenfunctions q ( t ) j as functions of the phase φ ( t ) as q ( φ ) j , where φ ( t ) = t (mod T ). Because weassume that the functional components associated with the eigenvalues λ i ( i ≥
2) decay quickly,we approximate the system state X ( t ) as X ( t ) ( σ ) ≃ X ( φ )0 ( σ ) + Rq ( φ )1 ( σ ) ( − τ ≤ σ ≤ X ( φ )0 is the system state with phase φ on the limit cycle and R is the amplitude of the eigencomponentcorresponding to λ . Substitution of this approximation into the functional differential equation(1) yields (cid:20) ddφ X ( φ )0 ( σ ) + R ddφ q ( φ )1 ( σ ) (cid:21) ˙ φ + q ( φ )1 ( σ ) ˙ R = ddσ X ( φ )0 ( σ ) + R ddσ q ( φ )1 ( σ ) , ( − τ ≤ σ < N ( X ( φ )0 ( · )) + R ˆ − τ dσ ′ ¯Ω ( φ ) ( σ ′ ) q ( φ )1 ( σ ′ ) + F nl ( φ, R ) + G ( φ, R, t ) , ( σ = 0) (A1)where F nl ( φ, R ) = N (cid:16) X ( φ )0 ( · ) + Rq ( φ )1 ( · ) (cid:17) − N (cid:16) X ( φ )0 ( · ) (cid:17) − R ˆ − τ dσ ′ ¯Ω ( φ ) ( σ ′ ) q ( φ )1 (cid:0) σ ′ (cid:1) . (A2)To derive the phase equation, we project both sides of Eq. (A1) onto the eigenfunction q ( φ )0 .Using the relations ddφ X ( φ )0 ( σ ) = q ( φ )0 ( σ ) , (A3)25 dφ q ( φ )1 ( σ ) = − λ q ( φ )1 ( σ ) + L ( φ ) (cid:16) q ( φ )1 (cid:17) ( σ ) , (A4)and N (cid:16) X ( φ )0 ( · ) (cid:17) = q ( φ )0 (0) , (A5)which follows from the definition q ( φ )0 (0) = dX /dt | t , we obtain h R D q ( φ ) ∗ , L ( φ ) (cid:16) q ( φ )1 (cid:17) ; φ Ei ˙ φ = 1 + R D q ( φ ) ∗ , L ( φ ) (cid:16) q ( φ )1 (cid:17) ; φ E + q ( φ ) ∗ (0) · ( F nl ( φ, R ) + G ( φ, R, t )) . (A6)The phase equation is thus given by˙ φ = 1 + q ( φ ) ∗ (0) · ( F nl ( φ, R ) + G ( φ, R, t ))1 + R D q ( φ ) ∗ , L ( φ ) (cid:16) q ( φ )1 (cid:17) ; φ E = 1 + q ( φ ) ∗ (0) · ( F nl ( φ, R ) + G ( φ, R, t )) − R D q ( φ ) ∗ , L ( φ ) (cid:16) q ( φ )1 (cid:17) ; φ E R D q ( φ ) ∗ , L ( φ ) (cid:16) q ( φ )1 (cid:17) ; φ E q ( φ ) ∗ (0) · ( F nl ( φ, R ) + G ( φ, R, t )) . (A7)Similarly, by projecting both sides of Eq. (A1) onto the eigenfunction q ( θ )1 , we obtain h R D q ( φ ) ∗ , (cid:16) − λ q ( φ )1 + L ( φ ) (cid:16) q ( φ )1 (cid:17)(cid:17) ; φ Ei ˙ φ + ˙ R = R D q ( φ ) ∗ , L ( φ ) (cid:16) q ( φ )1 (cid:17) ; φ E + q ( φ ) ∗ (0) · ( F nl ( φ, R ) + G ( φ, R, t )) . (A8)By substituting Eq. (A7) into Eq. (A8), the amplitude equation is derived as˙ R = λ R + q ( φ ) ∗ (0) · ( F nl ( φ, R ) + G ( φ, R, t )) − R hD q ( φ ) ∗ , L ( φ ) (cid:16) q ( φ )1 (cid:17) ; φ E − λ i R D q ( φ ) ∗ , L ( φ ) (cid:16) q ( φ )1 (cid:17) ; φ E q ( φ ) ∗ (0) · ( F nl ( φ, R ) + G ( φ, R, t )) (A9) Appendix B: Coefficients of the phase-amplitude equations
The expressions for the individual expansion coefficients in the phase and amplitude equa-tions (38) and (39) are as follows. a = 1 T ˆ T D q ( θ ) ∗ , L ( θ ) (cid:16) q ( θ )1 (cid:17) ; θ E dθ, (B1) a = 1 T ˆ T q ( θ ) ∗ (0) · F nl , ( θ ) dθ, (B2)26 = 1 T ˆ T q ( θ ) ∗ (0) · F nl , ( θ ) dθ, (B3) g ( ψ ) = 1 T ˆ T q ( θ ) ∗ (0) · G (cid:18) θ − ψr (cid:19) dθ, (B4) b = 1 T ˆ T D q ( θ ) ∗ , L ( θ ) (cid:16) q ( θ )1 (cid:17) ; θ E dθ, (B5) b = 1 T ˆ T q ( θ ) ∗ (0) · F nl , ( θ ) dθ, (B6) b = 1 T ˆ T q ( θ ) ∗ (0) · F nl , ( θ ) dθ, (B7)and g ( ψ ) = 1 T ˆ T q ( θ ) ∗ (0) · G (cid:18) θ − ψr (cid:19) dθ. (B8) Appendix C: Supplementary video
This supplementary video shows how the oscillators converge to either of the two stable oscilla-tory states. We numerically integrated the DDE (63) of 15 gene-regulatory oscillators subjected toa common sinusoidal external periodic force from several initial conditions from t = 0 to t = 10 .The parameters of the external force are G = 0 . r = 0 . R k = 0 and φ k = k T ( k = 1 , , ..., x, dx/dt )-plane areplotted. In panel (b) of the video, the time courses from two representative initial conditions areshown, where the magenta line is for φ ini = T and the blue one is for φ ini = T .27 x ’ ( t ) (b)(c) (d)(e) (f) tt f f t T T
0 0.5T Tq q q q * * e x p (- l t ) Y q q q q * *
0 0.5T T x ( t )
300 400 500 600−20−15−10−5 time (s) l n ( Y )
400 440−101−2 0 2−202 x(t) (a)(c) q ( f ) ( ) , q ( f ) * ( ) q ( f ) ( ) , q ( f ) * ( ) FIG. 1. (a) Time course of a scalar DDE with a cubic nonlinearity, Eq. (47), showing a slow convergence ofthe system state to the limit cycle. (b) Time course of the oscillator state projected onto the ( x, dx/dt )-plane.[(c), (d)] Extended adjoint method for calculating q ( t )1 ; (c) Peak heights of the time course of Y ( t = nT ) (0)measured at each period vs. t . The red squares are the data from which the Floquet exponent λ is evaluated.(d) Time evolution of exp ( − λ t ) Y ( t ) (0) after compensating the exponential decay. (e) Eigenfunctions andadjoint eigenfunctions associated with λ = 0 plotted as functions of φ . The functions q ( φ )0 (0) and q ( φ ) ∗ (0)are analytically derived, while ¯ q ( φ )0 (0) and ¯ q ( φ ) ∗ (0) are numerically obtained by the extended adjoint method.(f) Eigenfunctions q ( φ )1 (0) and adjoint eigenfunctions q ( φ ) ∗ (0) associated with λ . IG. 2. Evaluation of the asymptotic phase Φ from φ and R . (a) Difference φ − Φ between φ and Φplotted on the ( φ, R )-plane. The data are obtained by direct numerical integration from initial systemstates given by x ( t =0) ( σ ) = x ( φ )0 ( σ ) + Rq ( φ )1 ( σ ). (b) Difference φ − ˆΦ between φ and ˆΦ estimated by usingEq. (61). (c) Absolute difference | Φ − ˆΦ | between the asymptotic phase Φ measured directly by numericalintegration and ˆΦ estimated by using Eq. (61). (d) Asymptotic phase of the initial system states given by x ( t =0) ( σ ) = sin σ + p sin( σ/ φ evaluated using thelinearized isochrons, and the red line indicates the analytical estimation of asymptotic phase Φ. m a x a m p li t ude o f x m a x a m p li t ude o f x w (a) (b) linear phase-amplitude equations DDE linear approximationDDEapproximation r phase-amplitude equations FIG. 3. Maximum amplitude of the DDE (47) subjected to a periodic input. (a) Dependence of maximumamplitude on G at r = 1. (b) Dependence of maximum amplitude on r at G = 0 .
1. Blue points shownumerical results obtained by direct numerical integration of the original DDE (47). The red lines areanalytical predictions by the linear approximation, while black lines are numerical solutions of the phase-amplitude equations (38) and (39). .98 1.02-0.500.5 -3 -2.5 -2 -1.5 -1 -0.5 0-0.200.20.40.60.811.2 R x ( t ) sp1 sp2 w (a) (b)(c) (d) unstablestable w unstablestablesp1sp2 r r RR FIG. 4. (a) Stable and unstable fixed points of the amplitude R plotted against the frequency detuning r at G = 0 .
02. (b) Stable and unstable fixed points at G = 0 .
1. A bistable region exists near r = 1 . ψ , R )-plane at r = 1 .
052 and G = 0 .
1. The gray brokenlines show the nullclines satisfying ˙ ψ = 0 and the gray solid line shows ˙ R = 0. The orange and green dotsshow the stable fixed points (sp1 and sp2) and black lines show the trajectories started from ( − ,
0) and( − . , − . , . − . , − . G ( x, t ) = 0 . . t ). The orange line corresponds to the sp1, while the green line corresponds to thesp2 in panel (c). x ( t ) x ' ( t ) q ( f ) ( ) q ( f ) * ( ) q ( f ) ( ) q ( f ) * ( ) (a) (b)(c) (d) FIG. 5. (a) Time course of the gene-regulatory oscillator Eq. (63) without perturbation ( G = 0), showing aslow convergence to the limit cycle. (b) Trajectory of the system state projected onto the ( x, dx/dt )-plane.(c) Floquet and adjoint eigenfunctions q ( φ )0 (0) and q ( φ ) ∗ (0) associated with λ = 0. (d) Floquet and adjointeigenfunctions q ( φ )1 (0) and q ( φ ) ∗ (0) associated with λ . IG. 6. (a) Asymptotic phase values of initial functions x ( t =0) ( σ ) ≡ p far from the limit cycle. Theblack points indicate the asymptotic phase obtained by direct numerical integration, the blue line indicatesanalytical estimation of the phase ˆ φ , and the red line indicates analytical estimation of the asymptotic phaseˆΦ. (b) Stable and unstable fixed points plotted with respect to r at G = 0 .
05. (c) Fixed points at G = 0 . r = 0 . G = 0 . r = 0 . x ( t =0) ( σ ) = X ( φ ini =0 . T )0 ( σ ), while the blue shows theresult for x ( t =0) ( σ ) = X ( φ ini =0 . T )0 ( σ ) ( − τ ≤ σ ≤
1] L. Glass and M.C. Mackey,
From Clock to Chaos: The Rhythms of Life (Princeton University Press, NewJersey, 1988); B. Balachandran, T. Kalm´ar-Nagy, and D.E. Gilsinn (eds.),
Delay Differential Equations:Recent Advances and New Directions (Springer Verlag, New York, 2009); T. Erneux,
Applied DelayDifferential Equations (Springer, New York, 2009); F. M. Atay (ed.),
Complex Time-Delay Systems (Springer-Verlag, Berlin Heidelberg, 2010).[2] J. Lewis, Current Biology , 1398–1408 (2003); B. Nov´ak and J. J. Tyson, Nature Reviews MolecularCell Biology , 981–991 (2008); L. Herrgen, S. Ares, L.G. Morelli, C. Schr¨oter, F. J¨ulicher, and A.C.Oates, Current Biology , 1244–1253 (2010); M. Das, T. Drake, D. J. Wiley, P. Buchwald, D. Vavy-lonis, and F. Verde, Science , 239–243 (2012); I. Imayoshi, A. Isomura, Y. Harima, K. Kawaguchi,H. Kori, H. Miyachi, T. Fujiwara, F. Ishidate, and R. Kageyama, Science , 1203–1208 (2013).[3] W. Mather, M. R. Bennett, J. Hasty, and L. S. Tsimring, Physical Review Letters , 068105 (2009).[4] B. Doiron, B. Lindner, A. Longtin, L. Maler, and J. Bastian, Physical Review Letters, , 048101(2004); S. A. Campbell, ”Time Delays in Neural Systems”, in Handbook of Brain Connectivity (V. K.Jirsa and A. R. Mclntosh eds.), (Springer, Berlin Heidelberg, 2007); J.W. Kim and P.A. Robinson,Physical Review E , 031907 (2007); I. Yamaguchi, Y. Ogawa, Y. Jimbo, H. Nakao, and K. Kotani,PLoS ONE , e26497 (2011); M. C. Mackey, M. Tyran-Kami´nska, and H.-O. Walther, Journal ofMathematical Biology , 1139–1196 (2017).[5] R. J. Peterka and P. J. Loughlin, Journal of Neurophysiology
1, 410–423 (2004); C. Julien, Cardio-vascular Research , 12–21 (2006); J. J. Batzel and F. Kappel, Mathematical Biosciences , 61–74(2011).[6] M. C. Soriano, J. Garcıa-Ojalvo, C. R. Mirasso, and I. Fischer, Review of Modern Physics , 421–470(2013).[7] T. Kalm´ar-Nagy, G. St´ep´an, and F. C. Moon, Nonlinear Dynamics , 121–142 (2001); S. Chatterjee,Journal of Sound and Vibration , 1860–1876 (2011); E. Stone and S. A. Campbell, Journal ofNonlinear Science , 27–57 (2004).[8] M. Szyd lowski, A. Krawiec, and J. Tobo la, Chaos, Solitons and Fractals , 505–517 (2001); P.Brunovsk´y, A. Erd´elyi, and H.-O. Walther, Journal of Dynamics and Differential Equations , 393–432(2004).[9] J. J. Walker, J. R. Terry, and S. L. Lightman, Proceedings of the Royal Society B: Biological Sciences , 1627–1633 (2010); S. L. Lightman and B. L. Conway-Campbell, Nature Reviews Neuroscience ,710–718 (2010).[10] K. Horikawa, K. Ishimatsu, E. Yoshimoto, S. Kondo, and H. Takeda, Nature , 719–723 (2006).[11] M. A. Lema, D. A. Golombek, and J. Echave, Journal of Theoretical Biology , 565–573 (2000); M.Doi et al. , Nature Communications , 327 (2011); M. Ukai-Tadenuma, R. G. Yamada, H. Xu, J. A.Ripperger, A. C. Liu, and H. R. Ueda, Cell , 268–281 (2011).
12] D. M. Berson, F. A. Dunn, and M. Takao, Science , 1070–1073 (2002).[13] Sack, R. L., Auckley, D., Auger, R. R., Carskadon, M. A., Wright Jr, K. P., Vitiello, M. V., andZhdanova, I. V. Sleep, 30(11), 1460–1483. (2007).[14] H. Ukai, T. J. Kobayashi, M. Nagano, K.-h. Matsumoto, M. Sujino, T. Kondo, K. Yagita, Y. Shigeyoshi,and H. R. Ueda, Nature Cell Biology , 1327–1334 (2007); S. R. Pulivarthy, N. Tanaka, D. K. Welsh,L. De Haro, I. M. Verma, and S. Panda, Proceedings of the National Academy of Sciences of the UnitedStates of America , 20356–20361 (2007).[15] Michael J. Thorpy. Neurotherapeutics. 9(4), 687–701;(2012) Sack, R. L., Auckley, D., Auger, R. R.,Carskadon, M. A., Wright Jr, K. P., Vitiello, M. V., and Zhdanova, I. V. (2007). Sleep, 30(11), 1484-1501.[16] C. R. Jones, S. S. Campbell, S. E. Zone, F. Cooper, A. Desano, P. J. Murphy, B. Jones, L. Czajkowski,and L. J. Ptˇcek, Nature Medicine , 1062–1065 (1999).[17] Moldofsky, H, Musisi, S, and Phillipson, EA. Sleep 1986;9:61-5.[18] A.T. Winfree, The Geometry of Biological Time (Springer, Berlin Heidelberg, 1980).[19] Y. Kuramoto,
Chemical Oscillations, Waves, and Turbulence (Springer, Berlin Heidelberg, 1984).[20] A. Pikovsky, M. Rosenblum, and J. Kurths,
Synchronization: A Universal Concept in Nonlinear Sci-ences (Cambridge University press, Cambridge, 2001).[21] G. B. Ermentrout and D. H. Terman,
Mathematical Foundations of Neuroscience (Springer, New York,2010).[22] H. Nakao, Contemporary Physics , 188–214 (2016).[23] Y. Kuramoto and H. Nakao, ”On the concept of dynamical reduction: the case of coupled oscillators”,Phil. Trans. R. Soc. A , 20190041 (2019).[24] P. Ashwin, S. Coombes, and R. Nicks, Journal of Mathematical Neuroscience , 2 (2016).[25] K. Kotani, I. Yamaguchi, Y. Ogawa, Y. Jimbo, H. Nakao, and G. B. Ermentrout, Physical ReviewLetters , 044101 (2012).[26] V. Noviˇcenko and K. Pyragas, Physica D: Nonlinear Phenomena , 1090-1098 (2012).[27] C. C. Canavier and S. Achuthan, Mathematical Biosciences , 77–96 (2010); S. A. Oprisan and D.I. Austin, PLoS ONE , e0174304 (2017).[28] J. J. Rubin, J. E. Rubin, and G. B. Ermentrout, Physical Review Letters , 204101 (2013); W.Kurebayashi, S. Shirasaka, and H. Nakao, Physical Review Letters , 214101 (2013); K. Pyragasand V. Noviˇcenko, Physical Review E , 012910 (2015); Y. Park and G. B. Ermentrout, Journalof Computational Neuroscience , 269–281 (2016); V. Klinshov, S. Yanchuk, A. Stephan, and V.Nekorkin, Europhysics Letters , 50006 (2017).[29] O. Castej´on, A. Guillamon, and G. Huguet, Journal of Mathematical Neuroscience , 13 (2013); D.Wilson and J. Moehlis, Physical Review E , 052213 (2016); S. Shirasaka, W. Kurebayashi, and H.Nakao, Chaos , 023119 (2017); D. Wilson and G. B. Ermentrout, Journal of Mathematical Biology, , 37-66 (2018):
30] A. Stokes, Proceedings of the National Academy of Sciences of the United States of America ,1330-1334 (1962).[31] J. K. Hale and S. M. V. Lunel, Introduction to Functional Differential Equations (Springer-Verlag, NewYork, 1993).[32] C. Simmendinger, A. Wunderlin, and A. Pelster, Physical Review E , 5344 (1999).[33] W. Wischert, A. Wunderlin, A. Pelster, M. Olivier, and J. Groslambert, Physical Review E, , 203(1994).[34] G. Palm and T. Poggio, SIAM Journal on Applied Mathematics , 195–216 (1977).[35] J. K. Hale, Ordinary Differential Equations (Wiley-Interscience, New York, 1969).[36] A. Rabinovitch and M. Friedman, Physics Letters A , 319–325 (2006); K. Yoshimura and K. Arai,Physical Review Letters , 154101 (2008); H. Ben Amor, N. Glade, C. Lobos, and J. Demongeot,Acta Biotheoretica , 121–142 (2010); D. S. Goldobin, J.-N. Teramae, H. Nakao, and G. B. Ermen-trout, Physical Review Letters , 154101 (2010); H. M. Osinga and J. Moehlis, SIAM Journal onApplied Dynamical Systems , 1201-1228 (2010); G. Huguet and R. de la Llave, SIAM Journal onApplied Dynamical Systems , 1763-1802 (2013); A. Mauroy, B. Rhoads, J. Moehlis, and I. Mezic,SIAM Journal on Applied Dynamical Systems , 306–338 (2014); J. Hannam, B. Krauskopf and H.M. Osinga, The European Physical Journal Special Topics , 2645–2654 (2016); M. Bonnin, Inter-national Journal of Circuit Theory and Applications , 636–659 (2017).[37] K. C. Wedgwood, K. K. Lin, R. Thul, and S. Coombes, Journal of Mathematical Neuroscience , 2(2013).[38] G. S. Medvedev, Journal of Nonlinear Science , 441–464 (2011).[39] O. Diekmann, S. A. van Gils, S. M. V. Lunel, and H.-O. Walther, Delay Equations: Functional-,Complex-, and Nonlinear Analysis (Springer, New York, 1995).[40] The dual space of C is usually defined as a vector space consisting of all bounded linear functionalswith domain C . However, it is difficult to work directly with this space, because, in general, theadjoint system on this space is not described by a differential equation, i.e., the adjoint semigroup isnot strongly continuous [39, Ch. 2]. Hale introduced a bilinear form, under which the adjoint semigroupon the associated adjoint space is strongly continuous. In the celebrated work of Hale and Lunel, theconvenient system is called not the adjoint system but the transposed system to distinguish betweenthem. The transposed system can be viewed as the ”adjoint” system since their eigenfunctions areclosely related [31, Chs. 7 and 8]. Hence, we call the transposed system as the adjoint system in thispaper as in the literature [32].[41] B. Ermentrout, Y. Park, and D. Wilson, ”Recent advances in coupled oscillator theory”, Phil. Trans.R. Soc. A , 20190092 (2019).[42] The frequency of this subtraction procedure may not necessarily be once per time T . Generally speaking,the faster the q component grows relative to q , the more frequent subtraction would be required. Inthis study, at every 10 oscillation periods, whether Y ( t ) decays exponentially or not is judged by easuring 10 peaks of Y ( nT ) . If this is the case, the Floquet eigenvalue λ is evaluated using the peaks(Fig. 1(c)). The waveform of q ( t )1 ( σ ) is obtained by averaging q ( t )1 ( σ ) = e − λ t · Y ( t ) ( σ ) over the 10oscillation periods, where exponential decay is compensated (Fig. 1(d), (e)).[43] The eigenspace associated with the complex conjugate eigenvalues are numerically obtained by theconsecutive subtraction approach proposed in Section II B, but some care is required when evaluatingthe pair of complex conjugate eigenfunctions. Ding and Cvitanovi´c have developed an algorithm tocalculate the complex Floquet eigenvalues and eigenfunctions [44]. By introducing a moving coordinateframe spanned by the complex eigenfunctions as in Section II C, a reduced ODE form may be obtainedeven in such a situation. Though not generic, Floquet eigenvalues with multiplicities can also arise,whose treatment is beyond the scope of this paper. We here only mention that, in some mathematicalliterature, spectral projections of the time- T fundamental solution of linear periodic DDEs onto gener-alized eigenspaces associated with multiple eigenvalues has been discussed on the basis of the Dunfordintegral using residue calculus [39, 45].[44] X. Ding and P. Cvitanovi´c, SIAM Journal on Applied Dynamical Systems , 1434–1454 (2016).[45] M. V. S. Frasson and S. M. V. Lunel, Integral Equations and Operator Theory , 91–121 (2003).[46] H. Kori and Y. Morita, Dynamical System Approach to Biological Rhythms (Kyoritsu Syuppan, Tokyo,2011). (in Japanese)[47] A. Stokes, Journal of Differential Equations , 153–172 (1977).[48] J. K. Hale and M. Weedermann, Journal of Differential Equations , 219–246 (2004).[49] I. Le´on, and D. Paz´o, Physical Review E, , 012211 (2019).[50] M. Rosenblum and A. Pikovsky, Phil. Trans. R. Soc. A. , 20190093 (2019).[51] For simplicity, we use Y ( t =0)ini ( σ ) ≡ p when we naively evaluate it basedon an argument introduced on a plane ( x, ˙ x ) as tan − ( x/ ˙ x ).[53] B. Ermentrout, Simulating, Analyzing, and Animating Dynamical Systems: A Guide to XPPAUT forResearchers and Students (Society for Industrial and Applied Mathematics, Philadelphia, 2002).[54] The origin of the phase is set at t that satisfies x ( t ) = h x ( t ) i with dx ( t ) /dt >0.[55] See Supplemental Material at http://link.aps.org/supplemental/10.1103/PhysRevResearch.2.033106for a video that shows how the gene-regulatory oscillators subjected to a periodic external force convergeto either of the two stable oscillatory states.[56] Namaj¯unas, A., Pyragas, K., and Tamaˇseviˇcius, A. Physics Letters A, 201(1), 42-46 (1995).