Phase-amplitude reduction of transient dynamics far from attractors for limit-cycling systems

TL;DR

This paper develops a phase-amplitude reduction framework for analyzing transient dynamics of limit-cycling systems far from attractors, integrating isostables with Koopman operator analysis for improved accuracy and applicability.

Contribution

It introduces a fully consistent reduction method combining isostables and Koopman analysis, enabling transient dynamics study far from the limit cycle.

Findings

Avoids discontinuities in isostables

Enables analysis of states far from the limit cycle

Demonstrates efficient suppression of deviations in a biochemical oscillator

Abstract

Phase reduction framework for limit-cycling systems based on isochrons has been used as a powerful tool for analyzing rhythmic phenomena. Recently, the notion of isostables, which complements the isochrons by characterizing amplitudes of the system state, i.e., deviations from the limit-cycle attractor, has been introduced to describe transient dynamics around the limit cycle [Wilson and Moehlis, Phys. Rev. E 94, 052213 (2016)]. In this study, we introduce a framework for a reduced phase-amplitude description of transient dynamics of stable limit-cycling systems. In contrast to the preceding study, the isostables are treated in a fully consistent way with the Koopman operator analysis, which enables us to avoid discontinuities of the isostables and to apply the framework to system states far from the limit cycle. We also propose a new, convenient bi-orthogonalization method to obtain the response functions of the amplitudes, which can be interpreted as an extension of the adjoint covariant Lyapunov vector to transient dynamics in limit-cycling systems. We illustrate the utility of the proposed reduction framework by estimating optimal injection timing of external input that efficiently suppresses deviations of the system state from the limit cycle in a model of a biochemical oscillator.

Equations

$$\dot{\mathbf{X}}=\mathbf{F}(\mathbf{X}),\mathbf{X}\in {\mathbb{R}}^N,$$

$$U^t{s}_i(\mathbf{X})=e^{{\lambda }_{i}t}{s}_i(\mathbf{X}),$$

$$U^{\Delta t}{r}_i(\mathbf{X})=\mathrm{Re}(s_i(\varphi (\Delta t,\mathbf{X})))=e^{{\lambda }_i\Delta t}{r}_i(\mathbf{X}),$$

$$\dot{r_i}(\mathbf{X})=\underset{\Delta t\to 0}{\lim }\frac{U^{\Delta t}{r}_i(\mathbf{X})-r_i(\mathbf{X})}{\Delta t}={\lambda }_i{r}_i.$$

$$\dot{\theta }(\mathbf{X})=\omega .$$

$$\dot{\mathbf{X}}=\mathbf{F}(\mathbf{X})+ϵ\mathbf{p}(t).$$

$$\dot{\theta }$$

$$\dot{r_i}$$

$$\dot{\theta }$$

$$\dot{r_i}$$

$$U^T{r}_i({\mathbf{h}}_p)=e^{{\lambda }_{i}T}{r}_i(\mathbf{h}({\Theta }_1)+\delta \mathbf{x}).$$

$$U^T{r}_i({\mathbf{h}}_p)=r_i(\mathbf{h}({\Theta }_1)+\mathbf{M}(\mathbf{h}({\Theta }_1))\delta \mathbf{x}+O(||\delta \mathbf{x}|{|}^2)).$$

$$\nabla r_i^{†}(\mathbf{h}({\Theta }_1))\mathbf{M}(\mathbf{h}({\Theta }_1))=e^{{\lambda }_{i}T}\nabla r_i^{†}(\mathbf{h}({\Theta }_1)).$$

$$\nabla {\theta }^{†}(\mathbf{h}({\Theta }_1))\mathbf{M}(\mathbf{h}({\Theta }_1))=\nabla {\theta }^{†}(\mathbf{h}({\Theta }_1)).$$

$$\nabla r_i({\mathbf{X}}^{\ast }(t))\cdot \mathbf{F}({\mathbf{X}}^{\ast }(t))={\lambda }_i{r}_i,$$

$$\nabla \theta ({\mathbf{X}}^{\ast }(t))\cdot \mathbf{F}({\mathbf{X}}^{\ast }(t))=\omega .$$

$$\frac{\mathrm{d}\nabla {\psi }_i({\mathbf{X}}^{\ast }(t))}{\mathrm{d}t}\cdot \delta \mathbf{x}(t)$$

$$=-\nabla {\psi }_i({\mathbf{X}}^{\ast }(t))\cdot \mathrm{D}\mathbf{F}({\mathbf{X}}^{\ast }(t))\delta \mathbf{x}(t)$$

$$=-\mathrm{D}{\mathbf{F}}^{†}({\mathbf{X}}^{\ast }(t))\nabla {\psi }_i({\mathbf{X}}^{\ast }(t))\cdot \delta \mathbf{x}(t).$$

$$\frac{\mathrm{d}\nabla {\psi }_i({\mathbf{X}}^{\ast }(t))}{\mathrm{d}t}=-\mathrm{D}{\mathbf{F}}^{†}({\mathbf{X}}^{\ast }(t))\nabla {\psi }_i({\mathbf{X}}^{\ast }(t))$$

$${\mathbf{\gamma }}_i({\mathbf{X}}^{\ast }(t))\cdot \nabla {\psi }_j({\mathbf{X}}^{\ast }(t))={\delta }_{ij},$$

$$({\mathbf{\gamma }}_i({\mathbf{X}}^{\ast }(t))\cdot \mathbf{z}(t))\nabla {\psi }_i({\mathbf{X}}^{\ast }(t)).$$

$${\dot{\mathbf{\gamma }}}_i({\mathbf{X}}^{\ast }(t))=D\mathbf{F}({\mathbf{X}}^{\ast }(t)){\mathbf{\gamma }}_i({\mathbf{X}}^{\ast }(t)).$$

$$(\nabla {\psi }_i({\mathbf{X}}^{\ast }(t))\cdot \mathbf{y}(t)){\mathbf{\gamma }}_i({\mathbf{X}}^{\ast }(t)),$$

$$\dot{x}=\frac{\alpha }{1+z^n}-x,$$

$$\dot{y}=x-y,$$

$$\dot{z}=y-z.$$

$$s_{\ast }={\mathrm{argmin}}_{s\in {\mathcal{I}}_{\sigma }}|r_2({\mathbf{X}}_p^{\ast }(T_e))|,$$

$$\mathrm{sgn}(r_2({\mathbf{X}}^{\ast }(0))){\int }_0^{T_e}\mathbf{p}(t)\cdot \nabla r_2({\mathbf{X}}^{\ast }(t))e^{{\lambda }_2(T_e-t)}\mathrm{d}t$$

$$f_{{\lambda }_1}^{\ast }(\mathbf{X})=\underset{s\to \infty }{\lim }\frac{1}{s}{\int }_0^{s}f\circ \varphi (t,\mathbf{X})e^{-{\lambda }_{1}t}\mathrm{d}t.$$

$$f_{{\lambda }_i}^{\ast }(\mathbf{X})=\underset{s\to \infty }{\lim }\frac{1}{s}{\int }_0^s\left[f\circ \varphi (t,\mathbf{X})-\bar{f}-\sum _{k=1}^{i-1}{f}_{{\lambda }_k}^{\ast }(\mathbf{X})e^{{\lambda }_{k}t}\right]e^{-{\lambda }_{i}t}\mathrm{d}t,$$

$$g_i(\mathbf{X})=\nabla r_i({\mathbf{X}}_0({\theta }_{\ast }))\cdot (\mathbf{X}-{\mathbf{X}}_0({\theta }_{\ast })),$$

$$(g_i{)}_{{\lambda }_i}^{\ast }(\mathbf{X})=\underset{s\to \infty }{\lim }\frac{1}{s}{\int }_0^s{g}_i\circ \varphi (t,\mathbf{X})e^{-{\lambda }_{i}t}\mathrm{d}t.$$

Full text

The phase reduction theory provides a general framework to simplify a complex, multi-dimensional limit-cycling system describing a stable rhythmic activity to a one-dimensional phase equation evolving on a circle Winfree (2001); Kuramoto (2011); Hoppensteadt and Izhikevich (1997); Ermentrout and Terman (2010); Nakao (2016); Ashwin, Coombes, and Nicks (2016). It has been successfully used to understand synchronization phenomena of weakly interacting rhythmic elements in physical, chemical, biological and engineered systems Winfree (2001); Kuramoto (2011); Hoppensteadt and Izhikevich (1997); Ermentrout and Terman (2010); Nakao (2016); Ashwin, Coombes, and Nicks (2016); Pikovsky, Rosenblum, and Kurths (2003); Strogatz et al. (2005); Tinsley, Nkomo, and Showalter (2012); Tass (2007); Schultheiss, Prinz, and Butera (2011); Dörfler, Chertkov, and Bullo (2013). Methods to optimize and control synchronization of rhythmic elements have also been developed by using the phase reduction framework Moehlis, Shea-Brown, and Rabitz (2006); Kiss et al. (2007); Harada et al. (2010); Zlotnik et al. (2013, 2016). However, to describe the system dynamics far from the limit cycle, amplitude degrees of freedom should be taken into account. In this study, by extending preceding studies, we propose a phase-amplitude reduction framework that is applicable to transient dynamics far from the limit cycle.

I Introduction

The roles of amplitude degrees of freedom in limit-cycling systems, which represent deviations of the system states from the limit-cycle attractor and are eliminated in the phase-reduction framework, have been extensively studied because they are rich sources of intriguing oscillator dynamics at individual Pikovsky, Rosenblum, and Kurths (2003); Yoshimura and Arai (2008); Goldobin et al. (2010); Wedgwood et al. (2013); Hitczenko and Medvedev (2013); Mauroy et al. (2014); Ashwin, Coombes, and Nicks (2016) and ensemble Kuramoto (2011); Matthews, Mirollo, and Strogatz (1991); Cross and Hohenberg (1993); Aranson and Kramer (2002); Pikovsky, Rosenblum, and Kurths (2003); Nakao and Mikhailov (2009); Koseska, Volkov, and Kurths (2013); Ashwin, Coombes, and Nicks (2016) levels. In most studies, however, the analysis is restricted to the vicinity of a supercritical Hopf bifurcation, where a simple normal form (Stuart-Landau equation) of the oscillator dynamics is available Guckenheimer and Holmes (1983); Strogatz (2014). Some other studies use moving orthonormal frames along the limit cycle to define the amplitudes of the oscillator Wedgwood et al. (2013); Hitczenko and Medvedev (2013); Ashwin, Coombes, and Nicks (2016), which allow the quantitative study of the amplitude dynamics of oscillators far from bifurcation points. However, in general, those amplitude variables interact nonlinearly with each other, which hinders simplification of the system description. Thus, it is highly desirable to establish a framework for a quantitative reduced description of limit-cycling systems applicable to transient dynamics far from the limit cycle. Such a framework would facilitate in-depth studies of the roles of amplitude degrees of freedom of limit-cycling systems in realistic settings.

The key idea in the phase reduction is assigning the same phase value to the set of initial conditions that share the same asymptotic behavior. These sets of identical phase values are called isochrons Winfree (2001); Kuramoto (2011); Hoppensteadt and Izhikevich (1997); Ermentrout and Terman (2010); Nakao (2016); Pikovsky, Rosenblum, and Kurths (2003). Analogously, in a recent work Mauroy, Mezić, and Moehlis (2013), the notion of isostables is introduced by identifying the initial conditions that share the same relaxation property, i.e., the same decay rate toward the attractor. It has also been shown Mauroy, Mezić, and Moehlis (2013) that the isochrons and isostables can be understood from a unified point of view of the spectral properties of the Koopman (composition) operator Budišić, Mohr, and Mezić (2012). For each characteristic decay rate of the system state toward the attractor, a set of isostables representing an amplitude degree of freedom can be introduced, which is independent from the phase and the other amplitude degrees of freedom. By retaining a small number of amplitude variables representing dominant (slowly-decaying) part of the transient dynamics, reduced description of the system dynamics can be derived. The Koopman operator has attracted broad interest recently, because it is closely related to a rapidly developing data-driven approach to complex nonlinear systems, called the dynamic mode decomposition Budišić, Mohr, and Mezić (2012); Rowley et al. (2009); Schmid (2010); Williams, Kevrekidis, and Rowley (2015); N. B. Erichson and Kutz (2016); Mann and Kutz (2016); Proctor, Brunton, and Kutz (2016).

Amplitude reduction frameworks for a system near a stable equilibrium based on isostables have been established for multi-dimensional Mauroy, Mezić, and Moehlis (2013); Mauroy (2014); Wilson and Moehlis (2015) and infinite-dimensional systems Wilson and Moehlis (2016a) and have been used to formulate optimal control problems of moving the system state toward the equilibrium Mauroy, Mezić, and Moehlis (2013); Wilson and Moehlis (2015, 2016a). Recently, Wilson and Moehlis Wilson and Moehlis (2016b) have extended the isostable reduction framework to limit-cycling systems. However, the isostables introduced in their work have discontinuities on one leaf of the isochrons. To avoid this problem, it is assumed in Ref. Wilson and Moehlis (2016b) that the system evolves in a close-enough neighborhood of the limit cycle so that the discontinuities are negligible, and the amplitude response to perturbation in their reduced system involves the first order response evaluated only on the limit cycle. Therefore, their analysis is essentially equivalent to deriving a decoupled linear system preserving spectral properties of the original system in a vicinity of the limit-cycle attractor (called kinematically similar system in terms of Lyapunov transformations Lyapunov (1992); Adrianova (1995); Colonius and Kliemann (2014)) by making use of covariant properties of adjoint covariant Lyapunov vectors Kuptsov and Parlitz (2012) (also called adjoint Floquet vectors Traversa et al. (2015) or dual Lyapunov vectors Pikovsky and Politi (2015)). A method to analyze response functions of decoupled phase and amplitude variables in limit-cycling systems, which is based on the Lie symmetries formalism and is valid far from the attractors, has also been proposed Guillamon and Huguet (2009); Castejón, Guillamon, and Huguet (2013). However, the latter analysis is limited to two-dimensional dynamical systems and naive application of the method proposed in Ref. Castejón, Guillamon, and Huguet (2013), that is, solving adjoint equations to calculate the response functions, can yield flawed results numerically, as we discuss in this paper.

In this study, we introduce a phase-amplitude reduction framework to describe transient dynamics of stable limit-cycle oscillators, which is applicable to high-dimensional dynamics far from the limit-cycle attractor. We propose a systematic bi-orthogonalization method to numerically estimate the fundamental quantities for the reduction accurately, i.e., the first order response functions of the phase and amplitudes to perturbations along a given trajectory, which is not necessarily the limit cycle itself. These response functions can be interpreted as an extension of the adjoint covariant Lyapunov vectors to transient dynamics. We illustrate the utility of the proposed framework by estimating optimal injection timing of external input that realizes maximal suppression of the most persistent (least decaying) amplitude degree of freedom.

This paper is organized as follows: in Sec. II, phase and amplitudes in limit-cycling systems are introduced using the Koopman operator theory. In Sec. III, the phase-amplitude reduction framework for limit-cycling systems is introduced and the bi-orthogonalization method to obtain their response properties is developed. In Sec. IV, the theory is illustrated by analyzing the phase-amplitude response properties of a minimal chemical kinetic model of an oscillatory genetic circuit. Also, the optimal injection timing problem is introduced and analyzed. Section V summarizes the results.

II Phase, amplitudes and the Koopman operator

We consider a $N$-dimensional autonomous dynamical system

[TABLE]

where ${\bm{X}}(t)$ is a system state and ${\bm{F}}({\bm{X}})$ is a vector field. Suppose the system (1) has a periodic orbit $\chi:{\bm{X}}_{0}(t)$ with period $T$. Let $\phi:\mathbb{R}\times\mathbb{R}^{N}\to\mathbb{R}^{N}$ denote the flow induced by Eq. (1), i.e., $\phi(t,{\bm{X}})$ is the solution of Eq. (1) at the time $t$ with the initial condition ${\bm{X}}$ at $t=0$.

The stability of the periodic orbit $\chi$ is characterized by the characteristic multipliers Guckenheimer and Holmes (1983) $\Lambda_{i}\ (i=1,\cdots,N)$, which are the eigenvalues of the time-$T$ flow linearized around a point ${\bm{X}}_{0}(t_{*})$ on the orbit $\chi$ (also called the monodromy matrix): ${\bf M}({\bm{X}}_{0}(t_{*}))=\partial\phi(T,{\bm{X}})/\partial{\bm{X}}|_{{\bm{X}}={\bm{X}}_{0}(t_{*})}$. When the relation $1=\Lambda_{1}>|\Lambda_{2}|\geq\cdots\geq|\Lambda_{N}|$ holds, the periodic orbit $\chi$ is a stable limit cycle. For simplicity, we hereafter assume that the Floquet multipliers $\Lambda_{i}$ are positive, real, and simple. Extension to the case with complex conjugate multipliers can be performed in a parallel way to the analysis of stable equilibria Mauroy, Mezić, and Moehlis (2013); Wilson and Moehlis (2016a). We consider dynamics of the system in the basin of attraction ${\mathcal{B}}\subset\mathbb{R}^{N}$ of the stable limit cycle $\chi$.

The Koopman operator $U^{t}$ is a linear operator that describes the evolution of a function defined on the phase space, called an observable $f:\mathbb{R}^{N}\to\mathbb{C}$. It is defined as $U^{t}f({\bm{X}})=f\circ\phi(t,{\bm{X}})$, where $\circ$ represents composition of functions. The operator $U^{t}$ has eigenfunctions Mezić (2013); Lan and Mezić (2013) $s_{i}({\bm{X}})\ (i=1,\cdots,N)$ associated with eigenvalues $\lambda_{i}\ (i=1,\cdots,N)$, that is,

[TABLE]

where $\lambda_{1}=\sqrt{-1}\omega$, $\omega\equiv 2\pi/T$, and $\lambda_{i}=\mathrm{log}(\Lambda_{i})/T\ (i=2,\cdots,N)$. The eigenvalues correspond to the characteristic exponents of the limit cycle $\chi$ Guckenheimer and Holmes (1983), hence they reflect the spectral property of the limit-cycling system.

We hereafter assume that the vector field ${\bm{F}}$ is twice continuously differentiable so that the continuously differentiable eigenfunctions $s_{i}$ exist on the whole basin of attraction Lan and Mezić (2013), and we further assume the gradients of $s_{i}$ are Lipschitz continuous on $\mathcal{B}$, which is required for the perturbative analysis. Note that a non-resonant analyticity of ${\bm{F}}$, which holds generically in practical situations, is sufficient for the Lipschitz continuity, because this assures that $s_{i}$ is analytic.

Let us introduce amplitudes of the system state ${\bm{X}}$ by $r_{i}({\bm{X}})\equiv\mathrm{Re}(s_{i}({\bm{X}}))\ (i=2,\cdots,N)$, where $\mathrm{Re}(z)$ is the real part of a complex number $z$. Because

[TABLE]

each $r_{i}$ obeys

[TABLE]

We can also introduce a phase of ${\bm{X}}$ by $\theta({\bm{X}})\equiv\mathrm{arg}(s_{1}({\bm{X}}))$, where $\mathrm{arg}(z)$ is the argument of $z$, whose range is defined as the interval $[0,2\pi)$. Because $\lambda_{1}=\sqrt{-1}\omega$, $\theta$ obeys

[TABLE]

This definition of the phase coincides with that of the asymptotic phase used in the conventional phase reduction theory Winfree (2001); Kuramoto (2011); Hoppensteadt and Izhikevich (1997); Ermentrout and Terman (2010); Nakao (2016); Ashwin, Coombes, and Nicks (2016). Therefore, level sets of $\theta$ provide isochrons. Analogously, isostables are defined as level sets of $|r_{i}|$. Note that the linear form (4,5), which is valid in the entire basin of attraction Lan and Mezić (2013), is not necessarily derived by the perturbative power-series approach based on the Poincaré-Dulac normal form theory and its extensions Guckenheimer and Holmes (1983); Shilnilov et al. (1998); Gaeta (2002); Wiggins (2003); Sanders, Verhulst, and Murdock (2007). Hence we do not assume the non-resonance condition usually required for a complete linearization in the Poincaré-Dulac type scheme. See Sec. 3.2 of Lan and Mezić’s work Lan and Mezić (2013) for an example with resonance that can be linearized by using Koopman eigenfunctions including non-analytic (trans)monomials.

Because the sign of $r_{i}$ is neglected, each isostable is composed of two connected components corresponding to $+r_{i}$ and $-r_{i}$. These connected components of isostables, associated with one of the exponents $\lambda_{i}$, foliate the basin of attraction of the limit cycle, and each leaf of this foliation provides a level set of the amplitude associated with the exponent. From Eq. (4), we can see that initial conditions on the same isostable share the same decay rate toward the limit cycle. These phase and amplitudes defined above evolve independently under linear time invariant dynamics and thus provide simple description of the dynamics around the limit cycle.

Here, we note that the amplitudes can also be defined as $\tilde{r}_{i}({\bm{X}})\equiv|s_{i}({\bm{X}})|$, as in the preceding study Mauroy, Mezić, and Moehlis (2013). However, this definition makes a coordinate transformation ${\bm{X}}\mapsto(\theta,\tilde{r}_{2},\cdots,\tilde{r}_{N})^{{\dagger}}$ (${\dagger}$ denotes transpose) non-invertible, i.e., its inversion can be multi-valued in some region. The phase-amplitude expression may suffer from this ambiguity, particularly when we apply perturbations to the system. Therefore, we adopt the definition $r_{i}({\bm{X}})\equiv\mathrm{Re}(s_{i}({\bm{X}}))$ in this study.

III Reduction framework and a method to calculate the response functions of the phase and amplitudes

Suppose that perturbation $\epsilon{\bm{p}}(t)$, where $\epsilon>0$ characterizes its magnitude, is introduced to the oscillator (1) as

[TABLE]

We denote a coordinate transformation ${\bm{X}}\mapsto\Theta$ by ${\bm{X}}={\bm{h}}(\Theta)$, where $\Theta=(\theta,r_{2},\cdots,r_{N})^{\dagger}$. In this phase-amplitudes coordinate, the perturbed system (6) takes the following form:

[TABLE]

where $\nabla$ represents gradient and $\cdot$ is a dot product.

Consider a solution $\chi^{*}:{\bm{X}}^{*}(t)$ of the unperturbed system (1) with an initial condition ${\bm{X}}^{*}(0)$ taken arbitrarily in the basin of attraction $\mathcal{B}$, and let $\chi_{p}^{*}:{\bm{X}}_{p}^{*}(t)$ be a solution of the perturbed system (6) with the same initial condition ${\bm{X}}_{p}^{*}(0)={\bm{X}}^{*}(0)$ as the unperturbed system. As is known in a regular perturbation theory Sanders, Verhulst, and Murdock (2007); Coddington and Levinson (1955); Hoppensteadt (2000); Chicone (2006); Teschl (2012), we can show by the Grönwall-Bellman inequality that the magnitude of the error $||{\bm{X}}_{p}^{*}(t)-{\bm{X}}^{*}(t)||$, where $||\cdot||$ denotes the Euclidean norm, is bounded by $b\epsilon(e^{at}-1)/a$, where $a$ and $b$ are positive constants. This means that ${\bm{X}}_{p}^{*}(t)$ is in a neighborhood of radius $\epsilon$ of ${\bm{X}}^{*}(t)$ within a finite time interval of length $O(1)$. We here emphasize that this does not imply the breakdown of the continuous dependence of the solutions on $\epsilon$ within a specific, fixed finite time interval (as long as the unperturbed solution exists on an entire half line, which is the case here). In fact, once we fix an arbitrary large finite length interval $[0,T_{\rm f}]$, we can consider ${\bm{X}}_{p}^{*}(t)$ is in a neighborhood of radius $\epsilon$ of ${\bm{X}}^{*}(t)$ on this interval by taking appropriately small $\epsilon$, because $T_{\rm f}$ is independent of $\epsilon$, and this is sufficient for our argument. The fact that the length of this interval is $O(1)$ means that the convergence of ${\bm{X}}_{p}^{*}(t)$ to ${\bm{X}}^{*}(t)$ is non-uniform on an $\epsilon$-dependent interval $[0,\epsilon^{\beta})$ for any $\beta<0$, i.e., the limiting passages $t\to\epsilon^{\beta}$ and $\epsilon\to+0$ cannot be interchanged. This does not affect our analysis in this study, because no asymptotic properties of the perturbed dynamics are discussed. In this interval, we can expand the gradients using the Lipschitz continuity as $\nabla\theta({\bm{h}}(\Theta))=\nabla\theta({\bm{X}}^{*}(t))+O(\epsilon)$ and $\nabla r_{i}({\bm{h}}(\Theta))=\nabla r_{i}({\bm{X}}^{*}(t))+O(\epsilon)$ in Eqs. (7,8). Thus, we can approximate Eqs. (7,8) as

[TABLE]

by neglecting the terms of order $\epsilon^{2}$.

These equations are completely decoupled from each other and we can adopt combinations of these $N$ equations (9,10) as a reduced form of the system dynamics in the close-enough neighborhood of the transient trajectory $\chi_{*}$. In most cases, the first $K$ equations of (9,10) for some $K(\ll N)$ are of interest, because they describe relatively persistent, slowly decaying modes. Hereafter, we discuss a method to obtain the reduced $K$ equations. The phase and amplitude response functions to perturbation, $\nabla\theta({\bm{X}}^{*}(t))$ and $\nabla r_{i}({\bm{X}}^{*}(t))$, are the fundamental quantities for the proposed reduction framework.

First, we evaluate the gradients on the periodic orbit $\chi$. Consider an initial condition slightly deviated from the periodic orbit, ${\bm{h}}_{p}\equiv{\bm{h}}(\Theta_{1})+\delta{\bm{x}}$, where we defined $\Theta_{1}=(\theta,0,\cdots,0)^{{\dagger}}$. Then

[TABLE]

Using the time-$T$ flow, we can also express $U^{T}r_{i}({\bm{h}}_{p})$ as

[TABLE]

Equating the RHSs of Eqs. (11,12), Taylor expanding $r_{i}$ around ${\bm{h}}(\Theta_{1})$, considering that $r_{i}({\bm{h}}(\Theta_{1}))=0$ and that the direction of $\delta{\bm{x}}$ is arbitrary and taking the limit $||\delta{\bm{x}}||\to 0$, we can show that

[TABLE]

Similarly, we obtain

[TABLE]

Thus, the gradient vectors of the phase and amplitudes evaluated on $\chi$ are left eigenvectors of the monodromy matrix, which are called the adjoint covariant Lyapunov vectors Kuptsov and Parlitz (2012); Traversa et al. (2015); Pikovsky and Politi (2015). These vectors can be numerically obtained by the QR-decomposition based methods Kuptsov and Parlitz (2012); Pikovsky and Politi (2015) or by the spectral dichotomy approaches Froyland et al. (2013); Hüls .

Next, we seek the equations for the gradients of the phase and amplitudes on the transient trajectory $\chi^{*}:{\bm{X}}^{*}(t)$. Here, we introduce logarithmic amplitudes $\psi_{i}({\bm{X}})\equiv\mathrm{log}(|r_{i}({\bm{X}})|)\ (i=2,\cdots,N)$ in order to make the following treatment of the gradients of the amplitudes simple and parallel with the standard arguments in the conventional phase reduction theory. For convenience of notation, let $\psi_{1}({\bm{X}})=\theta({\bm{X}}).$ In the following, we evaluate the gradient vectors of $\psi_{i}$, whose directions coincide with those of $\theta$ and $r_{i}$. The gradients $\nabla\theta$ and $\nabla r_{i}$ can be calculated from $\nabla\psi_{i}$ by rescaling, where the following normalization conditions should be satisfied:

[TABLE]

These normalization conditions are equivalent to Eqs. (4,5).

We can derive adjoint equations for the gradients by using the same argument as the conventional derivation of the adjoint equation for the phase response curves, given by Brown et al. Brown, Moehlis, and Holmes (2004). It is well known that an infinitesimal error $\delta{\bm{x}}(0)$ introduced at $t=0$ between two unperturbed solutions ${\bm{X}}^{*}(t)+\delta{\bm{x}}(t)$ and ${\bm{X}}^{*}(t)$ satisfies the variational equation Guckenheimer and Holmes (1983); Shilnilov et al. (1998); Sanders, Verhulst, and Murdock (2007); Coddington and Levinson (1955); Chicone (2006); Teschl (2012) ${\rm d}(\delta{\bm{x}}(t))/{\rm d}t={\rm D}{\bm{F}}({\bm{X}}^{*}(t))\delta{\bm{x}}(t)$. Because each logarithmic amplitude $\psi_{i}$ increases constantly as $\dot{\psi_{i}}({\bm{X}}(t))=\nabla\psi_{i}({\bm{X}}(t))\cdot\dot{\bm{X}}(t)=\lambda_{i}$ in the absence of perturbation, the error in the logarithmic amplitude coordinate $\psi_{i}({\bm{X}}^{*}(t)+\delta{\bm{x}}(t))-\psi_{i}({\bm{X}}^{*}(t))=\nabla\psi_{i}({\bm{X}}^{*}(t))\cdot\delta{\bm{x}}(t)$ should be independent of time, i.e., ${\rm d}(\nabla\psi_{i}({\bm{X}}^{*}(t))\cdot\delta{\bm{x}}(t))/{\rm d}t=0$. This yields

[TABLE]

Here we used the variational equation and the definition of the adjoint matrix. We can take $N$ linearly independent initial errors $\delta{\bm{x}}_{i}(0)=\epsilon^{\prime}{\bm{e}}_{i}$, where $0<\epsilon^{\prime}\ll 1$ and ${\bm{e}}_{i}$ is the $i$th unit vector and define the fundamental solution matrix ${\bf L}(t)$ of the variational equation as ${\bf L}(t)=(\delta{\bm{x}}_{1}(t),\delta{\bm{x}}_{2}(t),\cdots,\delta{\bm{x}}_{N}(t))$. The sign of the determinant of the fundamental solution matrix, called the Wronskian, is time-invariant due to Liouville’s trace formula Adrianova (1995); Colonius and Kliemann (2014); Wiggins (2003); Coddington and Levinson (1955); Chicone (2006); Teschl (2012). Because ${\rm det}({\bf L}(0))=(\epsilon^{\prime})^{N}>0$, we obtain ${\rm det}({\bf L}(t))>0$ for all $t$, and thus the fundamental solution matrix is always invertible. Consider a matrix form of the Eq. (17), $({\rm d}(\nabla\psi_{i}({\bm{X}}^{*}(t))/{\rm d}t){\bf L}(t)=-{\rm D}{\bm{F}}^{{\dagger}}({\bm{X}}^{*}(t))\nabla\psi_{i}({\bm{X}}^{*}(t)){\bf L}(t)$. We can eliminate ${\bf L}(t)$ by multiplying its inverse from the right side on both sides of this equation. Therefore,

[TABLE]

should hold. Note that this equation should be solved with an appropriate end condition. Here, we can approximately take the end condition of Eq. (18) as $\nabla\psi_{i}({\bm{X}}^{*}(\tau))\parallel\nabla r_{i}({\bm{h}}(\Theta_{1}))|_{\theta=\theta_{*}}$ for some $t=\tau$ and $\theta=\theta_{*}$, because the gradient field $\nabla r_{i}({\bm{X}})$ is continuous and the transient trajectory eventually converges to the limit cycle. The adjoint tangent propagator $\mathcal{G}(t_{1},t_{2})\equiv{\bf N}(t_{2}){\bf N}^{-1}(t_{1})$, where ${\bf N}(t)$ is a fundamental solution matrix of the linear system given by Eq. (18), maps $\nabla\psi_{i}({\bm{X}}^{*}(t_{1}))$ to $\nabla\psi_{i}({\bm{X}}^{*}(t_{2}))$. Thus, $\nabla\theta({\bm{X}}^{*}(t_{2}))\parallel\mathcal{G}(t_{1},t_{2})\nabla\theta({\bm{X}}^{*}(t_{1}))$ and $\nabla r_{i}({\bm{X}}^{*}(t_{2}))\parallel\mathcal{G}(t_{1},t_{2})\nabla r_{i}({\bm{X}}^{*}(t_{1}))$ hold. Therefore, the gradient vectors of the phase and amplitudes are covariant with respect to the action of the propagator $\mathcal{G}$ and they can be interpreted as an extension of the adjoint covariant Lyapunov vectors to transient regimes (note that the adjoint covariant Lyapunov vectors evaluated on the limit cycle, given by Eqs. (13,14), are covariant w.r.t. the action of the adjoint of the monodromy matrix, which is the one period (time-$T$) propagator).

In the numerical estimation of $\nabla\theta$ (or $\nabla\psi_{1}$), a standard method is to integrate the adjoint equation backward in time, while renormalizing $\nabla\theta$ occasionally so that the normalization condition (16) is satisfied Ermentrout and Terman (2010). This is because $\nabla\theta$ corresponds to the neutrally stable component ($\mbox{Re}(\lambda_{1})=0$) while other components have negative growth rates ($\lambda_{2,...,N}<0$). However, in the present case, naive backward integration does not provide correct results for the amplitudes, $\psi_{2,...,N}$, because vector components caused by numerical errors in the relatively (backward-in-time) unstable covariant subspaces accumulate. Therefore, we have to develop a method to subtract them off. Note that the standard QR-decomposition based methods Kuptsov and Parlitz (2012); Pikovsky and Politi (2015) to obtain the covariant subspace require the ergodicity of the underlying dynamical process, hence they cannot be directly applied to the process far from attractors, and that the spectral dichotomy techniques Froyland et al. (2013); Hüls to evaluate them may not work well near the left boundary of the time evolution (See Secs. 2.6 and Sec. 2.7 of Hüls’s work Hüls ).

To develop a numerical method, we introduce dual vectors ${\bm{\gamma}}_{i}$ of $\nabla\psi_{i}$ that are bi-orthogonal to $\nabla\psi_{j}$ as

[TABLE]

where $\delta_{ij}$ is the Kronecker delta. By using ${\bm{\gamma}}_{i}({\bm{X}}^{*}(t))$, we can subtract the vector component in the covariant subspace $\nabla\psi_{i}({\bm{X}}^{*}(t))$ from the solution ${\bm{z}}(t)$ of Eq. (18), which is given by projecting ${\bm{z}}(t)$ onto this subspace as

[TABLE]

Differentiating Eq. (19) by $t$, we obtain $(\dot{{\bm{\gamma}}}_{i}({\bm{X}}^{*}(t))-D{\bm{F}}({\bm{X}}^{*}(t)){\bm{\gamma}}_{i}({\bm{X}}^{*}(t)))\cdot\nabla\psi_{j}({\bm{X}}^{*}(t))=0$. The sign of the Wronskian of Eq. (18) is time-invariant due to Liouville’s trace formula. By using this fact and linear independence of the left eigenvectors of the monodromy matrix, we can show linear independence of $\{\nabla\psi_{i}({\bm{X})}\}_{i=1}^{N}$ for every point ${\bm{X}}$ in the whole basin of attraction $\mathcal{B}$. Thus, we obtain

[TABLE]

The vectors ${\bm{\gamma}}_{i}$ are covariant w.r.t. the action of the propagator ${\mathcal{F}}(=({\mathcal{G}}^{{\dagger}})^{-1})$ of the linear system (21), hence they can be seen as covariant Lyapunov vectors extended to transient regimes. The relative stability relation of covariant subspace of Eq. (21) forward-in-time coincides with that of Eq. (18) backward-in-time. In order to subtract unstable components using the projection (20), the system (21) should be solved forward-in-time with an approximate initial condition ${\bm{\gamma}}_{i}({\bm{X}}^{*}(0))$. The vectors $\{\nabla\psi_{i}({\bm{X}^{*}(0))}\}_{i=1}^{N}$ can be approximated by direct numerical simulation of the dynamics, using the Fourier averages and the generalized Laplace averages Mauroy and Mezić (2012); Mezić (2013) (See Appendix A for details). Then, ${\bm{\gamma}}_{i}({\bm{X}}^{*}(0))$ can be obtained by using the bi-orthogonality relation (19).

Now, we introduce a bi-orthogonalization method to obtain the response functions of the phase and amplitudes up to the $K$th unstable mode. The procedure is as follows: (a) evaluate the adjoint Lyapunov vectors on the limit cycle $\chi$ and the characteristic exponents, (b) calculate $\{{\bm{\gamma}}_{i}({\bm{X}}^{*}(0))\}_{i=1}^{K}$ from $\{\nabla\psi_{i}({\bm{X}^{*}(0))}\}_{i=1}^{N}$ obtained by direct numerical simulation using the bi-orthogonality relation (19), (c) obtain $\nabla\psi_{1}({\bm{X}}^{*}(t))$ by backward integration of Eq. (18), (d) obtain ${\bm{\gamma}}_{1}({\bm{X}}^{*}(t))$ by forward integration of Eq. (21), (e) obtain $\nabla\psi_{2}({\bm{X}}^{*}(t))$ by backward integration of Eq. (18) while subtracting relatively unstable mode $\nabla\psi_{1}({\bm{X}}^{*}(t))$ by the projection (20), (f) obtain ${\bm{\gamma}}_{2}({\bm{X}}^{*}(t))$ by the forward integration of Eq. (21) while subtracting relatively unstable mode ${\bm{\gamma}}_{1}({\bm{X}}^{*}(t))$ by the projection

[TABLE]

where ${\bm{y}}(t)$ is a solution of Eq. (21), (g) perform (e) and (f) consecutively to obtain $\{\nabla\psi_{i}({\bm{X}^{*}(t))}\}_{i=3}^{K}$ and $\{{\bm{\gamma}}_{i}({\bm{X}}^{*}(t))\}_{i=3}^{K}$ (note that all relatively unstable modes should be subtracted during integration), (h) obtain $\nabla\theta$ and $\nabla r_{i}\ (i=2,\cdots,K)$ using the normalization conditions (15,16), where $r_{i}({\bm{X}}^{*}(t))$ on the transient orbit $\chi^{*}$ is evaluated using Eq. (4) with the initial condition $r_{i}({\bm{X}}^{*}(0))$, which is calculated in (b) by the direct numerical simulation.

This method has a significant computational advantages in evaluating the response functions. To calculate response functions $\{\nabla\psi_{i}\}_{i=1}^{K}$ at $m$ points on the transient orbit $\chi^{*}$, it is necessary to repeat long-time evolution $mK(N+1)$ times if we evaluate them directly by the direct numerical simulation. In contrast, we need only $K(N+1)+2K$ times long-time evolution in the proposed bi-orthogonalization method.

IV Examples

As an example, we analyze the Goodwin model, a minimal chemical kinetic model of an oscillatory genetic circuit Goodwin (1965); Gonze and Abou-Jaoudé (2013). The Goodwin model has a three-dimensional state ${\bm{X}}=(x,y,z)^{\dagger}\in\mathbb{R}^{3}$. The state variables $x,y$, and $z$ can be interpreted as concentrations of a given clock mRNA, the corresponding protein, and a transcriptional inhibitor, respectively. We use a simple dimensionless form of the Goodwin model Woller, Gonze, and Erneux (2014),

[TABLE]

The parameters are set as $\alpha=1.8$ and $n=20$. Figure 1(a) shows the stable periodic solution of the model. The period and Lyapunov exponents are estimated as $T=3.63$, $\lambda_{2}=-0.0766$, and $\lambda_{3}=-2.92$. We consider a transient solution ${\bm{X}}^{*}(t)$ with an initial condition ${\bm{X}}^{*}(0)=(1.30,0.900,0.800)^{{\dagger}}$. Figure 1(a) shows the transient solution. We set the end time $\tau=63.0$ for the backward integration in the following calculation.

In Fig. 1(b), the phase response function $\nabla\theta({\bm{X}}^{*}(t))$ obtained by the backward integration of the adjoint equation (18) is compared with the result of the direct numerical simulations. The results agree well, hence, along this transient solution ${\bm{X}}^{*}(t)$, $\nabla\theta({\bm{X}}^{*}(t))$ can always be considered as the most unstable covariant subspace.

Figure 1(c) shows the amplitude response functions $\nabla r_{2}({\bm{X}}^{*}(t))$, which is obtained by the proposed bi-orthogonalization method, by naive backward integration method, and by direct numerical simulations. All results are normalized using the condition (15). Note here that, in the close-enough neighborhood of the limit cycle orbit $\chi$, the vectors $\nabla r_{2}({\bm{X}}^{*}(t))$ and ${\bm{F}}({\bm{X}}^{*}(t))$ are nearly normal. Hence, the normalization procedure using (15) is very sensitive to tiny change in their directions. Therefore, not only the normalization condition (15) but the duality relation (19) must be carefully imposed on the results of the direct numerical simulation in order to make a reasonable comparison with those of the other methods. The results obtained by the naive backward integration considerably deviates from those obtained by direct numerical simulations, while those obtained by the proposed bi-orthogonalization method are in good agreement.

Next, we illustrate the utility of the reduced amplitude equation (10) by estimating the optimal injection timing of weak external input to suppress the most persistent component $r_{2}$ of the amplitudes. We apply a transient control input $\epsilon{\bm{p}}(t)$ of a fixed waveform ${\bm{w}}$ and a fixed duration $\tau_{*}$, i.e., ${\bm{p}}(t)={\bm{w}}(t-s)$ where ${\bm{w}}(\cdot)$ is nonzero only on $[0,\tau_{*}]$ and the time $s$ determines the injection timing of the input. In the spirit of Mauroy’s preceding study Mauroy (2014), we introduce a finite-horizon optimal control problem of minimizing the amplitude $|r_{2}|$ at a given time $T_{e}$. This control problem can be formulated as follows: find the injection timing $s_{*}$ such that

[TABLE]

where $\mathcal{I}_{\sigma}\equiv[0,T_{e}-\tau_{*}]$ and ${\bm{X}}^{*}_{p}(t)$ is the solution of Eq. (6). When the magnitude of the input $\epsilon$ is sufficiently small, the evolution of the amplitude $r_{2}$ is approximated by the reduced equation (10). Then, using an analytical solution of the linear one-dimensional non-homogeneous differential equation (10) of $r_{2}$, the optimal control problem (23) can be approximated to the problem of finding $s_{*}$ such that

[TABLE]

is minimized.

Figure 2 shows the effect of the control input on the amplitude $r_{2}({\bm{X}}_{p}^{*}(T_{e}))$ at time $T_{e}=5$. The control input is assumed as ${\bm{w}}(t)=(0,0,-1)^{{\dagger}}$ and $\tau_{*}=0.25$. The results obtained by the analytical solution of the reduced amplitude equation (10) is compared with the result of direct numerical simulations, showing good agreement for sufficiently weak input ($\epsilon=0.01,0.1$). This verifies the validity of the approximate amplitude equation in the present situation. Thus, the optimal injection timing of sufficiently weak input can be theoretically predicted using the formula (24), because it is essentially equivalent to solving the approximate amplitude equation (10) directly. In this case, the initial value of the amplitude is negative, i.e., $r_{2}({\bm{X}}^{*}(0))<0$. Hence, the optimal injection timing $s_{*}$ of the sufficiently weak input can be estimated by finding the maximum of the waveform in Fig. 2, which gives $s_{*}=2.08$ in this particular case. Finally, we note that when the magnitude becomes large $(\epsilon=1.0)$, the approximation (10) fails and then the results considerably deviate from each other.

V Conclusion

We formulated a phase-amplitude reduction framework for stable limit-cycling systems, which can be applied to transient dynamical regimes far from attractors in high-dimensional systems. We also developed a bi-orthogonalization method for numerical estimation of the response function of the phase and amplitudes, which provides accurate phase-amplitude response functions. As an application, we illustrated that the response functions accurately predicts the optimal injection timing of external input which efficiently suppress deviations from attractors. The proposed theory would be useful in analyzing and controlling response properties of high-dimensional rhythmic systems.

Acknowledgements.

We thank Yoshiyuki Yamaguchi for useful comments on this work. We are also grateful to anonymous reviewers for helpful suggestions and comments. S. S. acknowledges financial support from Japan Society for the Promotion of Science (JSPS) KAKENHI Grant No. 15J12045. W. K. acknowledges financial support from JSPS KAKENHI Grant No. 16K16125. H. N. acknowledges financial support from JSPS KAKENHI Grants No. 16H01538 and No. 16K13847.

Appendix A The Fourier averages and the generalized Laplace averages

In this section, we introduce methods to obtain the phase and amplitudes by direct numerical simulation of the dynamics.

The phase variable $\theta({\bm{X}})$ is evaluated as $\theta({\bm{X}})=\mathrm{arg}(f_{\lambda_{1}}^{*}({\bm{X}}))$, where the Fourier average Mauroy and Mezić (2012) $f_{\lambda_{1}}^{*}({\bm{X}})$ of an observable $f$ is given by

[TABLE]

The amplitude variable $r_{i}({\bm{X}})$ is obtained by $r_{i}({\bm{X}})=\mathrm{Re}(f_{\lambda_{i}}^{*}({\bm{X}}))$, where the generalized Laplace average Mezić (2013) $f_{\lambda_{i}}^{*}({\bm{X}})$ of $f$ is given by

[TABLE]

where $\bar{f}$ is an averaged observable along the periodic orbit $\chi$: $\bar{f}=(1/T)\int_{0}^{T}{f\circ\phi(t,{\bm{X}}_{0}(t_{*}))\mathrm{d}t}$.

We can simplify the generalized Laplace averages using convenient observables $g_{i}\ (i=2,\cdots,N)$ defined as

[TABLE]

where $\theta_{*}=\theta({\bm{X}})$. Here, the adjoint covariant Lyapunov vectors $\nabla r_{i}({\bm{X}}_{0}(\theta_{*}))$ are normalized so that they are dual to the unitized covariant Lyapunov vectors ${\bm{\gamma}}_{i}({\bm{X}}_{0}(\theta_{*}))$. Each of these observables evolves with its corresponding characteristic exponent asymptotically, because, in the close-enough neighborhood of the periodic orbit $\chi$, $g_{i}$ coincides with the $i$th amplitude variable $r_{i}$. Hence, we can show that $\bar{g}_{i}=0$ and $(g_{i})_{\lambda_{k}}^{*}({\bm{X}})=0\ (k=1,\cdots,i-1)$ for any ${\bm{X}}$ in the basin of attraction $\mathcal{B}$. Thus, we can replace the generalized Laplace average with the Laplace average:

[TABLE]