Nonlinear phase-amplitude reduction of delay-induced oscillations

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 Kotani ; Pyragas , 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)\in\mathbb{R}^{N}$, represented as a column vector, with a maximum delay time $\tau>0$. To construct a solution of the DDE, we have to take into account the history of $X(t)$ from $t-\tau$ to $t$. Thus, we introduce its history-function representation, $X^{(t)}({\sigma})\equiv X(t+{\sigma})$ $(-\tau\leq\sigma\leq 0)$ Stokes ; Hale ; Simmendinger . Here, $X^{(t)}(\cdot)\in C_{0}$ and $C_{0}=C([-\tau,0]\to\mathbb{R}^{N})$ is a Banach space of (column) vector-valued continuous functions mapping $[-\tau,0]$ into $\mathbb{R}^{N}$, which is equipped with a norm $||x||_{C_{0}}=\sup_{\theta\in[-\tau,0]}||x(\theta)||$, where $||\cdot||$ is the usual Euclidean norm on $\mathbb{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 $\mathcal{N}:C_{0}\to{\mathbb{R}}^{N}$ represents the system dynamics and $G:C_{0}\times\mathbb{R}\to{\mathbb{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 Wischert ; see Ref. Palm 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 (4) 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_{0}^{(t)}(\cdot)\in C_{0}$ satisfying $X_{0}^{(t+T)}=X_{0}^{(t)}$, where

In what follows, we also denote the limit cycle as $X_{0}^{(\phi)}$, where we use the phase $\phi$ ($0\leq\phi<T$) in place of the time $t$ to parametrize system state on the limit cycle. The phase $\phi$ increases from $0$ to $T$, where the origin $\phi=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 $\phi$ increases with a constant frequency $1$, i.e., $\phi=t\ (\mbox{mod}\ T)$. Similarly, we also denote $T$-periodic history functions, such as the Floquet eigenfunctions, using the phase $\phi$ 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 $\phi$ to the set of system states $\{X^{(t)}\}$ that asymptotically converge to the same system state as $X_{0}^{(\phi)}$ when the system evolves without perturbation Kotani ; Pyragas , i.e., $\lim_{t\to\infty}\|X^{(t)}-X_{0}^{(\phi+t)}\|_{C_{0}}=0$, yielding the notion of asymptotic phase $\Phi(X^{(t)}):C_{0}\to[0,T)$ that maps a system state $X^{(t)}$ in the basin to a phase value. The asymptotic phase $\Phi$ satisfies

when the system state evolves in the basin of the limit cycle without perturbation. The isosurfaces of $\Phi$, 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 Winfree2 ; Kuramoto ; Pikovsky ; Ermentrout2 ; Hale_ODE and provides in-depth insights into strongly-perturbed oscillatory dynamics Ashwin ; Canavier+ ; Rubin+ ; Castejon+ ; Rabinovitch+ . Recently, it has also been defined for DDEs and other non-conventional oscillatory systems Kotani ; Pyragas .

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. Hale_ODE ; Wedgwood ; Medvedev , 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. KuramotoPTRSA 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.

We first describe the Floquet theory for the DDE (4) 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 history-function representation $Y^{(t)}(\cdot)\in C_{0}$ with $Y^{(t)}({\sigma})\equiv Y(t+{\sigma})$ $(-\tau\leq\sigma\leq 0)$ as

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 $\mathcal{N}$ with respect to $X^{(t)}$ evaluated at the system state $X^{(t)}=X_{0}^{(t)}$ on the limit cycle. Note that Eq. (8) gives a periodically driven linear system because $X_{0}^{(t)}$ is $T$-periodic. In what follows, we expand $\mathcal{N}$ in a functional Taylor series in $Y^{(t)}$ as

where $L^{(t)}\left(Y^{(t)}\right)(0)$ represents a linear functional of $Y^{(t)}$ defined in Eq. (12) with $\sigma=0$ and $F_{\mathrm{nl}}(Y^{(t)}(\cdot))$ represents the remaining nonlinear functional of $Y^{(t)}$, respectively, and both of these functionals are evaluated at $X^{(t)}=X_{0}^{(t)}$.

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 ${\mathcal{N}}$ as

and

using Dirac’s delta function $\delta(\cdot)$, we obtain

where $\mathcal{N}_{j}(t)\equiv\partial_{x_{j}}\mathcal{N}(x_{1},x_{2})$ $(j=1,2)$ is evaluated at $(x_{1},x_{2})=(X_{0}^{(t)}(0),X_{0}^{(t)}(-\tau))$. The linearized dynamics for the deviation $Y(t)$ can then be written as

Let us introduce a time-periodic linear operator $\hat{L}$ of period $T$, which acts on a complexified Banach space $(C_{0})_{\mathbb{C}}$ (Diekmann, , Sec. III.7) as

and rewrite Eq. (8) as $(\hat{L}Y^{(t)})(\sigma)=0$. Because $Y^{(t)}$ obeys a periodically driven linear system, by the Floquet theorem for linear DDEs Stokes ; Hale ; Simmendinger , the spectrum of $\hat{L}$ is at most countable and

is satisfied, where $\lambda_{i}\in{\mathbb{C}}$ is the $i$-th Floquet eigenvalue and $q_{i}^{(t)}\in(C_{0})_{\mathbb{C}}$ is the corresponding $T$-periodic Floquet eigenfunction ($i=0,1,2,...$). Here, the largest eigenvalue, which is $0$ and simple by the Floquet theorem, is denoted as $\lambda_{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 bilinear . Following Refs. Stokes ; Hale ; Simmendinger , we define a bilinear form of two functions, $A\in(C_{0})_{\mathbb{C}}$ and $B\in(C_{0})_{\mathbb{C}}^{*}$, as

Here, $(C_{0})_{\mathbb{C}}^{*}=C([0,\tau]\to\mathbb{C}^{N*})$ is the dual space of $(C_{0})_{\mathbb{C}}$ with respect to the bilinear form, consisting of (row) vector-valued functions that map the interval $[0,\tau]$ to $\mathbb{C}^{N*}$, and $[\cdot,\cdot]$ denotes the Hermitian scalar product of $V\in\mathbb{C}^{N*}$ and $U\in\mathbb{C}^{N}$ defined as $[V,U]=\sum_{k=1}^{N}V_{k}\overline{U_{k}}$ where $V_{k}$ and $U_{k}$ are vector components of $V$ and $U$, respectively. An adjoint operator $\hat{L}^{*}$ of $\hat{L}$ with respect to this bilinear form can then be derived as

where

Here, $Y^{*}(t)\ \in\mathbb{C}^{N*}$ is a row vector of $N$ complex components and $Y^{(t)*}(s)\equiv Y^{*}(t+s)$ $(0\leq s\leq\tau)\in(C_{0})_{\mathbb{C}}^{*}$ is its history-function representation.

The adjoint eigenfunction $q_{i}^{(t)*}\in(C_{0})^{*}_{\mathbb{C}}$ of $q_{i}^{(t)}$, which is also $T$-periodic, satisfies

where $\bar{\lambda}_{i}$ is the complex conjugate of ${\lambda}_{i}$. If $\lambda_{i}\neq\lambda_{j}$, $q_{i}^{(t)}$ is orthogonal to $q_{j}^{(t)*}$ with respect to the bilinear form Eq. (29), and hence they can be normalized to satisfy the biorthogonal relation $\langle{q_{i}^{(t)*}},q_{j}^{(t)};t\rangle=\delta_{i,j}$. The zero eigenfunction of the linear operator $\hat{L}$ can be chosen as $q_{0}^{(t)}(\sigma)=dX_{0}/dt|_{t+\sigma}$ ($-\tau\leq\sigma\leq 0$), which can be confirmed by differentiating Eq. (4) with respect to $t$ at $X^{(t)}=X_{0}^{(t)}$ on the periodic orbit Kotani . Note that this definition specifies the normalization of $q_{0}^{(t)}$. For the other eigenfunctions $q_{i}^{(t)}$ ($i=1,2,...$), we normalize them such that $\max_{0\leq t\leq T}\left(q_{i}^{(t)}(0)\right)=1$. We use this convention for the normalization throughout this study. We note that the zero eigenfunction of the linear operator $\hat{L}$ corresponds to the tangential component along the limit cycle, namely, the phase direction. Moreover, the zero eigenfunction $q_{0}^{(t)*}$ of the adjoint operator $\hat{L}^{*}$ gives the phase sensitivity function of the limit cycle, which characterizes linear response property of the oscillator phase to weak perturbations Kotani ; Pyragas . Similarly, the other eigenfunctions $q_{i}^{(t)*}$ ($i=1,2,...$) characterize linear response properties of the amplitudes and called isostable response curves for the case of ODEs ErmentroutPTRSA .

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 of $\hat{L}$ Kotani ; Pyragas . 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_{0}^{(t)*}$ developed in the previous studies is that we calculate the adjoint eigenfunctions also for $\lambda_{i}\>(i\geq 1)$. Therefore, during numerical integration, we need to remove unnecessary functional components in the directions of the lower-order 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\geq 1$, we also need to renormalize the solutions of the equations by a factor $e^{\lambda_{i}t}$ determined by the Floquet exponent in order to obtain the correct eigenfunctions.

To numerically calculate the $i$-th eigenfunction $q_{i}^{(t)}$, we integrate the linearized equation

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 $\lambda_{i}$ is numerically evaluated from the exponential decay rate of $Y^{(t)}$. Then the eigenfunction $q_{i}^{\left(t\right)}\left(\sigma\right)$ is obtained by compensating the exponential decay of $Y^{(t)}(\sigma)$ as $q_{i}^{\left(t\right)}\left(\sigma\right)=e^{-\lambda_{i}t}\ Y^{(t)}(\sigma)$ ($-\tau\leq\sigma\leq 0$). See Sec. III.C and Ref. footnotesec3c for further details. In a similar way, the $i$-th adjoint eigenfunction $q_{i}^{(t)^{*}}$ 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^{-\lambda_{i}t}$. We call this procedure the extended adjoint method in this study.

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_{0}^{(t)}$ 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 lambda_complex , as

where $X_{0}^{(\phi)}$ is a system state on the limit cycle parametrized by the phase $\phi\in\left[0,T\right)$, $q_{i}^{(\phi)}$ ($i=1,...,M$) is the Floquet eigenfunction associated with $\lambda_{i}$ and denoted as a function of $\phi$ rather than $t$, and $\{\rho_{i}(t)\}$ are real expansion coefficients representing amplitudes of the Floquet eigenmodes. The symbol $\simeq$ indicates that we approximate $X^{(t)}(\sigma)$ by its projection on the space spanned by the $M$ eigenfunctions $\{q_{1}^{(\phi)},...,q_{M}^{(\phi)}\}$. 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 $\phi$ for a given state $X^{(t)}$ is determined in such a way that the state difference $X^{(t)}-X_{0}^{(\phi)}$ does not have a tangential functional component $q_{0}^{(\phi)}$ along the limit cycle. Thus, we assume the following orthogonality condition:

namely, the difference $X^{(t)}-X_{0}^{(\phi)}$ is on the hyperplane that is tangent to the isochron on the limit cycle at $X_{0}^{(\phi)}$. 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}^{(\phi)}\}$ ($i=1,...,M$), nonlinear interactions between different eigenmodes generally arise. Specifically, when the eigenvalue $\lambda_{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 Wedgwood ; Morita . Such transformation methods have also been developed for DDEs in Refs. Stokes2 ; Hale2 , 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 $\lambda_{1}$ takes a negative real value near zero and $\mbox{Re}\{\lambda_{2}\}\ll\lambda_{1}$ for simplicity. To derive the phase and amplitude equations, we retain only the slowest two modes associated with $\lambda_{0}$ and $\lambda_{1}$ and approximate $X^{(t)}({\sigma})$ as

where $R(t)=\rho_{1}(t)$ is the amplitude of the eigenmode corresponding to $\lambda_{1}$. The symbol $\simeq$ here indicates that we further approximate $X^{(t)}(\sigma)$ by its projection on a one-dimensional space spanned by $q_{1}^{(\phi)}$. We substitute this expression into Eq. (4) and then project both sides of Eq. (4) onto the eigenfunctions $q_{0}^{(\phi)}$ and $q_{1}^{(\phi)}$, respectively, by using biorthogonality of the eigenfunctions and derive the equations for the phase $\phi$ 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 $\mathcal{N}$ in Eq. (14) is approximated by an ordinary function of $\phi$ and $R$,

and the external perturbation is also approximated as

In Eq. (43) and Eq. (45), both the second and third terms on the right-hand side depend on $F_{\mathrm{nl}}$ and $G$. Note that $F_{\mathrm{nl}}(\phi,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. (46).

Thus, by projecting the DDE onto the first two eigenfunctions, a set of two-dimensional coupled ordinary differential equations for the phase $\phi$ 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. (43) and Eq. (45) in a series of $R$ and hence the dynamics of $\phi$ and $R$ are nonlinearly coupled at the second and higher orders in $R$. This nonlinearity can be a source of intriguing oscillatory dynamics Canavier+ ; Rubin+ ; Castejon+ ; Rabinovitch+ . We also note that the lowest-order phase-amplitude equations (see Refs. Castejon+ ; ErmentroutPTRSA ; KuramotoPTRSA for the case of ODEs)

are obtained at the lowest-order approximation in $R$, where $F_{\mathrm{nl}}\left(\phi,R\right)$ 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 phase-amplitude reduced equations for analyzing higher-order effects of perturbations on limit cycles described by ODEs, such as non-pairwise phase interactions KuramotoPTRSA , higher-order phase reduction Pazo , nonlinear phase coupling function RosenblumPTRSA , and higher-order approximations of coupling functions ErmentroutPTRSA , which can capture more detailed aspects of synchronization than the lowest-order phase equation.

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^{\prime}=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 $\Delta\omega=1-r$.

We introduce a slow phase variable $\psi\equiv\phi-rt$. The equations for $\psi$ and $R$ can be written as

and

We also expand the nonlinear term $F_{\mathrm{nl}}$ in Taylor series of $R$ up to $R^{N}$ as

where $\{F_{\mathrm{nl},\ell}\}$ ($\ell=2,3,...$) are expansion coefficients. Note that the series for $F_{\mathrm{nl}}$ starts from $O(R^{2})$.

Considering that $\psi$ evolves only slowly while $rt$ rapidly increases, we approximate the terms with $\psi+rt$ in Eqs. (51) and (54) by their one-period average, for example, as

and

where $\psi$ is kept constant during the integration. Expanding $F_{\mathrm{nl}}(\psi+rt,R)$ up to $O(R^{3})$ and averaging the coefficients, we obtain approximate equations for $\psi$ 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.

Numerically, the values of the phase $\phi$ 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 $\phi$ so that it satisfies the orthogonality condition Eq. (39). Numerically, we find the value $\hat{\phi}$ 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 (Hale2, , Lemma1). However, in general, there can exist multiple values of $\hat{\phi}$ that satisfy Eq. (39) in the range $0\leq\hat{\phi}<T$. To choose the appropriate value from them, for each candidate of $\hat{\phi}$, we evaluate the corresponding $q_{1}$ component as

and adopt the value of $\hat{\phi}$ that has the smallest $\left|\hat{R}\right|$ as the estimate of $\phi$, and the smallest $\hat{R}$ as the estimate of $R$.

The phase $\phi$ defined by the Floquet eigenfunction, which we use in the present study for the phase-amplitude description, is different from the asymptotic phase $\Phi$; the isosurface of $\Phi$ is generally curved and tangent to the isophase plane of $\phi$ at each point on the limit cycle. Since the asymptotic phase $\Phi$ provides useful information on the nonlinear dynamical properties of the oscillator, it is convenient if we can approximate $\Phi$ using $\phi$ and $R$. In this subsection, we propose a method to approximately evaluate the asymptotic phase of an unperturbed oscillator from $\phi$ and $R$ defined by the Floquet eigenfunctions, which is valid when $R$ is sufficiently small.

When the perturbation is absent ($G=0$), Eq. (41) for $\phi$ can be written as

where

The asymptotic phase $\Phi$ of the system state $X^{(t_{0})}$ at time $t_{0}$ with phase $\phi_{0}$ and amplitude $R_{0}$ can approximately be obtained by integrating $d(\phi(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 $\phi$ and $R$ and assume that $\phi$ increases constantly with frequency $1$ and $R$ decays exponentially with rate $\lambda_{1}$ as

at the lowest-order approximation. The asymptotic phase $\Phi$ of the system state $X^{(t_{0})}$ can then be approximately evaluated as