Nonlinear phase-amplitude reduction of delay-induced oscillations

Spontaneous oscillations induced by time delays are observed in many real-world systems. Phase reduction theory for limit-cycle oscillators described by delay-differential equations (DDEs) has been developed to analyze their synchronization properties, but it is applicable only when the perturbation applied to the oscillator is sufficiently weak. In this study, we formulate a nonlinear phase-amplitude reduction theory for limit-cycle oscillators described by DDEs on the basis of the Floquet theorem, which is applicable when the oscillator is subjected to perturbations of moderate intensity. We propose a numerical method to evaluate the leading Floquet eigenvalues, eigenfunctions, and adjoint eigenfunctions necessary for the reduction and derive a set of low-dimensional nonlinear phase-amplitude equations approximately describing the oscillator dynamics. By analyzing an analytically tractable oscillator model with a cubic nonlinearity, we show that the asymptotic phase of the oscillator state in an infinite-dimensional state space can be approximately evaluated and non-trivial bistability of the oscillation amplitude caused by moderately strong periodic perturbations can 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 complex dynamics under non-weak perturbations, which may be relevant to real physiological phenomena such as circadian rhythm sleep disorders.


I. INTRODUCTION
Time-delayed feedback can break continuous time-translational symmetry and lead to oscillatory behavior in many physical, biological, social, and engineered systems [1][2][3][4][5][6][7][8]. In biology, for example, ultradian oscillations in the hypothalamic-pituitary-adrenal (HPA) system are induced by timedelayed synthesis of hormones in the adrenal cortex [9]. Also, somite segmentation in zebrafish is regulated by oscillatory dynamics induced by time delays in the synthesis of proteins [10], and mammalian circadian rhythm is generated by feedback regulations of clock genes in suprachiasmatic nucleus (SCN) [11]. Such oscillatory dynamics in systems with time delays can be described as stable limit-cycle orbits of delay-differential equations (DDEs).
In many of such systems, each oscillatory unit, or oscillator, is not isolated but perturbed by external forcing or by mutual coupling with other oscillators, and the state of each oscillator may deviate from the unperturbed limit cycle of an isolated oscillator when the perturbation is not sufficiently weak. Therefore, it is important to understand how perturbations of moderate intensity 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 couplings between cells are necessary for the spatio-temporal oscillatory dynamics of her1 (zebrafish hairyrelated gene1) expression [10]. In the case of circadian clock genes, the oscillatory period in the freerunning condition is known to be slightly different from 24 h, but they are entrained by the periodic external day-and-night lights through retinal ganglion cells [12,13]. Strong light stimulation can further induce large modulation in the activities of the clock genes [14]. Since irregular dynamics of circadian rhythms manifest themselves as diseases such as sleep disorders [13,[15][16][17], understanding of the dynamics of circadian clock genes under strong perturbations may facilitate therapies for sleep disorders.
The phase reduction theory is a standard mathematical framework for characterizing response properties of weakly perturbed limit-cycle oscillators and analyzing their synchronization dynamics via dimensionality reduction [18][19][20][21][22][23][24]. Recently, the phase reduction theory has been extended also to DDEs exhibiting stable limit-cycle oscillations, which requires non-trivial mathematical generalization because DDEs are infinite-dimensional dynamical systems [25,26]. However, the phase reduction has a strong limitation in that it is applicable only when the oscillator state remains sufficiently close to the unperturbed limit cycle. Specifically, when non-weak perturbations are applied or relaxation time of the system state to the limit cycle is not sufficiently small, the amplitude degrees of freedom may no longer be enslaved by the phase, leading to the breakdown of the lowest-order phase-only description. In such cases, the nonlinear interaction of the phase and 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 oscillatory systems described by ordinary differential equations (ODEs), such as higher-order phase resetting curves [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 their infinite-dimensional nature. Thus, a general framework for dimensionality reduction of limit-cycle oscillators described by DDEs that can analyze the effect of moderately strong perturbations is needed. 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 interact significantly in a nonlinear manner. We develop a nonlinear phase-amplitude reduction theory for DDEs, which gives a general mathematical framework for reducing DDEs describing limit-cycle oscillators to low-dimensional ordinary differential equations on the basis of the Floquet theory [30][31][32]. We also propose a practical numerical method, which we call the extended adjoint method, to evaluate the Floquet eigenvalues, eigenfunctions, and their adjoint functions, which are necessary for the reduction. By using biorthogonality of the Floquet eigenfunctions and their adjoints, we project the oscillator state onto an eigenspace spanned by a few slowly-decaying Floquet eigenfunctions and derive a set of phase-amplitude equations which takes into account nonlinear interactions between the slowly-decaying Floquet eigenmodes. In contrast to the standard lowest-order phase reduction, 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 is strongly perturbed.
We confirm the validity of the theory using an analytically-tractable DDE with a cubic nonlinearity by showing that the reduced phase-amplitude equations accurately predict the amplitude of the phase-locked oscillations under a periodic force, which exhibits non-trivial bistable response induced by the non-weak amplitude effects. We then apply the theory to a model of a generegulatory oscillator under moderately strong forcing and analyze its synchronization dynamics.
We show that the reduced phase-amplitude equations can also predict nontrivial bistable dynamics of the system, which is analogous to a circadian disorder called advanced sleep-phase syndrome (ASPS) [15][16][17].

II. THEORY
In this section, we derive a set of reduced nonlinear phase-amplitude equations for limit-cycle oscillators described by DDEs on the basis of the Floquet theory and propose a practical numerical method to calculate the Floquet eigenvalues, eigenfunctions and their adjoints that are necessary for the reduction. We also derive approximate phase-amplitude equations for the oscillators subjected to 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 when they are periodically perturbed, is not easy because they describe nonlinear infinite-dimensional dynamical systems on Banach spaces. Our aim is to derive simpler tractable equations by reducing them to low-dimensional ODEs while preserving their essential quantitative properties and to analyze synchronization dynamics of nonlinear oscillators described by such DDEs under moderately strong external perturbations. In previous studies [25,26], phase reduction methods for stable limit-cycle solutions of DDEs have been developed, which are applicable when the perturbations given to the system is sufficiently weak. In this study, we develop a nonlinear phase-amplitude reduction theory for DDEs.
We consider a DDE for X(t) ∈ R N , represented as a column vector, with a maximum delay time τ > 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, [30][31][32]. Here, where || · || is the usual Euclidean norm on R N . This history function X (t) represents the state of the dynamical system described by the DDE at time t, where the state space of the system is given by the infinite-dimensional Banach space C 0 .
Using the above notation, a DDE can generally be written as Here, the vector-valued functional N : C 0 → R N represents the system dynamics and G : C 0 × R → 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 systems with discrete delays but also systems with distributed delays [33]; see Ref. [34] for the relation between nonlinear functionals and their kernel representations, which are widely used for systems with distributed delays described by integro-differential equations.
We consider a situation in which the DDE (1) without the external perturbation (G = 0) has a stable limit-cycle solution X 0 (t) whose period is T , i.e., X 0 (t + T ) = X 0 (t), and represent it as a history function X In what follows, we also denote the limit cycle as X (φ) 0 , where we use the phase φ (0 ≤ φ < T ) in place of the time t to parametrize system state on the limit cycle. The phase φ increases from 0 to T , where the origin φ = 0 can be chosen as a specific system state on the limit cycle. When the system state evolves along the limit cycle without perturbation, the phase φ increases with a constant 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 limit cycle by assigning the same phase value φ to the set of system states {X (t) } that asymptotically converge to the same system state as X (φ) 0 when the system evolves without perturbation [25,26], i.e., lim t→∞ X (t) − X (φ+t) 0 C 0 = 0, yielding the notion of asymptotic phase Φ(X (t) ) : C 0 → [0, T ) that maps a system state X (t) in the basin to a phase value. The asymptotic phase Φ satisfies when the system state evolves in the basin of the limit cycle without perturbation. The isosurfaces of Φ, called isochrons, are not simply hyperplanes in general. For ordinary differential equations, the asymptotic phase has been used as a canonical representation of rhythms of stable oscillatory dynamics [18][19][20][21]35] and provides in-depth insights into strongly-perturbed oscillatory dynamics [24,[27][28][29]36]. Recently, it has also been defined for DDEs and other non-conventional oscillatory systems [25,26].
We assume that the relaxation dynamics of the system state to the limit cycle can be decomposed into a few slow modes and remaining faster modes, which are well separated in time scale from each other. In this case, a rectangular coordinate frame moving along the periodic orbit, which was used in Refs. [35,37,38], is not useful for reducing the dynamics to low-dimensional ODEs, because fast and slow components interact already at the lowest order in this coordinate frame. It is also not easy to proceed with the asymptotic phase and associated amplitudes, because they are generally given by highly nonlinear functionals of the system state X(t). We therefore use a coordinate frame defined by the Floquet eigenfunctions to decompose the system state as discussed in Ref. [23] for ODEs. The space spanned by the Floquet eigenfunctions with non-vanishing relaxation rates is tangent to the isochron at each point on the limit cycle. For this purpose, we need to calculate the Floquet 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 0 (t) as Y (t) = X(t) − X 0 (t), and introduce its The linearized variational equation for Y (t) is given by where L (t) (Y (t) ) is a history representation of a linear functional defined by Here,Ω 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 In what follows, we expand N in a functional Taylor series in Y (t) as where L (t) Y (t) (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 these functionals are evaluated at X (t) = X (t) 0 .
As an example, let us consider a simple DDE, which is equivalent to in the history-function representation. By using the chain rule for functional differentiation and representing the terms in N as and using Dirac's delta function δ(·), we obtain where 0 (−τ )). The linearized dynamics for the deviation Y (t) can then be written as Let us introduce a time-periodic linear operatorL of period T , which acts on a complexified and rewrite Eq. (5) as (LY (t) )(σ) = 0. Because Y (t) obeys a periodically driven linear system, by the Floquet theorem for linear DDEs [30][31][32], the spectrum ofL is at most countable and is satisfied, where λ i ∈ C is the i-th Floquet eigenvalue and q (t) i ∈ (C 0 ) C is the corresponding Tperiodic Floquet eigenfunction (i = 0, 1, 2, ...). Here, the largest eigenvalue, which is 0 and simple by the Floquet theorem, is denoted as λ 0 = 0 and the other eigenvalues are arranged in descending order of the real part.
We also introduce adjoint eigenfunctions with respect to a bilinear form appropriate for DDEs [40]. Following Refs. [30][31][32], we define a bilinear form of two functions, A ∈ (C 0 ) C and B ∈ (C 0 ) * C , as Here, is the dual space of (C 0 ) C with respect to the bilinear form, consisting of (row) vector-valued functions that map the interval [0, τ ] to C N * , and [·, ·] denotes the Hermitian scalar product of V ∈ C N * and U ∈ C N defined as [V, U ] = N k=1 V k U k where V k and U k are vector components of V and U , respectively. An adjoint operatorL * ofL with respect to this bilinear form can then be derived as where Here, Y * (t) ∈ C N * is a row vector of N complex components and with respect to the bilinear form Eq. (17), and hence they can be normalized to satisfy the biorthogonal relation The zero eigenfunction of the linear operatorL can be chosen as q , which can be confirmed by differentiating Eq. (1) with respect to t at 0 on the periodic orbit [25]. Note that this definition specifies the normalization of q  i (i = 1, 2, ...), we normalize them such that max 0≤t≤T q (t) i (0) = 1. We use this convention for the normalization throughout this study. We note that the zero eigenfunction of the linear operatorL corresponds to the tangential component along the limit cycle, namely, the phase direction. Moreover, the zero eigenfunction q (t) * 0 of the adjoint operatorL * gives the phase sensitivity function of the limit cycle, which characterizes linear response property of the oscillator phase to weak perturbations [25,26]. Similarly, the other eigenfunctions q (t) * i (i = 1, 2, ...) characterize linear response properties of the amplitudes and called isostable response curves for the case of ODEs [41].
The adjoint eigenfunctions can numerically be obtained by an extension of the adjoint method for DDEs, which was previously used to calculate the adjoint zero eigenfunction ofL [25,26]. That is, we numerically integrate the linearized and its adjoint equations while subtracting unnecessary functional components by using the biorthogonality between the eigenfunctions and adjoint eigenfunctions. The main difference from the adjoint method for q (t) * 0 developed in the previous studies is that we calculate the adjoint eigenfunctions also for λ i (i ≥ 1). Therefore, during numerical integration, we need to remove unnecessary functional components in the directions of the lowerorder eigenfunctions from 0-th to (i − 1)-th, which grow faster than the i-th component in order to calculate the i-th eigenfunction precisely. For i ≥ 1, we also need to renormalize the solutions of the equations by a factor e λ i t determined by the Floquet exponent in order to obtain the correct eigenfunctions.
To numerically calculate the i-th eigenfunction q forward in time. During the calculation, we subtract the 0-th to (i−1)-th eigencomponents from the numerical 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 . See Sec. III.C and Ref. [42] for further details. In a similar way, the i-th adjoint eigenfunction q (t) * i is calculated by numerically integrating the adjoint linear equation backward in time while subtracting unnecessary eigencomponents and then compensating the numerical 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 by projecting the system state onto a moving coordinate frame spanned by a small number of Floquet eigenfunctions. That is, we decompose the deviation of the system state X (t) from that on the limit cycle X (t) 0 by using the eigenfunctions associated with the leading M eigenvalues other than 0, which are assumed to be real and simple for the sake of simplicity [43], as where X (φ) 0 is a system state on the limit cycle parametrized by the phase φ ∈ [0, T ), q The symbol ≃ indicates that we approximate X (t) (σ) by its projection on the space spanned by 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 state along the Floquet eigenfunction corresponding to the direction transversal to the limit cycle and represents the deviation of the system state from the limit cycle. Here, the phase value φ for a given state X (t) is determined in such a way that the state difference along the limit cycle. Thus, we assume the following orthogonality condition: is on the hyperplane that is tangent to the isochron on the limit cycle 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, when the eigenvalue λ 1 with the largest non-zero part is close to 0, the perturbed system state does not go back to the limit cycle quickly, and hence nonlinear interactions between the phase eigenmode and the slowest-decaying amplitude eigenmode should be taken into account for better description of the system.
For ordinary differential equations, such coupled nonlinear phase-amplitude equations have been derived 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 the treatments of oscillatory dynamics in these studies are rather abstract. Quantitative analysis of synchronization dynamics of DDEs using the coordinate transform proposed therein have not been very fruitful despite their potential advantages, mainly due to the lack of practical methods for numerically evaluating the Floquet eigenfunctions.
We hereafter restrict ourselves to the case in which λ 1 takes a negative real value near zero and Re{λ 2 } ≪ λ 1 for simplicity. To derive the phase and amplitude equations, we retain only the slowest two modes associated with λ 0 and λ 1 and approximate X (t) (σ) as where R(t) = ρ 1 (t) is the amplitude of the eigenmode corresponding to λ 1 . The symbol ≃ here indicates that we further approximate X (t) (σ) by its projection on a one-dimensional space spanned by q (φ) 1 . We substitute this expression into Eq. (1) and then project both sides of Eq. (1) onto the eigenfunctions q (φ) 0 and q (φ) 1 , respectively, by using biorthogonality of the eigenfunctions and derive the equations for the phase φ and the amplitude R.
As explained in Appendix A, the phase equation can be derived as or, by rewriting the right-hand side, and the amplitude equation can similarly be derived as where the nonlinear functional N in Eq. (8) is approximated by an ordinary function of φ and R, (29) and the external perturbation is also approximated as 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 2 ) or higher, because the constant terms and linear terms in R have already been subtracted in Eq. (29).
Thus, by projecting the DDE onto the first two eigenfunctions, a set of two-dimensional coupled ordinary differential equations for the phase φ and amplitude R is obtained. In order to consider the higher-order effects of the amplitude deviations, we have not expanded the third-order terms in Eq. (27) and Eq. (28) in a series of R and hence the dynamics of φ and R are nonlinearly coupled at the second and higher orders in R. This nonlinearity can be a source of intriguing oscillatory dynamics [27][28][29]36]. We also note that the lowest-order phase-amplitude equations (see Refs. [23,29,41] for the case of ODEs) are obtained at the lowest-order approximation in R, where F nl (φ, R) is O(R 2 ) and does not appear at the lowest order.
Finally, before proceeding, we note that there are also other formulations of phase or phaseamplitude reduced equations for analyzing higher-order effects of perturbations on limit cycles described 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 frequency is close to the natural frequency of the oscillator, we may further derive simpler, approximate phase-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 and periodic in t with period T ′ = T /r (frequency r), i.e., We also assume that the detuning between the natural frequency of the oscillator and the periodic force is small and denote it as ∆ω = 1 − r.
We introduce a slow phase variable ψ ≡ φ − rt. The equations for ψ and R can be written as and We also expand the nonlinear term F nl in Taylor series of R up to R N as where {F nl,ℓ } (ℓ = 2, 3, ...) are expansion coefficients. Note that the series for F nl starts from O(R 2 ).
Considering that ψ evolves only slowly while rt rapidly increases, we approximate the terms with ψ + rt in Eqs. (33) and (34) by their one-period average, for example, as and q (ψ+rt) * 0 where ψ is kept constant during the integration. Expanding F nl (ψ + rt, R) up to O(R 3 ) and averaging the coefficients, we obtain approximate equations for ψ and R as and where the equations for the individual coefficients are given in Appendix B. We check the validity of 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 choosing the phase value φ so that it satisfies the orthogonality condition Eq. (24). Numerically, we find the valueφ that minimizes the mean squared error, There exists a neighborhood of the periodic orbit where the phase and amplitude components defined by using the Floquet eigenfunctions are uniquely determined [48,Lemma1]. However, in general, there can exist multiple values ofφ that satisfy Eq. (24) in the range 0 ≤φ < T . To choose the appropriate value from them, for each candidate ofφ, we evaluate the corresponding q 1 component asR and adopt the value ofφ that has the smallest R as the estimate of φ, and the smallestR as the estimate of R.

F. Approximate evaluation of the asymptotic phase
The phase φ defined by the Floquet eigenfunction, which we use in the present study for the phase-amplitude description, is different from the asymptotic phase Φ; the isosurface of Φ is generally curved and tangent to the isophase plane of φ at each point on the limit cycle. Since the asymptotic phase Φ provides useful information on the nonlinear dynamical properties of the oscillator, it is convenient if we can approximate Φ using φ and R. In this subsection, we propose a 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 where The asymptotic phase Φ of the system state X (t 0 ) at time t 0 with phase φ 0 and amplitude R 0 can approximately be obtained by integrating d(φ(s), R(s)) until the system state goes back sufficiently close to the limit cycle as When 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 with rate λ 1 as at the lowest-order approximation. The asymptotic phase Φ of the system state X (t 0 ) can then be approximately evaluated aŝ In Sec. III E and Sec. IV, we use the above method to estimate the asymptotic phase Φ of the oscillator 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 oscillator described by a scalar DDE with a cubic nonlinearity, for which approximate expressions of the Floquet eigenfunctions and their adjoints can be analytically derived, and analyze the effect of a periodic force on the dynamics.

A. Model
The model is represented as where x(t) ∈ R, ǫ = 0.05 is a small constant, and the external periodic force is described by where G 0 is the intensity of the periodic force and r is the ratio of the natural frequency 2π/T of the 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 0 (t) = sin t, or 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 large as compared to the oscillation period as shown in Figs. 1(a) and (b).
We denote the small deviation of the system state from the limit cycle as . The linear operatorL of this system is given by Eq. (6) with where δ is Dirac's delta function. By retaining the first two leading eigenvalues, the nonlinear phase-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 of the model Eq. (47) without the external force (G = 0) analytically. In what follows, we consider the case in which the relaxation of the system state to the limit cycle is slow and assume that λ 1 is small and O(ǫ). First, the zero eigenfunction ofL is given exactly as and the adjoint eigenfunction is To find the exponent λ 1 with the second largest real part, we introduce an ansatz where l is a constant and plug this into Eqs. (5) and (6). We then obtain the approximate eigenvalue and the associated eigenfunction up to O(ǫ) as and respectively. Similarly, for the corresponding adjoint eigenfunction, we approximately obtain sin(t + s) + π 2 cos(t + s) where the constant C * 1 is determined from the normalization condition q 1 ; t = 1.

C. Numerical evaluation of the eigenvalues and eigenfunctions
To confirm the validity of the approximate analytical results for the Floquet eigenvalues, eigenfunctions, and adjoint eigenfunctions obtained in the previous subsection, we numerically evaluate these quantities by the extended adjoint method and compare with the approximate analytical results.
First, as in the conventional adjoint method for DDEs [25,26], we compute q (t) 0 (σ) (−τ ≤ σ ≤ 0), which is simply dX 0 /dt| t+σ , and then q     (0) are plotted with respect to φ in Fig. 1(f). We can confirm a good agreement between the numerical results and approximate analytical results for the eigenfunctions. The numerical value of the largest negative exponent λ 1 is approximately evaluated as −0.030, which is also close to the theoretical value −8ǫ/(π 2 + 4) = −0.029.

D. Phase-amplitude equations
We now derive a set of nonlinear phase-amplitude equations from Eq. (47) with the periodic sinusoidal force. The nonlinear term F nl (φ, R) in Eq. (29) is explicitly given by 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 approximate equations for the phase difference ψ = φ − rt and R are given in the form of Eqs. (38) and (39) with and Using numerically evaluated eigenvalues and eigenfunctions, the coefficients in Eqs. (38) and (39) can be calculated as λ 1 = −0.029, a 0 = 1.8418,a 2 = 0.0436, a 3 = 0.0415; b 0 = 1.5353,b 2 = −0.0053, b 3 = 0.0212; and g 01 = −0.9622, g 02 = 0, g 11 = −0.8239, and g 12 = 0.5245. From these coefficients, the equations for the phase difference ψ and the amplitude R are obtained as Thus, we have approximately reduced an infinite-dimensional dynamical system described by a DDE 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 phase derived in Sec. II F by evolving the present model from initial conditions far from the limit cycle. We also directly evaluate the asymptotic phase Φ for several initial conditions by numerically integrating the system and measuring the time necessary for the system state to converge sufficiently close to the limit cycle for comparison.
As the first example, we try to estimate Φ when the initial function is on the φ-R plane, that Figure 2 (a) shows Φ − φ for given initial values of φ and R obtained by direct numerical integration of the DDE, and Fig. 2 (b) shows analytical results ofΦ − φ obtained from Eq. (61). Figure 2 (c) shows the absolute difference between Φ and its analytical estimationΦ. We can confirm a good agreement between the approximate analytical curve 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 the initial functions as x (t=0) (σ) = sin σ + p sin(σ/2) (−τ ≤ σ ≤ 0) with varying values of p [52], and evaluated their asymptotic phase Φ by direct numerical integration of the DDE. Figure 2 (d) shows the phase φ, the asymptotic phase Φ estimated by Eq. (61), and the asymptotic phase Φ obtained by direct numerical integration. We can confirm that the approximate analytical estimate of the asymptotic phase given by Eq. (61) gives reasonable agreement with the direct numerical results even 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 with small frequency detuning. In particular, we focus on the average effect of the periodic force on the amplitude R in the phase-locked state, which cannot be analyzed without the amplitude equation.
Since g 02 = 0 in Eq. (58), the ψ-nullcline on whichψ = 0 is obtained from the averaged equation (38) as where F s represents the right-hand side of Eq. (39). The effect of the intensity of the periodic force G 0 and the detuning ∆ω on the stationary amplitude R of the oscillation in the steady state can be evaluated from the partial derivatives of F s (R, G 0 , ∆ω) by the implicit function theorem. Figure 3 shows the predicted amplitude of the oscillation. The dependence of the amplitude on G 0 at r = 1 is plotted in Fig. 3(a), where the stationary amplitude obtained from the averaged phase-amplitude equations (38) and (39) are compared with the linear approximation of the stationary amplitude with a slope ∂R/∂G 0 | R=0,G 0 =0,∆ω=0 = 18.2. Similarly, Fig. 3 (b) shows the dependence of the amplitude on r at G 0 = 0.1, where the result of the phase-amplitude equations are compared with linear approximation of the amplitude with a slope ∂R/∂r | R=0.72,G 0 =0.1,∆ω=0 = 9.91.
We can confirm that the linear approximation appropriately predicts the changes in the stationary amplitude of the delay-induced oscillator subjected to a non-weak external periodic force when it is slightly modulated. Moreover, the nonlinear phase-amplitude equations can predict the amplitude more precisely than the linear approximation in the given parameter range.

G. Bistable response of delay-induced oscillation to a periodic force
In this subsection, we demonstrate that the present model can exhibit a nontrivial bistable response to a periodic force by a bifurcation analysis of Eq. (38) and Eq. (39). Such a phenomenon results from higher-order amplitude effects and cannot be described by the phase-only equation nor the lowest-order phase-amplitude equations. Using XPP-AUTO [53], we numerically find stationary solutions in the range R > −0.5 where the inverse 1/(1 + Ra 0 ) exists (note that a 0 = 1.8418).
Depending on the parameters G 0 and r, we observe quantitatively different behaviors of the system state as shown in Fig. 4.  (1.052t). We can clearly observe the bistable dynamics of the oscillator caused by moderately strong periodic forcing.

IV. GENE-REGULATORY OSCILLATOR
In this section, as a more complex, biologically-motivated example of DDEs, we investigate a model of gene regulation [3] under a periodic sinusoidal force given by where x(t) ∈ R is the state variable representing protein concentration and α, β, γ, C 0 , R 0 , and the delay time τ are real parameters. The first term of the right-hand side represents protein synthesis with time delay for transcription and translation, while the second and the third terms represent degradation and dilution of the protein, respectively. Following the previous research [3], we set β = 0.1, C 0 = 10, and τ = 1. The external periodic force is G(t) = G 0 sin 2π T rt with intensity G 0 and frequency mismatch r. We set the rate constant of synthesis as α = 100, degradation as γ = 100, and Michaelis constant of degradation as R 0 = 10 so that the system exhibits a slow convergence 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 only numerically. Figures 5 shows the system state x (t) converging toward the limit-cycle attractor; For this model, the T -periodic linear operatorL is given by Eq. (15) with We first evaluate the validity of the approximate expression of the asymptotic phase in the absence of the external force (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 integration of 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 seen that the approximate analytical result reproduces the result of direct numerical measurement of the asymptotic phase.
We next consider how the gene-regulatory oscillator behaves when it is subjected to a periodic external force. We conduct bifurcation analysis for different values of G 0 and r in the same way as that for Eq. (47) using XPP-AUTO. When the external periodic force is weak (G 0 = 0.05) and the frequency mismatch is small enough, the system is synchronized to the periodic force with a single stable amplitude as shown in Fig. 6(b), namely, the amplitude response is monostable. When we apply a stronger force,  tions. Therefore, the reduction is practically manageable even though the dynamical system to be reduced is an infinite-dimensional DDE.
Despite the simplicity, the resulting reduced equations convey richer information than simply linearizing the system state around the periodic orbit. To illustrate this, we first studied an analytically tractable DDE with a cubic nonlinearity. We derived an approximate expression of the nonlinear asymptotic phase in terms of the phase and amplitude and verified its validity using direct numerical integration of the original system. Moreover, we revealed nontrivial bistable synchronization of the system with a periodic external forcing, where the amplitude can take two different stable values depending on the initial condition, which cannot be analyzed within the conventional phase-only or the lowest-order phase-amplitude equations. We also analyzed a model of gene-regulatory oscillator and showed that the reduced phase-amplitude equations also enabled us 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 weak attraction of the limit cycle and nonlinear interactions between the phase and amplitude can alter the 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 with the day-and-night lights is stabilized in a similar manner to that is shown in Fig. 6(d) of the second model. It has also been reported that the free-running period of circadian oscillation in ASPS 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 theoretical results imply that weak attraction of the limit cycle and nonlinear interactions between the phase and amplitude could induce small-amplitude oscillations and bistability of the out-of phase and in-phase synchronized states. If this is the case, the rate of attraction of the system state to the limit cycle, the Floquet exponent λ 1 in our study, could be used as another effective index to understand circadian rhythm disorders in addition to conventional indices like the free-running periods and amplitudes of oscillation [11,[14][15][16]. Thus, the phase-amplitude analysis of delayinduced 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 oscillations due to the effect of time delay [1][2][3][4][5][6][7][8]. For example, breathing of chronic heart failure patients is a typical example of such natural systems [1]. The present study would provide further insights into nontrivial breathing dynamics. An example of artificial systems is the Mackey-Glass electrical circuit [56] that can be modeled by a DDE, for which the present theory is readily applicable to analyze the synchronization dynamics. The present framework for reducing such time-delayed systems to a set of nonlinear phase-amplitude equations can be useful as a general analytical method to elucidate the origin of complex synchronization properties under the effect of non-weak perturbations or fluctuations. Further investigation on the nonlinear phase-amplitude equations would 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. Bard j , where φ (t) = t (mod T ). Because we assume that the functional components associated with the eigenvalues λ i (i ≥ 2) decay quickly, we approximate the system state X (t) as is the system state with phase φ on the limit cycle and R is the amplitude of the eigencomponent corresponding to λ 1 . Substitution of this approximation into the functional differential equation To derive the phase equation, we project both sides of Eq. (A1) onto the eigenfunction q (φ) 0 . Using the relations and N X which follows from the definition q (φ) 0 (0) = dX 0 /dt| t , we obtain The phase equation is thus given bẏ Similarly, by projecting both sides of Eq. (A1) onto the eigenfunction q By substituting Eq. (A7) into Eq. (A8), the amplitude equation is derived aṡ and     (σ) (−τ ≤ σ ≤ 0). The green line represents the external force. measuring 10 peaks of Y (nT ) . If this is the case, the Floquet eigenvalue λ 1 is evaluated using the peaks ( Fig. 1(c)). The waveform of q (t) 1 (σ) is obtained by averaging q (t) 1 (σ) = e −λ1t · Y (t) (σ) over the 10 oscillation periods, where exponential decay is compensated (Fig. 1(d), (e)).
[43] The eigenspace associated with the complex conjugate eigenvalues are numerically obtained by the consecutive subtraction approach proposed in Section II B, but some care is required when evaluating the pair of complex conjugate eigenfunctions. Ding and Cvitanović have developed an algorithm to calculate the complex Floquet eigenvalues and eigenfunctions [44]. By introducing a moving coordinate frame spanned by the complex eigenfunctions as in Section II C, a reduced ODE form may be obtained even 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 mathematical literature, spectral projections of the time-T fundamental solution of linear periodic DDEs onto generalized eigenspaces associated with multiple eigenvalues has been discussed on the basis of the Dunford integral using residue calculus [39,45].