††thanks: Corresponding author. E-mail<EMAIL_ADDRESS>
# Estimating asymptotic phase and amplitude functions of limit-cycle
oscillators from time series data

Norihisa Namura  Department of Systems and Control Engineering, Tokyo
Institute of Technology, Tokyo 152-8552, Japan  Shohei Takata  Department of
Systems and Control Engineering, Tokyo Institute of Technology, Tokyo
152-8552, Japan  Katsunori Yamaguchi  Department of Systems and Control
Engineering, Tokyo Institute of Technology, Tokyo 152-8552, Japan  Ryota
Kobayashi  Graduate School of Frontier Sciences, The University of Tokyo,
Chiba 277-8561, Japan,  Mathematics and Informatics Center, The University of
Tokyo, Tokyo 113-8656, Japan,  JST, PRESTO, Saitama 332-0012, Japan  Hiroya
Nakao  Department of Systems and Control Engineering, Tokyo Institute of
Technology, Tokyo 152-8552, Japan

###### Abstract

We propose a method for estimating the asymptotic phase and amplitude
functions of limit-cycle oscillators using observed time series data without
prior knowledge of their dynamical equations. The estimation is performed by
polynomial regression and can be solved as a convex optimization problem. The
validity of the proposed method is numerically illustrated by using two-
dimensional limit-cycle oscillators as examples. As an application, we
demonstrate data-driven fast entrainment with amplitude suppression using the
optimal periodic input derived from the estimated phase and amplitude
functions.

##  I INTRODUCTION

There are various nonlinear rhythmic phenomena in the real world, including
brain waves Stankovski  _et al._ (2015, 2017a), animal gaits Borgius  _et al._
(2014); Collins and Stewart (1993); Kobayashi  _et al._ (2016); Funato  _et
al._ (2016), heartbeats and respiration Kralemann  _et al._ (2013), passive
walking Hackel (2007); Garcia  _et al._ (1998), many of which can be modeled
mathematically as limit-cycle oscillators Strogatz (2015). The phase reduction
method Kuramoto (1984); Hoppensteadt and Izhikevich (1997); Winfree (2001);
Ermentrout and Terman (2010); Nakao (2016); Monga  _et al._ (2019); Kuramoto
and Nakao (2019); Ermentrout  _et al._ (2019) is useful for analyzing
synchronization properties of limit-cycle oscillators subjected to weak
perturbations, which represents the state of a multidimensional nonlinear
oscillator using only the phase variable introduced along its limit cycle and
describes the dynamics approximately by a one-dimensional phase equation.
Recently, the phase reduction method has been extended to the phase-amplitude
reduction method Wedgwood  _et al._ (2013); Wilson and Moehlis (2016); Mauroy
and Mezić (2016); Shirasaka  _et al._ (2017); Mauroy and Mezić (2018); Monga
_et al._ (2019); Kotani  _et al._ (2020); Shirasaka  _et al._ (2020); Nakao
(2021), which incorporates the amplitude variable characterizing the distance
of the system state from the limit cycle. The resulting phase-amplitude
equation can be used, for example, in deriving the optimal periodic force for
stable entrainment that suppresses amplitude deviations Takata  _et al._
(2021); Shirasaka  _et al._ (2020); Monga and Moehlis (2019).

The phase reduction method is based on the notion of the asymptotic phase
Winfree (2001), but it is generally not possible to analytically obtain the
phase function that gives the asymptotic phase of the system state even if the
explicit mathematical model of the oscillator is available. Similarly, it is
not possible to analytically determine the amplitude functions of the
oscillator in general. Moreover, if the mathematical model of the target
oscillator is unknown, the phase and amplitude functions should be determined
from observed data.

In this study, we propose a simple method for estimating the phase and
amplitude functions from the time series data observed from limit-cycle
oscillators without relying on their mathematical models. Rather than
measuring the values of the phase and amplitude by evolving the system state
until it converges to the limit cycle, we estimate them by polynomial
regression from the differenced time series of the system state started from
various initial conditions. Our method gives a convex optimization problem,
which is computationally inexpensive and can be solved globally. We show that
the phase and amplitude functions are estimated reasonably well around the
limit cycle, including moderately nonlinear regimes, by the proposed method
using known models of limit-cycle oscillators as examples.

For estimating the phase and amplitude functions from observed data, Extended
Dynamic Mode Decomposition (EDMD) Schmid (2010); Kutz  _et al._ (2016);
Williams  _et al._ (2015) and other system-identification methods Lusch  _et
al._ (2018); Wilson (2020, 2021) can also be employed. In particular, EDMD can
estimate the Koopman eigenvalues and eigenfunctions of the system from
observed data Williams  _et al._ (2015); Klus  _et al._ (2020); Li  _et al._
(2017); Takeishi  _et al._ (2017); Alford-Lago  _et al._ (2021); Terao  _et
al._ (2021); Gulina and Mauroy (2021); Kevrekidis  _et al._ (2016); Kutz  _et
al._ (2016); Héas  _et al._ (2020). The natural frequency and Floquet
exponents of the oscillator can then be evaluated from the eigenvalues, and
the phase and amplitude functions can be obtained from the associated
eigenfunctions Mauroy and Mezić (2016, 2018); Mauroy  _et al._ (2020, 2013);
Shirasaka  _et al._ (2017). In contrast to EDMD, our proposed method estimates
the natural frequency and the dominant Floquet exponent separately from the
observed data by using a conventional method for Lyapunov exponents and then
uses them to estimate the phase and amplitude functions by polynomial
regression. Our method thus gives an alternative to EDMD and can more robustly
estimate the phase and amplitude functions depending on the data.

As closely related but different problems, methods for estimating the phase
response properties of limit-cycle oscillators Galán  _et al._ (2005); Ota
_et al._ (2009); Aonishi and Ota (2006); Imai  _et al._ (2017); Cestnik and
Rosenblum (2018); Funato  _et al._ (2016); Nakae  _et al._ (2010); Ota  _et
al._ (2011) and for estimating phase coupling functions between interacting
oscillators from observed data have been proposed and applied to experimental
data Tokuda  _et al._ (2007); Kralemann  _et al._ (2008); Stankovski  _et al._
(2015); Duggento  _et al._ (2012); Stankovski  _et al._ (2017b, a); Ota and
Aoyagi (2014); Onojima  _et al._ (2018); Suzuki  _et al._ (2018); Arai  _et
al._ (2021). In these studies, the phase values of the oscillators are
extracted from the time series, e.g., by extracting the spike timing of neural
oscillators or by using the Hilbert transform, and then the phase response
property of the oscillator or the phase coupling functions between the
oscillators are estimated under the assumption that the oscillators are
described by the phase model. The phase of the oscillator is typically
estimated in the vicinity of the limit cycle and deviation from the limit
cycle is not considered; the phase and amplitude functions are generally not
introduced explicitly.

This paper is organized as follows. We first outline the notions of asymptotic
phase and amplitude in Sec. II and then describe the method for estimating the
phase and amplitude functions by polynomial regression in Sec. III. The
validity of the proposed method is illustrated by numerical simulations using
two types of oscillators and how to evaluate the estimated results is
described in Sec. IV. Section V demonstrates data-driven fast entrainment of
the oscillator with amplitude suppression using the estimated phase and
amplitude functions, and Sec. VI gives concluding remarks.

##  II PHASE AND AMPLITUDE of limit-cycle oscillators

###  II.1 Asymptotic phase and amplitude functions

We first introduce the asymptotic phase and amplitude functions Wedgwood  _et
al._ (2013); Wilson and Moehlis (2016); Shirasaka  _et al._ (2017); Nakao
(2021). We consider a limit-cycle oscillator described by

$\displaystyle\frac{d}{dt}\bm{X}(t)=\bm{F}(\bm{X}(t)),$ (1)

where $\bm{X}(t)\in\mathbb{R}^{N}$ is the system state at time $t$. We assume
that the system has an exponentially stable limit cycle trajectory
$\bm{X}_{0}(t)$ with a natural period $T$ and frequency $\omega=2\pi/T$, which
is a $T$-periodic function of $t$ satisfying $\bm{X}_{0}(t+T)=\bm{X}_{0}(t)$.

First, we assign a phase $\theta\in[0,2\pi)$ for each state on the limit
cycle, where $0$ and $2\pi$ are considered identical, by choosing a state
$\bm{X}_{0}(0)$ on the limit cycle as the phase origin, $\theta=0$, and
defining the phase of the state $\bm{X}_{0}(t)$ at $t>0$ as $\theta=\omega
t\;\;(\mathrm{mod}\;\;2\pi)$. We denote by $\bm{\mathcal{X}}_{0}(\theta)$ the
state on the limit cycle with phase $\theta$. Next, we extend the definition
of the phase to the basin of the limit cycle. The phase of a system state
$\bm{X}_{A}$ in the basin is defined as $\theta$ if
$\lim_{t\to\infty}\left\|\Psi_{t,t_{0}}\bm{X}_{A}-\Psi_{t,t_{0}}\bm{\mathcal{X}}_{0}(\theta)\right\|=0$
holds, where $\Psi_{t,t_{0}}$ is the flow of Eq. (1) and $\|\cdot\|$ is the
$L_{2}$ norm. That is, if the state $\Psi_{t,t_{0}}\bm{X}_{A}$ started from
$\bm{X}_{A}$ at $t_{0}$ converges to the state
$\Psi_{t,t_{0}}\bm{\mathcal{X}}_{0}(\theta)$ started from
$\bm{\mathcal{X}}_{0}(\theta)$ at $t_{0}$ as $t\to\infty$, we consider that
the phase state $\bm{X}_{A}$ has the same phase $\theta$ as
$\bm{\mathcal{X}}_{0}(\theta)$. This defines a phase function $\Theta$ for all
states $\bm{X}$ in the basin, which assigns a phase value
$\theta(t)=\Theta(\bm{X}(t))$ to the state $\bm{X}(t)$ and satisfies

$\displaystyle\frac{d}{dt}\theta(t)=\frac{d}{dt}\Theta(\bm{X}(t))=\omega.$ (2)

The phase defined in this way is called the asymptotic phase, and the level
sets of the phase function is called “isochrons” Kuramoto (1984); Hoppensteadt
and Izhikevich (1997); Winfree (2001); Ermentrout and Terman (2010); Nakao
(2016). The asymptotic phase can also be interpreted as the argument of the
eigenfunction of the Koopman generator $\bm{F}(\bm{X})\cdot\nabla$ of Eq. (1)
associated with the eigenvalue $i\omega$ Mauroy and Mezić (2016, 2018); Mauroy
_et al._ (2020, 2013); Shirasaka  _et al._ (2017).

Next, we introduce the (dominant) amplitude function. In a similar way to the
phase $\theta(t)$, we assign a scalar amplitude $r(t)=R(\bm{X}(t))$ for all
states $\bm{X}$ the basin, where the amplitude function $R$ satisfies

$\displaystyle\frac{d}{dt}r(t)=\frac{d}{dt}R(\bm{X}(t))=\lambda
R(\bm{X}(t))=\lambda r(t).$ (3)

Here, the coefficient $\lambda<0$ is the Floquet exponent of the limit cycle
with the largest non-zero real part, which is assumed to be real and simple.
As the state $\bm{X}(t)$ converges to the limit cycle, the amplitude $r(t)$
decays exponentially and satisfies $R(\bm{\mathcal{X}}_{0}(\theta))=0$ when
the state is on the limit cycle. The amplitude function can also be considered
an eigenfunction of the Koopman generator $\bm{F}(\bm{X})\cdot\nabla$ of Eq.
(1) associated with the eigenvalue $\lambda$. The level sets of the amplitude
function is called “isostables” in the framework of the Koopman operator
analysis of limit-cycling systems Mauroy and Mezić (2016, 2018); Mauroy  _et
al._ (2020, 2013); Shirasaka  _et al._ (2017).

By using the phase and amplitude functions, we can reduce the dimensionality
of the oscillator and globally linearize the dynamics in the basin of the
limit cycle, which yields a simple description of the oscillator useful for
the analysis. Moreover, near the limit cycle, we can write down phase-
amplitude equations describing weakly perturbed oscillatory dynamics. The
simplicity of the phase equation has played a prominent role in the
theoretical analysis of synchronization phenomena in various types of limit-
cycling systems Kuramoto (1984); Hoppensteadt and Izhikevich (1997); Winfree
(2001); Ermentrout and Terman (2010).

###  II.2 Response and sensitivity functions

We next introduce the response and sensitivity functions of the phase and
amplitude, which are used for the validation of our estimation method.

The phase response function (PRF, also known as the phase response or
resetting curve, PRC Kuramoto (1984); Winfree (2001); Brown  _et al._ (2004);
Ermentrout and Terman (2010)) gives the phase difference of the oscillator
caused by an impulse perturbation applied to the oscillator at phase $\theta$
and is defined as

$\displaystyle g(\theta;\bm{\varsigma})$
$\displaystyle=\Theta(\bm{\mathcal{X}}_{0}(\theta)+\bm{\varsigma})-\Theta(\bm{\mathcal{X}}_{0}(\theta))=\Theta(\bm{\mathcal{X}}_{0}(\theta)+\bm{\varsigma})-\theta,$
(4)

where $\bm{\varsigma}\in\mathbb{R}^{N}$ represents the direction and intensity
of the impulse Winfree (2001); Nakao  _et al._ (2005); Ermentrout and Terman
(2010); Nakao (2016); Arai and Nakao (2008). When $|\bm{\varsigma}|$ is
sufficiently small, we can approximate $g$ by using the Taylor expansion of
$\Theta$ near $\bm{\varsigma}=\bm{0}$,
$\Theta(\bm{\mathcal{X}}_{0}(\theta)~{}+~{}\bm{\varsigma})=\Theta(\bm{\mathcal{X}}_{0}(\theta))+\nabla\Theta(\bm{X})|_{\bm{X}=\bm{\mathcal{X}}_{0}(\theta)}~{}\cdot~{}\bm{\varsigma}+~{}O(|\bm{\varsigma}|^{2})$,
as
$g(\theta;\bm{\varsigma})\simeq\nabla\Theta(\bm{X})|_{\bm{X}=\bm{\mathcal{X}}_{0}(\theta)}\cdot\bm{\varsigma}$.
We denote the gradient of $\Theta(\bm{X})$ evaluated at the state
$\bm{X}=\bm{\mathcal{X}}_{0}(\theta)$ by

$\displaystyle\bm{Z}(\theta)=\nabla\Theta(\bm{X})|_{\bm{X}=\bm{\mathcal{X}}_{0}(\theta)}$
(5)

and call it the phase sensitivity function (PSF, also known as the
infinitesimal phase resetting curve, iPRC). This $\bm{Z}(\theta)$
characterizes the linear response of the oscillator phase to a weak
perturbation given at phase $\theta$.

Similarly, we introduce the amplitude response function (ARF) characterizing
the amplitude difference of the oscillator caused by an impulse
$\bm{\varsigma}$ given at phase $\theta$ as

$\displaystyle h(\theta;\bm{\varsigma})$
$\displaystyle=R(\bm{\mathcal{X}}_{0}(\theta)+\bm{\varsigma})-R(\bm{\mathcal{X}}_{0}(\theta))=R(\bm{\mathcal{X}}_{0}(\theta)+\bm{\varsigma}).$
(6)

When $|\bm{\varsigma}|$ is sufficiently small, we can approximate $h$ by using
the Taylor expansion of $R(\bm{\mathcal{X}}_{0}(\theta)+\bm{\varsigma})$ near
$\bm{\varsigma}=\bm{0}$,
$R(\bm{\mathcal{X}}_{0}(\theta)+\bm{\varsigma})=R(\bm{\mathcal{X}}_{0}(\theta))+\nabla
R(\bm{X})|_{\bm{X}=\bm{\mathcal{X}}_{0}(\theta)}\cdot\bm{\varsigma}+O(|\bm{\varsigma}|^{2})$,
as $h(\theta;\bm{\varsigma})\simeq{\nabla
R(\bm{X})|_{\bm{X}=\bm{\mathcal{X}}_{0}(\theta)}}\cdot\bm{\varsigma}$. Here,
we denote the gradient of the amplitude function $R(\bm{X})$ at the state
$\bm{X}=\bm{\mathcal{X}}_{0}(\theta)$ by

$\displaystyle\bm{I}(\theta)=\nabla
R(\bm{X})|_{\bm{X}=\bm{\mathcal{X}}_{0}(\theta)}$ (7)

and call it the amplitude sensitivity function (ASF, also known as isostable
response function in the Koopman operator analysis of limit-cycling systems).
This $\bm{I}(\theta)$ characterizes the linear response of the oscillator
amplitude to a weak perturbation given at $\theta$.

If the mathematical model of the oscillator is known, the PSF $\bm{Z}(\theta)$
and ASF $\bm{I}(\theta)$ can be determined by numerically calculating the
$2\pi$-periodic solutions of the following adjoint equations Brown  _et al._
(2004); Ermentrout and Terman (2010); Kuramoto and Nakao (2019); Shirasaka
_et al._ (2017); Takata  _et al._ (2021):

$\displaystyle\omega\frac{d}{d\theta}\bm{Z}(\theta)$
$\displaystyle=-\bm{J}(\theta)^{\top}\bm{Z}(\theta),$ (8)
$\displaystyle\omega\frac{d}{d\theta}\bm{I}(\theta)$
$\displaystyle=-\left(\bm{J}(\theta)^{\top}-\lambda\right)\bm{I}(\theta),$ (9)

where $\bm{J}(\theta)$ denotes the Jacobian matrix of the vector field
$\bm{F}$ at $\bm{X}=\bm{\mathcal{X}}_{0}(\theta)$. The PSF $\bm{Z}(\theta)$
should satisfy
$\bm{Z}(\theta)\cdot\bm{F}\left(\bm{\mathcal{X}}_{0}(\theta)\right)=\omega$ as
a normalization condition. The normalization condition for the ASF can be
chosen arbitrary and will be specified later.

###  II.3 Direct numerical calculation of phase and amplitude functions

It is generally difficult to obtain the phase and amplitude functions
analytically from a mathematical model, but they can be obtained by direct
numerical calculation of the mathematical model as follows.

For the phase function Nakao (2016), we first choose a point on the limit
cycle as the origin of the phase, then assign the phase values to the states
on the limit cycle. For the states in the basin of the limit cycle, we let
them evolve for an integer multiple of the period $T$ until they converge to
the limit cycle and then determine their phase values.

Similarly, for the amplitude function, we let the state $\bm{X}$ in the basin
evolve until it becomes sufficiently close to (but not completely on) the
limit cycle, and find the point $\bm{\mathcal{X}}(\theta)$ on the limit cycle
with the same phase $\theta$ as $\bm{X}$. The amplitude of the point
sufficiently close to the limit cycle can be obtained by taking the inner
product of the vector $\Delta\bm{X}=\bm{X}-\bm{\mathcal{X}}(\theta)$ between
the two points and the amplitude sensitivity function $\bm{I}(\theta)$ as
$R(\bm{X})\simeq\bm{I}(\theta)\cdot\Delta\bm{X}$, which follows from the
Taylor expansion of the amplitude function.

Note that these methods are exhaustive and computationally required. We will
use the phase and amplitude functions directly calculated from the
mathematical model by the above method as ‘true’ functions to characterize the
accuracy of the estimated functions from the observed data.

##  III Proposed method of estimation

###  III.1 Estimation of the frequency and Floquet exponent

In our method, the natural frequency $\omega$ and the largest non-zero Floquet
exponent $\lambda$ of the oscillator should be measured before the estimation
of the phase and amplitude functions. The natural frequency $\omega$ can be
estimated from the average period $T$ of the system state to perform
rotations. To estimate the sum of the Floquet exponent from time series data,
we use the method of Wolf et al. Wolf  _et al._ (1985) for the Lyapunov
exponents of dynamical systems. For two-dimensional limit-cycle oscillators,
the Floquet exponents are real and exactly the same as the Lyapunov exponents.
As one Lyapunov exponent $\lambda_{1}$ corresponding to the phase direction is
$0$, it is enough to estimate the sum of the Lyapunov exponents,
$\lambda_{1}+\lambda_{2}$, which can be performed by measuring the evolution
of a small phase-plane area for a long time. The same method can also be used
for the two largest Lyapunov exponents of higher-dimensional oscillators when
$\lambda$ is real.

We assume that a discretized time sequence of the system state (with
observation noise) from an initial condition not on the limit cycle is
obtained. By taking the logarithmic ratio of the area $S(t_{k})$ of the
triangle formed by three nearby states in the time series at time $t_{k}$ to
the area $S(t_{k+1})$ at $t_{k+1}$ and summing over a long period of time, we
can obtain $\lambda$ as the sum of the Lyapunov exponents:

$\displaystyle\lambda=\lambda_{1}+\lambda_{2}=\frac{1}{t_{L}-t_{0}}\sum_{k=0}^{L-1}\log\frac{S(t_{k+1})}{S(t_{k})},$
(10)

where $L$ is the length of the data and the area $S(t_{k})$ is reset to $0.1$
at each $t_{k}$ so that it does not vanish due to numerical errors.

###  III.2 Estimation of the phase function

We first propose a method for estimating the phase function by polynomial
regression from observed data of an unknown limit-cycle oscillator. We do not
measure the absolute phase value of each state, which is difficult with the
time series; rather, we use the phase differences between two consecutive
states in the time series and use Eq. (2) for the estimation.

We assume that we can obtain discretized time sequences of
$\bm{X}(t)\in\mathbb{R}^{N\times 1}$ from various initial conditions (with
observation noise). Each sequence is sampled at equal intervals of $\Delta t$
with $M$ sampling points, where $\Delta t$ is sufficiently smaller than the
natural period $T$ of the oscillator. We construct a standardized row vector
of polynomials up to order $p\geq 1$ from $\bm{X}$,

$\displaystyle\bm{U}(\bm{X},p)=$
$\displaystyle\quad\begin{bmatrix}1&\bar{x}_{1}&\cdots&\bar{x}_{N}&\bar{x}^{2}_{1}&\overline{x_{1}x_{2}}&\cdots&\bar{x}_{N}^{2}&\cdots&\bar{x}_{N}^{p}\end{bmatrix}\in\mathbb{R}^{1\times
P},$ (11)

where $x_{k}$ is the $k$-th component of $\bm{X}$ and all cross terms like
$\overline{x_{1}^{2}x_{2}^{3}}$ up to order $p$ are included. The overline
indicates standardization, i.e., $\bar{z}=:(z-\mu)/\sigma$ is the standard
score of $z$, where $\mu$ and $\sigma$ are the mean and the standard deviation
of $z$, respectively. We consider approximating the sine and cosine of the
phase function $\Theta(\bm{X})$ by polynomials of order $p$ as

$\displaystyle\begin{aligned}
\cos(\Theta(\bm{X}))\simeq\bm{U}(\bm{X},p)\bm{z}_{1},\quad\sin(\Theta(\bm{X}))\simeq\bm{U}(\bm{X},p)\bm{z}_{2},\end{aligned}$
(12)

where $\bm{z}_{1}\in\mathbb{R}^{P\times 1}$ and
$\bm{z}_{2}\in\mathbb{R}^{P\times 1}$ are the coefficient vectors of the
polynomials. Hereafter, $\bm{U}(\bm{X},p)$ is denoted by $\bm{U}(\bm{X})$. Our
aim is to find the best polynomial approximations of $\sin\Theta$ and
$\cos\Theta$ that are consistent with the definition of $\Theta$, Eq. (2),
from the time series of $\bm{X}(t)$.

Taking the time derivatives of $\sin\Theta(\bm{X})$ and $\cos\Theta(\bm{X})$,
we obtain

$\displaystyle\begin{aligned}
\frac{d}{dt}\cos(\Theta(\bm{X}(t)))&=-\omega\sin(\Theta(\bm{X}(t))),\\\
\frac{d}{dt}\sin(\Theta(\bm{X}(t)))&=\omega\cos(\Theta(\bm{X}(t))),\end{aligned}$
(13)

and using Eq. (12), we obtain

$\displaystyle\begin{aligned}
\frac{d}{dt}\bm{U}(\bm{X}(t))\bm{z}_{1}&=\left(\left(\left.\frac{\partial\bm{U}(\bm{X})}{\partial\bm{X}}\right|_{\bm{X}=\bm{X}(t)}\right)^{\top}\frac{d\bm{X}(t)}{dt}\right)\cdot\bm{z}_{1}\\\
&\simeq-\omega\bm{U}(\bm{X}(t))\bm{z}_{2},\\\
\frac{d}{dt}\bm{U}(\bm{X}(t))\bm{z}_{2}&=\left(\left(\left.\frac{\partial\bm{U}(\bm{X})}{\partial\bm{X}}\right|_{\bm{X}=\bm{X}(t)}\right)^{\top}\frac{d\bm{X}(t)}{dt}\right)\cdot\bm{z}_{2}\\\
&\simeq\omega\bm{U}(\bm{X}(t))\bm{z}_{1},\end{aligned}$ (14)

where $\partial{\bm{U}(\bm{X})}/\partial{\bm{X}}\in\mathbb{R}^{N\times P}$ is
a matrix whose $(j,k)$-component is given by $\partial U_{k}/\partial X_{j}$
($j=1,\dots,N,k=1,\dots,P$) and “$\cdot$” indicates the dot product of two
vectors.

We discretize the time as $t=i\Delta t$ with a small time step $\Delta t$ and
represent the time series $\bm{X}(t)$ as $\bm{X}_{i}=\bm{X}(i\Delta t)$. We
also denote numerical derivative of $\bm{X}(t)$ at $t=i\Delta t$ as
$\dot{\bm{X}}_{i}$, which is calculated by a linear regression of several
consecutive data points in the time series since observation noise is
overlapped on the data Knowles and Renka (2014). From Eq. (14), the following
equations should hold approximately for each $i=1,\dots,M$: Introducing
$\bm{z}=\left[\bm{z}_{1}^{\top},\bm{z}_{2}^{\top}\right]^{\top}\in\mathbb{R}^{2P\times
1}$, these equations can be expressed as

$\displaystyle\begin{aligned}
\begin{bmatrix}(\bm{U}^{\prime}_{i})^{\top}\dot{\bm{X}}_{i}&\omega\bm{U}_{i}\\\
-\omega\bm{U}_{i}&(\bm{U}^{\prime}_{i})^{\top}\dot{\bm{X}}_{i}\end{bmatrix}\bm{z}\coloneqq{}&\bm{c}_{i}\bm{z}={}\begin{bmatrix}0\\\
0\end{bmatrix},\end{aligned}$ (15)

where $\bm{U}_{i}=\bm{U}(\bm{X}_{i})$ and
$\bm{U}^{\prime}(\bm{X}_{i})=\partial
U(\bm{X})/\partial\bm{X}|_{\bm{X}=\bm{X}_{i}}$. We try to find the best
coefficient vector $\bm{z}$ that satisfies the above equation as much as
possible for the given time series. Defining a matrix $\bm{C}$ by
$\bm{C}=\left[\bm{c}_{1}^{\top},\bm{c}_{2}^{\top},\cdots,\bm{c}_{M}^{\top}\right]^{\top}\in\mathbb{R}^{2M\times
2P}$, this gives a minimization problem of the following objective function:

$\displaystyle E_{\Theta}=\|\bm{C}\bm{z}\|^{2}$ (16)

where $\|\bm{C}\bm{z}\|^{2}=\sum_{k=1}^{2M}\\{(\bm{C}\bm{z})_{k}\\}^{2}$.

To fix the origin of the estimated phase function, we also require that the
absolute phase of a given state $\tilde{\bm{X}}$ is zero, i.e.,

$\displaystyle\Theta\left(\tilde{\bm{X}}\right)=0.$ (17)

Thus, the polynomial approximations of $\cos\Theta\left(\tilde{\bm{X}}\right)$
and $\sin\Theta\left(\tilde{\bm{X}}\right)$ should satisfy

$\displaystyle\begin{aligned}
\bm{U}\left(\tilde{\bm{X}}\right)\bm{z}_{1}=1,\quad\bm{U}\left(\tilde{\bm{X}}\right)\bm{z}_{2}=0,\end{aligned}$
(18)

and using $\bm{z}=\left[\bm{z}_{1}^{\top},\bm{z}_{2}^{\top}\right]^{\top}$,

$\displaystyle\begin{bmatrix}\bm{U}\left(\tilde{\bm{X}}\right)&\bm{0}\\\
\bm{0}&\bm{U}\left(\tilde{\bm{X}}\right)\end{bmatrix}\bm{z}=\begin{bmatrix}1\\\
0\end{bmatrix}.$ (19)

Considering the above objective function and the constraints, the proposed
estimation method can be formulated as an optimization problem for the
coefficient $\bm{z}$ as follows:

$\displaystyle\begin{aligned} &&\hat{\bm{z}}={}&\mathop{\rm
argmin}\limits_{\bm{z}}\|\bm{C}\bm{z}\|^{2}\\\
&&&\mathrm{s.t.}\begin{bmatrix}\bm{U}\left(\tilde{\bm{X}}\right)&\bm{0}\\\
\bm{0}&\bm{U}\left(\tilde{\bm{X}}\right)\end{bmatrix}\bm{z}=\begin{bmatrix}1\\\
0\end{bmatrix}.\end{aligned}$ (20)

This optimization problem is convex with respect to the coefficient vector
$\bm{z}$.

###  III.3 Estimation of the amplitude function

We next develop a method for estimating the amplitude function by polynomial
regression using the same observed data as in the previous subsection. For
simplicity, we consider the amplitude function of a two-dimensional
oscillator; the method can also be applied straightforwardly to higher-
dimensional oscillators whose Floquet exponent with the largest non-zero part
is real.

Since the Floquet exponents of a two-dimensional oscillator are real, i.e., a
zero exponent associated with the phase (tangential) direction and another
negative exponent associated with the amplitude direction, the amplitude
function is also real. Therefore, we can approximate the amplitude function
$R(\bm{X})$ by a polynomial as

$\displaystyle R(\bm{X})\simeq\bm{U}(\bm{X})\bm{\zeta},$ (21)

where $\bm{\zeta}\in\mathbb{R}^{P\times 1}$ is the coefficient vector.
Plugging Eq. (21) into Eq. (3), we obtain

$\displaystyle\left(\left(\left.\frac{\partial\bm{U}(\bm{X})}{\partial\bm{X}}\right|_{\bm{X}=\bm{X}(t)}\right)^{\top}\frac{d\bm{X}(t)}{dt}\right)\cdot\bm{\zeta}=\lambda\bm{U}(\bm{X})\bm{\zeta}.$
(22)

Discretizing the time and plugging the time series data, we obtain

$\displaystyle\begin{bmatrix}\left(\left.\frac{\partial\bm{U}(\bm{X})}{\partial\bm{X}}\right|_{\bm{X}=\bm{X}_{i}}\right)^{\top}\dot{\bm{X}}_{i}-\lambda\bm{U}(\bm{X}_{i})\end{bmatrix}\bm{\zeta}\coloneqq\bm{d}_{i}\bm{\zeta}=0$
(23)

for $i=1,...,M$, where $\bm{d}_{i}\in\mathbb{R}^{1\times P}$. As in the case
for the phase function, the numerical differentiation of $\bm{X}(t)$ is
evaluated by the polynomial interpolation. Defining a matrix $\bm{D}$ by
$\bm{D}=\left[\bm{d}_{1}^{\top},\bm{d}_{2}^{\top},\cdots,\bm{d}_{M}^{\top}\right]^{\top}\in\mathbb{R}^{M\times
P}$, we obtain the following objective function to minimize:

$\displaystyle E_{R}=\|\bm{D}\bm{\zeta}\|^{2}.$ (24)

On the limit cycle, the value of the amplitude function should take $0$, i.e.,
$R(\bm{\mathcal{X}}_{0}(\theta))=0$. This condition is given as the constraint
for the polynomial approximation of $R$ as $\bm{U}(\bm{X}_{0,i})\bm{\zeta}=0$,
where $\bm{X}_{0,i}$ $(i=1,...,M_{0})$ denotes $M_{0}$ oscillator states on
the limit cycle, which can be chosen at equal time intervals. Defining a
matrix $\bm{U}_{0}$ by
$\bm{U}_{0}=\left[\bm{U}(\bm{X}_{0,1})^{\top},\bm{U}(\bm{X}_{0,2})^{\top},\cdots,\bm{U}(\bm{X}_{0,M_{0}})^{\top}\right]^{\top}$,
we obtain

$\displaystyle\bm{U}_{0}\bm{\zeta}=\bm{0}.$ (25)

As stated previously, the scale of the amplitude can be arbitrarily chosen. To
fix the scale of the amplitude function (this also fixes the scale of the
amplitude sensitivity function $\bm{I}(\theta)$ introduced in Sec. II.2, we
specify the absolute amplitude value of a given state $\tilde{\bm{X}}$ as
$R\left(\tilde{\bm{X}}\right)=r_{0}$ ($r_{0}>0$) as in the case of the phase
function. This point $\tilde{\bm{X}}$ should not be on the limit cycle and
should be chosen near the fixed point of the system in order to capture the
shape of the amplitude function appropriately. This constraint can be
expressed as

$\displaystyle\bm{U}\left(\tilde{\bm{X}}\right)\bm{\zeta}\coloneqq\tilde{\bm{U}}\bm{\zeta}=r_{0}.$
(26)

Moreover, we require that the norm of the coefficient vector $\bm{\zeta}$ does
not become too large so that the estimated amplitude function does not diverge
outside of the limit cycle due to the above two constraints near the fixed
point and on the limit cycle. We thus include an additional objective function
proportional to the norm $\|\bm{\zeta}\|^{2}$.

Summarizing, the proposed estimation method can be formulated as an
optimization problem for the coefficient vector $\bm{\zeta}$ as follows:

$\displaystyle\begin{aligned} &&\hat{\bm{\zeta}}={}&\mathop{\rm
argmin}\limits_{\bm{\zeta}}\left\|\bm{D}\bm{\zeta}\right\|^{2}+\gamma\|\bm{\zeta}\|^{2}\\\
&&&\mathrm{s.t.}\begin{bmatrix}\tilde{\bm{U}}\\\
\bm{U}_{0}\end{bmatrix}\bm{\zeta}=\begin{bmatrix}r_{0}\\\
\bm{0}\end{bmatrix},\end{aligned}$ (27)

where $\gamma$ is a parameter for the penalty for the solution norm. This
optimization problem is convex with respect to the coefficient vector
$\bm{\zeta}$.

The second term of the objective function in Eq. (27) takes the same form as
the regularization term in the ridge regression. The hyper parameter $\gamma$
was determined by using the L-curve Hansen (1992). The L-curve is defined by a
log-log plot of the error $\|D\bm{\zeta}\|^{2}$ and the solution norm
$\|\bm{\zeta}\|^{2}$. We empirically choose the value of $\gamma$ where the
absolute value of the slope falls below $6$.

##  IV Verification of the proposed method

###  IV.1 Data used for estimation

In this section, we verify the validity of the proposed method by numerical
simulations using two-dimensional limit-cycle oscillators, the Stuart-Landau
oscillator Kuramoto (1984); Nakao (2016) and the van der Pol oscillator van
der Pol (1927); van der Pol and van der Mark (1927); van der Pol (1926). The
time series data are produced by numerical integration of the dynamical system
from $n$ initial states that are taken uniformly at random in a certain
region, where each initial state is evolved until it converges to the limit
cycle and $M$ data points are sampled from the trajectory. Here, we took $n$
initial points to increase the number of data points outside the limit cycle,
while we assumed a single time series in our explanation of the methods in
Sec. III. These $n$ time series are concatenated and treated as one time
series vector, but the regression and differentiation are performed only
within the individual time series.

###  IV.2 Evaluation of estimated results

To examine the performance of the proposed method, we consider limit-cycle
oscillators whose mathematical models are known and compare the estimated
results with the true results obtained by direct numerical calculations of the
mathematical models. We evaluate the estimation accuracy by comparing the
PSFs, PRFs, ASFs, and ARFs obtained from the estimated and the true phase and
amplitude functions. In what follows, we use normalized PRFs (nPRFs) and
normalized ARFs (nARFs) defined as the PRFs and ARFs divided by the impulse
intensity, respectively.

The $j$-th component of the PSF $\bm{Z}(\theta)$ and ASF $\bm{I}(\theta)$ are
obtained by numerical derivative of the estimated phase and amplitude
functions as

$\displaystyle\hat{Z}_{j}(\theta_{i})=\frac{\hat{\Theta}(\bm{X}_{0,i}+\epsilon\bm{e}_{j})-\hat{\Theta}(\bm{X}_{0,i}-\epsilon\bm{e}_{j})}{2\epsilon}$
(28)

and

$\displaystyle\hat{I}_{j}(\theta_{i})=\frac{\hat{R}(\bm{X}_{0,i}+\epsilon\bm{e}_{j})-\hat{R}(\bm{X}_{0,i}-\epsilon\bm{e}_{j})}{2\epsilon},$
(29)

for sufficiently small $\epsilon$, where $\bm{e}_{j}$ is the unit vector in
the $j$-th direction and $\bm{X}_{0,i}$ is the system state on the limit cycle
with phase $\theta_{i}$. We use $\epsilon=10^{-5}$ in what follows.

The nPRF ${G}_{\varsigma,j}(\theta)\coloneqq
g(\theta,\varsigma\bm{e}_{j})/\varsigma$ and the nARF
$H_{\varsigma,j}(\theta)\coloneqq h(\theta,\varsigma\bm{e}_{j})/\varsigma$ of
the oscillator with respect to the impulse of finite intensity $\varsigma$
given to the $j$-th component of the system state can be calculated from the
estimated phase and amplitude functions as

$\displaystyle\hat{G}_{\varsigma,j}(\theta_{i})\coloneqq\frac{\hat{\Theta}(\bm{X}_{0,i}+\varsigma\bm{e}_{j})-\hat{\Theta}(\bm{X}_{0,i})}{\varsigma},$
(30)

and

$\displaystyle\hat{H}_{\varsigma,j}(\theta_{i})\coloneqq\frac{\hat{R}(\bm{X}_{0,i}+\varsigma\bm{e}_{j})-\hat{R}(\bm{X}_{0,i})}{\varsigma},$
(31)

respectively.

We evaluate the accuracy of the PSFs, ASFs, nPRFs, and nARFs by using the
coefficient of determination as follows:

$\displaystyle
R_{Y}^{2}=1-\frac{\sum_{i}\left(Y_{j}(\theta_{i})-\hat{Y}_{j}(\theta_{i})\right)^{2}}{\sum_{i}\left(Y_{j}(\theta_{i})-\overline{Y}_{j}\right)^{2}},$
(32)

where $Y_{j}$ is either of $Z_{j}$, $I_{j}$, $G_{\varsigma,j}$, or
$H_{\varsigma,j}$, and $\overline{Y}_{j}$ represents the mean of
$\\{Y_{j}(\theta_{i})\\}$ over $i$. The closer the coefficient of
determination is to $1$, the higher the accuracy of the estimation.

Figure 1:  Phase and amplitude functions of the SL oscillator used for the
estimation. (a) Time series data used for the estimation. (b) Estimated non-
zero Floquet exponent. (c) Estimated phase function. (d) True phase function.
(e) Estimated amplitude function. (f) True amplitude function. In (c)-(f), the
red bold circle shows the limit cycle.  Figure 2:  Phase sensitivity functions
(PSFs) of the SL oscillator. Each graph shows the estimated PSF (red solid
curve) and the true PSF (black dotted curve). (a): $x_{1}$ component. (b):
$x_{2}$ component. Estimated PSFs are obtained by numerical differentiation
with $\epsilon=10^{-5}$. True PSFs are calculated by solving the adjoint
equation.  Figure 3:  Normalized phase response functions (nPRFs) of the SL
oscillator. Each graph shows the estimated nPRF (red solid curve) and the true
nPRF (black dotted curve). (a,c): Impulses applied in the $x_{1}$ direction.
(a) $\varsigma=-0.2$, (c) $\varsigma=+0.2$. (b,d): Impulses applied in the
$x_{2}$ direction. (b) $\varsigma=-0.2$, (d) $\varsigma=+0.2$.

###  IV.3 Stuart-Landau oscillator

First, we verify the validity of the proposed method using the Stuart-Landau
(SL) oscillator, for which the analytical solutions of the phase and amplitude
functions are obtained Kato  _et al._ (2021). The SL oscillator is described
as follows:

$\displaystyle\frac{d}{dt}\begin{bmatrix}x_{1}\\\
x_{2}\end{bmatrix}=\begin{bmatrix}x_{1}-\alpha x_{2}-(x_{1}-\beta
x_{2})\left(x_{1}^{2}+x_{2}^{2}\right)\\\ \alpha x_{1}+x_{2}-(\beta
x_{1}+x_{2})\left(x_{1}^{2}+x_{2}^{2}\right)\end{bmatrix},$ (33)

where $x_{1},x_{2}$ are the variables and $\alpha=2,\;\beta=1$ are parameters.
We set the sampling interval as $\Delta t=0.005$, the number of initial points
as $n=1200$, and the number of data points as $M=500$. We estimated the phase
function for the SL oscillator under the observation noise, which is given by
Gaussian noise with mean $0$ and standard deviation $5\times 10^{-3}$ added to
the time series. The data used for the estimation is shown in Fig. 1(a). The
Floquet exponents are estimated with $L=1000$ and $t_{k+1}-t_{k}=0.25$ in Eq.
(10). The natural frequency is estimated as $\omega\simeq 0.9997$ and the non-
zero Floquet exponent is estimated as $\lambda\simeq-2.0457$, while their
theoretical values are $\omega=1$ and $\lambda=-2$, respectively. The
estimated value of the Floquet exponent is plotted with respect to the length
$L$ in Fig. 1(b). The maximum degree of the polynomial is set to $p=18$.

We estimated the phase function ${\Theta}$ by the proposed method (Sec. III.2)
from the time series data (Fig. 1(a)). The estimated phase function (Fig.
1(c)) is compared with the true one (Fig. 1(d)). The phase function is
accurately estimated by the proposed method except for the central and
peripheral regions far from the limit cycle. The discrepancy is due to the
lack of data in those regions (trajectories stay near the limit-cycle
attractor most of the time) and also due to the limited representation power
of the polynomials.

We evaluated the estimation accuracy of the PSFs and PRFs. Figure 2 compares
the PSFs $Z_{1}$ and $Z_{2}$ obtained from the estimated phase function using
Eq. (28) with the true function obtained by solving the adjoint equation (8).
Figure 3 compares the nPRFs $G_{\varsigma,j}$ for finite-intensity impulses
obtained from the estimated phase function using Eq. (30) with the true
functions obtained from the analytical solution. The estimated functions are
well consistent with the true ones and the estimation error is small; the
accuracy of the estimation is $R_{Z}^{2}\geq 0.9859$ for the PSFs and
$R_{G}^{2}\geq 0.9903$ for the nPRFs, respectively (Table 2). These results
indicate that the proposed method accurately estimates the phase function in
the vicinity of the limit-cycle and also a bit farther out from the limit
cycle where moderate nonlinearity comes in.

Next, we estimated the amplitude function $R$ by the proposed method (Sec.
III.3) from the same time series data (Fig. 1(a)). The hyper parameter
$\gamma$ is determined to be $\gamma=1.0\times 10^{6}$ from the L-curve. The
estimated amplitude function (Fig. 1(e)) is compared with the true amplitude
function (Fig. 1(f)). Our method can accurately estimate the amplitude
function near the limit cycle, but the estimate deviates in the central and
peripheral regions. In particular, the estimate cannot reproduce the
divergence of the amplitude function near the fixed point. This discrepancy is
due to the lack of data and the singularity near the fixed point of the
system.

Figure 4 compares the ASFs $I_{1}$ and $I_{2}$ obtained from the estimated
amplitude function using Eq. (29) with the true function obtained by solving
the adjoint equation (9). Figure 5 compares the nARFs
$H_{\varsigma,j}(\theta)$ obtained from the estimated amplitude function using
Eq. (31) with the true results. Despite the discrepancy in the amplitude
function far from the limit cycle, the estimated results agree reasonably well
with the true results for the impulse intensities used here. The accuracy of
estimation is $R_{I}^{2}\geq 0.9998$ for the ASFs and $R_{H}^{2}\geq 0.9725$
for the nARFs (Table 2). These results confirm that the amplitude function is
estimated appropriately around the limit cycle including moderately nonlinear
regimes.

Figure 4:  Amplitude sensitivity functions (ASFs) of the SL oscillator. Each
graph shows the estimated ASF (red solid curve) and the true ASF (black dotted
curve). (a) $x_{1}$ component. (b) $x_{2}$ component. Estimated ASFs are
obtained by numerical differentiation with $\epsilon=10^{-5}$. True ASFs are
calculated by solving the adjoint equation.  Figure 5:  Normalized amplitude
response functions (nARFs) of the SL oscillator. Each graph shows the
estimated nARF (red solid curve) and the true nARF (black dotted curve).
(a,c): Impulses applied in the $x_{1}$ direction. (a) $\varsigma=-0.1$, (c)
$\varsigma=+0.1$. (b,d): Impulses applied in the $x_{2}$ direction. (b)
$\varsigma=-0.1$, (d) $\varsigma=+0.1$.  Table 1: Coefficients of
determination of the nPRFs and PSFs of the SL oscillator.

$Z_{1}$ | $Z_{2}$ | $G_{-0.2,1}$ | $G_{-0.2,2}$ | $G_{0.2,1}$ | $G_{0.2,2}$
---|---|---|---|---|---
$0.9869$ | $0.9859$ | $0.9912$ | $0.9903$ | $0.9927$ | $0.9929$

Table 2: Coefficients of determination of the nARFs and ASFs of the SL
oscillator.

$I_{1}$ | $I_{2}$ | $H_{-0.1,1}$ | $H_{-0.1,2}$ | $H_{0.1,1}$ | $H_{0.1,2}$
---|---|---|---|---|---
$0.9998$ | $0.9999$ | $0.9779$ | $0.9727$ | $0.9725$ | $0.9742$

###  IV.4 van der Pol oscillator

Next, we verify the validity of the proposed method using the van der Pol
(vdP) oscillator. The vdP oscillator van der Pol (1926, 1927); van der Pol and
van der Mark (1927) is described as follows:

$\displaystyle\frac{d}{dt}\begin{bmatrix}x_{1}\\\
x_{2}\end{bmatrix}=\begin{bmatrix}x_{2}\\\
\nu\left(1-x_{1}^{2}\right)x_{2}-x_{1}\end{bmatrix},$ (38)

where $x_{1},x_{2}$ are the variables and the parameter is chosen as $\nu=1$.
We set the sampling interval as $\Delta t=0.005$, the number of initial points
as $n=1200$, and the number of data points as $M=1000$. We estimated the phase
function for the vdP oscillator under the observation noise, given by Gaussian
noise with mean $0$ and standard deviation $5\times 10^{-3}$ to the time
series as the observation noise. The data used for the estimation is shown in
Fig. 6(a). The Floquet exponents are estimated with $L=1000$ and
$t_{k+1}-t_{k}=0.25$ in Eq. (10). The natural frequency is estimated as
$\omega\simeq 0.9434$ and the non-zero Floquet exponent is estimated as
$\lambda\simeq-1.0885$. The latter value agree well with the theoretical value
$\lambda\simeq-1.0581$ that is evaluated from the monodromy matrix of the
original model Takata  _et al._ (2021). The estimated value of the Floquet
exponent with respect to the length $L$ is shown in Fig. 6(b). The maximum
degree of the polynomial is set to $p=18$.

We estimated the phase function ${\Theta}$ by the proposed method (Sec. III.2)
from the time series data (Fig. 6(a)). The estimated result (Fig. 6(c)) is
compared with the true phase function (Fig. 6(d)) obtained directly from the
mathematical model, Eq. (38). As in the case with the SL oscillator, the
estimated function is close to the true one except for the central and
peripheral regions far from the limit cycle.

We examined the estimation accuracy of the PSFs and PRFs. Figure 7 compares
the estimated PSFs $Z_{1},$ and $Z_{2}$ with the true ones. The estimated PSFs
agree well with the true PSFs, implying the accurate estimation of the phase
function in the vicinity of the limit cycle. Figure 8 compares the estimated
nPRFs $G_{\varsigma,j}(\theta)$ for finite-intensity impulses with the true
results. Again, the estimated nPRFs agree well with the true results. Overall,
the proposed method can estimate the PSFs and nPRFs accurately (Table 4):
$R_{Z}^{2}\geq 0.9896$ and $R_{G}^{2}\geq 0.9892$ for PSF and nPRFs,
respectively. These results indicate that the proposed method accurately
estimates the phase function in the moderately nonlinear regime around the
limit cycle.

Next, we estimated the amplitude function $R$ by the proposed method (Sec.
III.3) from the same data (Fig. 6(a)). The hyper parameter $\gamma$ is
determined to be $\gamma=1.0\times 10^{3}$ from the L-curve. The estimation
result (Fig. 6(e)) is compared with the true amplitude function (Fig. 6(f))
obtained directly from the mathematical model, Eq. (38). As in the SL case,
the estimate agrees with the true function near the limit cycle, but a
discrepancy arises in the central and peripheral regions far from the limit
cycle.

We examined the estimation accuracy of the ASFs and ARFs. Figure 9 compares
the estimated ASFs with the true ASFs, showing a reasonable agreement. Figure
10 compares the estimated nARFs with the true nARFs for finite-intensity
impulses. The estimates also agree well with the true nARFs. Again, the
proposed method can estimate the ASF and nARFs accurately (Table 4):
$R_{I}^{2}\geq 0.9792$ and $R_{H}^{2}\geq 0.9694$ for ASF and nARFs,
respectively. These results suggest that the proposed method estimates the
amplitude function around the limit cycle accurately.

Figure 6:  Phase and amplitude functions of the vdP oscillator. (a) Time
series data of the vdP oscillator used for the estimation. (b) Estimated non-
zero Floquet exponent. (c) Estimated phase function. (d) True phase function.
(e) Estimated amplitude function. (f) True amplitude function. In (c)-(f), the
red curve shows the limit cycle.  Figure 7:  Phase sensitivity functions
(PSFs) of the vdP oscillator. Each graph shows the estimated PSF (red solid
curve) and the true PSF (black dotted curve). (a): $x_{1}$ component. (b):
$x_{2}$ component. Estimated PSFs are obtained by numerical differentiation
with $\epsilon=10^{-5}$. True PSFs are calculated by solving the adjoint
equation.  Figure 8:  Normalized phase response functions (nPRFs) of the vdP
oscillator. Each graph shows the estimated nPRF (red solid curve) and the true
nPRF (black dotted curve). (a,c): Impulses applied in the $x_{1}$ direction.
(a) $\varsigma=-0.4$, (c) $\varsigma=+0.4$. (b,d): Impulses applied in the
$x_{2}$ direction. (b) $\varsigma=-0.4$, (d) $\varsigma=-0.4$.  Figure 9:
Amplitude sensitivity functions (ASFs) of the vdP oscillator. Each graph shows
the estimated ASF (red solid curve) and the true ASF (black dotted curve).
(a): $x_{1}$ component. (b): $x_{2}$ component. Estimated ASFs are obtained by
numerical differentiation with $\epsilon=10^{-5}$. True ASFs are calculated by
solving the adjoint equation.  Figure 10:  Normalized amplitude response
functions (nARFs) of the vdP oscillator. Each graph shows the estimated nARF
(red solid curve) and the true nARF (black dotted curve). (a,c): Impulses
applied in the $x_{1}$ direction. (a) $\varsigma=-0.2$, (c) $\varsigma=+0.2$.
(b,d): Impulses applied in the $x_{2}$ direction. (b) $\varsigma=-0.2$, (d)
$\varsigma=+0.2$.  Table 3: Coefficients of determination of the nPRFs and
PSFs of the vdP oscillator.

$Z_{1}$ | $Z_{2}$ | $G_{-0.4,1}$ | $G_{-0.4,2}$ | $G_{0.4,1}$ | $G_{0.4,2}$
---|---|---|---|---|---
$0.9971$ | $0.9892$ | $0.9980$ | $0.9896$ | $0.9979$ | $0.9914$

Table 4: Coefficients of determination of the nARFs and ASFs of the vdP
oscillator.

$I_{1}$ | $I_{2}$ | $H_{-0.2,1}$ | $H_{-0.2,2}$ | $H_{0.2,1}$ | $H_{0.2,2}$
---|---|---|---|---|---
$0.9795$ | $0.9792$ | $0.9753$ | $0.9694$ | $0.9707$ | $0.9736$

##  V Data-driven optimal entrainment with amplitude suppression

In this section, we apply the proposed method to the optimal entrainment of
the oscillator with amplitude suppression using the PSF and ASF developed in
Ref. Takata  _et al._ (2021). We consider a limit-cycle oscillator driven by
an external periodic input. The natural frequency of the oscillator is
$\omega$ and the frequency of the periodic input is $\Omega$, where $\Omega$
is close to $\omega$. When the input is sufficiently weak, the reduced and
averaged phase equation for the phase difference $\phi$ between the oscillator
and the external periodic input is given by

$\displaystyle\frac{d}{dt}\phi=\Delta+\Gamma(\phi),\quad\Gamma(\phi)=\left[\bm{Z}(\phi+\Omega
t)\cdot\bm{q}(\Omega t)\right]_{t},$ (43)

where $\Delta=\omega-\Omega$ is the frequency mismatch between the periodic
input and the oscillator, $\Gamma$ is the phase coupling function, $\bm{q}$ is
the periodic input, and
$\left[f(t)\right]_{t}=(\Omega/2\pi)\int_{0}^{2\pi/\Omega}f(s)ds$ denotes the
average of $f$ over one period of the external input. The optimal periodic
input $\bm{q}$ can be obtained by solving the following optimization problem
as discussed in Ref. Zlotnik  _et al._ (2013):

$\displaystyle\begin{aligned} &&\bm{q}={}&\mathop{\rm
argmax}\limits_{\bm{q}}\;-\Gamma^{\prime}(\phi^{*})\\\
&&&\mathrm{s.t.}\quad\Delta+\Gamma(\phi^{*})=0,\;\left[\|\bm{q}\|^{2}\right]_{t}=Q,\end{aligned}$
(44)

where $\Gamma^{\prime}$ is the derivative of the phase coupling function at
the stable phase locking point $\phi^{*}$; the first constraint is for the
system to have a phase-locked solution at $\phi^{*}$, and the second
constraint is for the power of the periodic input to be $Q>0$. However, when
the periodic input is not sufficiently weak, the trajectory deviates from the
limit cycle and the above method based on the phase reduction may not work
appropriately.

In the amplitude-penalty method Takata  _et al._ (2021), the amplitude
equation is also used to suppress the deviations of the trajectory from the
limit cycle. It is formulated by adding a penalty term characterizing the
effect of the input on the amplitude variable of the oscillator to the
optimization problem (44) as follows:

$\displaystyle\begin{aligned} &&\bm{q}={}&\mathop{\rm
argmax}\limits_{\bm{q}}\;-\Gamma^{\prime}(\phi^{*})-k\left[\left(\bm{I}(\phi^{*}+\Omega
t)\cdot\bm{q}(\Omega t)\right)^{2}\right]_{t}\\\
&&&\mathrm{s.t.}\;\Delta+\Gamma(\phi^{*})=0,\;\left[\|\bm{q}\|^{2}\right]_{t}=Q,\end{aligned}$
(45)

where the second term in the objective function represents the penalty for the
amplitude deviation and $k>0$ is the weight of the penalty. By solving the
modified optimization problem (45), the optimal input $\bm{q}$ that
simultaneously improves the stability of the phase-locked solution $\phi^{*}$
and suppresses the deviation of the amplitude from the limit cycle can be
obtained.

The results of the amplitude penalty method for the vdP oscillator using the
PSF and ASF estimated under the observation noise are shown in Fig. 11. The
parameters are set to $k=2\times 10^{4},Q=0.05$, $\phi^{*}=0$, and $\Delta=0$.

Figure 11(a) shows the evolution of the phase differences with and without the
amplitude penalty, where results using the PSF and ASF calculated by the
adjoint method (theoretical) and those estimated from the observed data (data-
driven) are compared. The amplitude penalty method leads to the correct phase-
locking point $\phi^{*}=0$, where the results using both PSFs and ASFs are
almost indistinguishable, while the conventional method without amplitude
suppression leads to incorrect phase-locking point different from $\phi^{*}=0$
due to amplitude deviations.

Figure 11(b) shows the trajectories driven by the periodic inputs with and
without amplitude penalty, where the estimated PSF and ASF are used. It can be
seen that the trajectory driven by the periodic input with amplitude penalty
is almost the same as the limit cycle without perturbation, while the
trajectory without amplitude penalty deviates from the limit cycle. This
result suggests that the proposed method is a promising approach to realize
the optimal entrainment of the oscillator in a data-driven manner without the
knowledge of the mathematical model.

Figure 11:  Data-driven optimal entrainment of the vdP oscillator. (a)
Evolution of the phase differences for the case with periodic input obtained
from the estimated PSF and ASF (red solid), for the case with input obtained
from true PSF and ASF (green dotted), and for the case without amplitude
suppression (blue dashed). (b) Trajectories of the oscillator state; the case
with periodic input obtained from the estimated PSF and ASF (red solid) and
the case without the amplitude penalty (blue dashed) are compared. Black
dotted curve shows the limit cycle without periodic input.

##  VI Concluding remarks

In this study, we proposed a method for estimating the phase and amplitude
functions of limit-cycle oscillators only from observed time series. Our
method is based on the polynomial regression, which is simple and
computationally efficient. The resulting optimization problems are convex and
can be easily solved. The numerical results based on two limit-cycle
oscillators (SL and vdP oscillator) showed that our method can estimate the
phase and amplitude function accurately around the limit cycle under the
observation noise. Furthermore, we demonstrated that the proposed method
enables us to achieve the optimal entrainment with amplitude suppression in a
data-driven manner, and the performance was comparable to the model-based
case.

The proposed method cannot accurately estimate the phase and amplitude
functions in the phase-plane regions away from the limit cycle (Fig. 1 and 6).
In particular, it cannot estimate the singular behavior of the amplitude
function near the fixed point of the system. However, our method accurately
estimates the PSF and nPRFs (Figs. 2, 3, 7, and 8) and ASF and nARFs (Figs. 4,
5, 9, and 10). In the phase-amplitude description, the PSF, nPRFs, ASF, and
nARFs are essential functions for the phase-amplitude description of limit-
cycle oscillators. Thus, the discrepancy of the phase and amplitude functions
away from the limit cycle is not a serious problem for applications. Indeed,
we demonstrated that the proposed method is useful for data-driven optimal
entrainment of the limit-cycle oscillators (Fig. 11). The proposed method is
potentially helpful also for data-driven analysis and control of limit-cycle
oscillators subjected to pulsatile forcing or coupling Ermentrout and Terman
(2010); Nakao  _et al._ (2005); Arai and Nakao (2008).

As mentioned in the introduction, a major alternative to our proposed method
is the EDMD Williams  _et al._ (2015), which estimates Koopman eigenvalues and
eigenfunctions from the observed data. The asymptotic phase function $\Theta$
and amplitude function $R$ are related to the Koopman eigenfunctions with
Koopman eigenvalues $i\omega$ and $\lambda$ respectively, namely, they satisfy
$(d/dt)e^{i\Theta}=i\omega e^{i\Theta}$ and $(d/dt)R=\lambda R$ Mauroy and
Mezić (2016, 2018); Mauroy  _et al._ (2020, 2013); Shirasaka  _et al._ (2017).
The primary difference of the present study from those based on EDMD is that
the natural frequency $\omega$ and the Floquet exponent $\lambda$ are given as
external parameters (measured separately) and then the phase and amplitude
functions are estimated, while in the EDMD, the frequency and Floquet
exponents are estimated from the time series together with the eigenfunctions.
Therefore, our method can be more stable, because $\omega$ and $\lambda$ are
measured separately by using the classical, well-established method for the
estimation of Lyapunov exponents Wolf  _et al._ (1985). We thus think that our
proposed method is complementary to the EDMD methods in the analysis of real-
world oscillatory signals.

In this study, we applied the proposed method to limit-cycle oscillators with
known mathematical models in order to evaluate its accuracy. Our future work
is to apply the method to real-world observed data such as ECGs Kralemann  _et
al._ (2013). To this end, we plan to further develop a method of preprocessing
for the time series to mitigate the effect of stronger noise and non-
stationarity.

###### Acknowledgements.

H. N. thanks financial support from JSPS KAKENHI (Nos. JP17H03279, JP18H03287,
and JPJSBP120202201) and JST CREST (No. JP-MJCR1913). R. K. thanks JSPS
KAKENHI (Nos. 18K11560, 19H01133, 21H03559, 21H04571), JST PRESTO (No.
JPMJPR1925), and AMED (No. JP21wm0525004).

## References

  * Stankovski  _et al._ (2015) Tomislav Stankovski, Valentina Ticcinelli, Peter VE McClintock, and Aneta Stefanovska, “Coupling functions in networks of oscillators,” New Journal of Physics 17, 035002 (2015).
  * Stankovski  _et al._ (2017a) Tomislav Stankovski, Valentina Ticcinelli, Peter VE McClintock, and Aneta Stefanovska, “Neural cross-frequency coupling functions,” Frontiers in Systems Neuroscience 11, 33 (2017a).
  * Borgius  _et al._ (2014) Lotta Borgius, Hiroshi Nishimaru, Vanessa Caldeira, Yuka Kunugise, Peter Löw, Ramon Reig, Shigeyoshi Itohara, Takuji Iwasato, and Ole Kiehn, “Spinal glutamatergic neurons defined by EphA4 signaling are essential components of normal locomotor circuits,” Journal of Neuroscience 34, 3841–3853 (2014).
  * Collins and Stewart (1993) J. J. Collins and I. N. Stewart, “Coupled nonlinear oscillators and the symmetries of animal gaits,” Journal of Nonlinear Science 3, 349–392 (1993).
  * Kobayashi  _et al._ (2016) Ryota Kobayashi, Hiroshi Nishimaru, and Hisao Nishijo, “Estimation of excitatory and inhibitory synaptic conductance variations in motoneurons during locomotor-like rhythmic activity,” Neuroscience 335, 72–81 (2016).
  * Funato  _et al._ (2016) Tetsuro Funato, Yuki Yamamoto, Shinya Aoi, Takashi Imai, Toshio Aoyagi, Nozomi Tomita, and Kazuo Tsuchiya, “Evaluation of the phase-dependent rhythm control of human walking using phase response curves,” PLoS Computational Biology 12, e1004950 (2016).
  * Kralemann  _et al._ (2013) Björn Kralemann, Matthias Frühwirth, Arkady Pikovsky, Michael Rosenblum, Thomas Kenner, Jochen Schaefer, and Maximilian Moser, “In vivo cardiac phase response curve elucidates human respiratory heart rate variability,” Nature Communications 4, 2418 (2013).
  * Hackel (2007) Matthias Hackel,  _Humanoid Robots_ (IntechOpen, Rijeka, 2007).
  * Garcia  _et al._ (1998) Mariano Garcia, Anindya Chatterjee, Andy Ruina, and Michael Coleman, “The simplest walking model: Stability, complexity, and scaling,” Journal of Biomechanical Engineering 120, 281–288 (1998).
  * Strogatz (2015) Steven H. Strogatz,  _Nonlinear Dynamics and Chaos_ (CRC Press, Boca Raton, 2015).
  * Kuramoto (1984) Yoshiki Kuramoto,  _Chemical Oscillations, Waves, and Turbulence_ (Springer, Berlin, 1984).
  * Hoppensteadt and Izhikevich (1997) Frank C. Hoppensteadt and Eugene M. Izhikevich,  _Weakly Connected Neural Networks_ (Springer, New York, 1997).
  * Winfree (2001) Arthur T. Winfree,  _The Geometry of Biological Time_ (Springer, New York, 2001).
  * Ermentrout and Terman (2010) G. Bard Ermentrout and David H. Terman,  _Mathematical Foundations of Neuroscience_ (Springer, New York, 2010).
  * Nakao (2016) Hiroya Nakao, “Phase reduction approach to synchronisation of nonlinear oscillators,” Contemporary Physics 57, 188–214 (2016).
  * Monga  _et al._ (2019) Bharat Monga, Dan Wilson, Tim Matchen, and Jeff Moehlis, “Phase reduction and phase-based optimal control for biological systems: a tutorial,” Biological cybernetics 113, 11–46 (2019).
  * Kuramoto and Nakao (2019) Yoshiki Kuramoto and Hiroya Nakao, “On the concept of dynamical reduction: the case of coupled oscillators,” Philosophical Transactions of the Royal Society A 377, 20190041 (2019).
  * Ermentrout  _et al._ (2019) Bard Ermentrout, Youngmin Park, and Dan Wilson, “Recent advances in coupled oscillator theory,” Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 377, 20190092 (2019).
  * Wedgwood  _et al._ (2013) Kyle CA Wedgwood, Kevin K Lin, Ruediger Thul, and Stephen Coombes, “Phase-amplitude descriptions of neural oscillator models,” The Journal of Mathematical Neuroscience 3, 1–22 (2013).
  * Wilson and Moehlis (2016) Dan Wilson and Jeff Moehlis, “Isostable reduction of periodic orbits,” Phys. Rev. E 94, 052213 (2016).
  * Mauroy and Mezić (2016) Alexandre Mauroy and Igor Mezić, “Global stability analysis using the eigenfunctions of the Koopman operator,” IEEE Transactions on Automatic Control 61, 3356–3369 (2016).
  * Shirasaka  _et al._ (2017) Sho Shirasaka, Wataru Kurebayashi, and Hiroya Nakao, “Phase-amplitude reduction of transient dynamics far from attractors for limit-cycling systems,” Chaos: An Interdisciplinary Journal of Nonlinear Science 27, 023119 (2017).
  * Mauroy and Mezić (2018) Alexandre Mauroy and Igor Mezić, “Global computation of phase-amplitude reduction for limit-cycle dynamics,” Chaos: An Interdisciplinary Journal of Nonlinear Science 28, 073108 (2018).
  * Kotani  _et al._ (2020) Kiyoshi Kotani, Yutaro Ogawa, Sho Shirasaka, Akihiko Akao, Yasuhiko Jimbo, and Hiroya Nakao, “Nonlinear phase-amplitude reduction of delay-induced oscillations,” Physical Review Research 2, 033106 (2020).
  * Shirasaka  _et al._ (2020) Sho Shirasaka, Wataru Kurebayashi, and Hiroya Nakao, “Phase-amplitude reduction of limit cycling systems,” in  _The Koopman Operator in Systems and Control_ (Springer, 2020) pp. 383–417.
  * Nakao (2021) Hiroya Nakao, “Phase and amplitude description of complex oscillatory patterns in reaction-diffusion systems,” in  _Physics of Biological Oscillators: New Insights into Non-Equilibrium and Non-Autonomous Systems_ , edited by Aneta Stefanovska and Peter V. E. McClintock (Springer International Publishing, Cham, 2021) pp. 11–27.
  * Takata  _et al._ (2021) Shohei Takata, Yuzuru Kato, and Hiroya Nakao, “Fast optimal entrainment of limit-cycle oscillators by strong periodic inputs via phase-amplitude reduction and Floquet theory,” Chaos: An Interdisciplinary Journal of Nonlinear Science 31, 093124 (2021).
  * Monga and Moehlis (2019) Bharat Monga and Jeff Moehlis, “Optimal phase control of biological oscillators using augmented phase reduction,” Biological Cybernetics 113, 161–178 (2019).
  * Schmid (2010) Peter J. Schmid, “Dynamic mode decomposition of numerical and experimental data,” ‘ 656, 5–28 (2010).
  * Kutz  _et al._ (2016) J. Nathan Kutz, Steven L. Brunton, Bingni W. Brunton, and Joshua L. Proctor,  _Dynamic mode decomposition: data-driven modeling of complex systems_ (SIAM, 2016).
  * Williams  _et al._ (2015) Matthew O. Williams, Ioannis G. Kevrekidis, and Clarence W. Rowley, “A data-driven approximation of the Koopman operator: Extending dynamic mode decomposition,” Journal of Nonlinear Science 25, 1307–1346 (2015).
  * Lusch  _et al._ (2018) Bethany Lusch, J. Nathan Kutz, and Steven L. Brunton, “Deep learning for universal linear embeddings of nonlinear dynamics,” Nature communications 9, 1–10 (2018).
  * Wilson (2020) Dan Wilson, “A data-driven phase and isostable reduced modeling framework for oscillatory dynamical systems,” Chaos: An Interdisciplinary Journal of Nonlinear Science 30, 013121 (2020).
  * Wilson (2021) Dan Wilson, “Data-driven inference of high-accuracy isostable-based dynamical models in response to external inputs,” Chaos: An Interdisciplinary Journal of Nonlinear Science 31, 063137 (2021).
  * Klus  _et al._ (2020) Stefan Klus, Feliks Nüske, Sebastian Peitz, Jan-Hendrik Niemann, Cecilia Clementi, and Christof Schütte, “Data-driven approximation of the Koopman generator: Model reduction, system identification, and control,” Physica D: Nonlinear Phenomena 406, 132416 (2020).
  * Li  _et al._ (2017) Qianxiao Li, Felix Dietrich, Erik M. Bollt, and Ioannis G. Kevrekidis, “Extended dynamic mode decomposition with dictionary learning: A data-driven adaptive spectral decomposition of the Koopman operator,” Chaos: An Interdisciplinary Journal of Nonlinear Science 27, 103111 (2017).
  * Takeishi  _et al._ (2017) Naoya Takeishi, Yoshinobu Kawahara, and Takehisa Yairi, “Learning Koopman invariant subspaces for dynamic mode decomposition,” Advances in Neural Information Processing Systems 30 (2017).
  * Alford-Lago  _et al._ (2021) Daniel J Alford-Lago, Christopher W Curtis, Alexander Ihler, and Opal Issan, “Deep learning enhanced dynamic mode decomposition,” arXiv preprint arXiv:2108.04433 (2021).
  * Terao  _et al._ (2021) Hiroaki Terao, Sho Shirasaka, and Hideyuki Suzuki, “Extended dynamic mode decomposition with dictionary learning using neural ordinary differential equations,” Nonlinear Theory and Its Applications, IEICE 12, 626–638 (2021).
  * Gulina and Mauroy (2021) Marvyn Gulina and Alexandre Mauroy, “Two methods to approximate the Koopman operator with a reservoir computer,” Chaos: An Interdisciplinary Journal of Nonlinear Science 31, 023116 (2021).
  * Kevrekidis  _et al._ (2016) I Kevrekidis, Clarence W Rowley, and M Williams, “A kernel-based method for data-driven Koopman spectral analysis,” Journal of Computational Dynamics 2, 247–265 (2016).
  * Héas  _et al._ (2020) Patrick Héas, Cédric Herzet, and Benoit Combes, “Generalized kernel-based dynamic mode decomposition,” in  _ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ (IEEE, 2020) pp. 3877–3881.
  * Mauroy  _et al._ (2020) Alexandre Mauroy, Igor Mezić, and Yoshihiko Susuki,  _The Koopman Operator in Systems and Control: Concepts, Methodologies, and Applications_ , Vol. 484 (Springer Nature, 2020).
  * Mauroy  _et al._ (2013) Alexandre Mauroy, Igor Mezić, and Jeff Moehlis, “Isostables, isochrons, and Koopman spectrum for the action-angle representation of stable fixed point dynamics,” Physica D: Nonlinear Phenomena 261, 19–30 (2013).
  * Galán  _et al._ (2005) Roberto F. Galán, G. Bard Ermentrout, and Nathaniel N. Urban, “Efficient estimation of phase-resetting curves in real neurons and its significance for neural-network modeling,” Phys. Rev. Lett. 94, 158101 (2005).
  * Ota  _et al._ (2009) Kaiichiro Ota, Masaki Nomura, and Toshio Aoyagi, “Weighted spike-triggered average of a fluctuating stimulus yielding the phase response curve,” Phys. Rev. Lett. 103, 024101 (2009).
  * Aonishi and Ota (2006) Toru Aonishi and Keisuke Ota, “Statistical estimation algorithm for phase response curves,” Journal of the Physical Society of Japan 75, 114802 (2006).
  * Imai  _et al._ (2017) Takashi Imai, Kaiichiro Ota, and Toshio Aoyagi, “Robust measurements of phase response curves realized via multicycle weighted spike-triggered averages,” Journal of the Physical Society of Japan 86, 024009 (2017).
  * Cestnik and Rosenblum (2018) Rok Cestnik and Michael Rosenblum, “Inferring the phase response curve from observation of a continuously perturbed oscillator,” Scientific reports 8, 1–10 (2018).
  * Nakae  _et al._ (2010) Ken Nakae, Yukito Iba, Yasuhiro Tsubo, Tomoki Fukai, and Toshio Aoyagi, “Bayesian estimation of phase response curves,” Neural Networks 23, 752–763 (2010).
  * Ota  _et al._ (2011) Keisuke Ota, Toshiaki Omori, Shigeo Watanabe, Hiroyoshi Miyakawa, Masato Okada, and Toru Aonishi, “Measurement of infinitesimal phase response curves from noisy real neurons,” Physical Review E 84, 041902 (2011).
  * Tokuda  _et al._ (2007) Isao T Tokuda, Swati Jain, István Z Kiss, and John L Hudson, “Inferring phase equations from multivariate time series,” Physical review letters 99, 064101 (2007).
  * Kralemann  _et al._ (2008) Björn Kralemann, Laura Cimponeriu, Michael Rosenblum, Arkady Pikovsky, and Ralf Mrowka, “Phase dynamics of coupled oscillators reconstructed from data,” Physical Review E 77, 066205 (2008).
  * Duggento  _et al._ (2012) Andrea Duggento, Tomislav Stankovski, Peter VE McClintock, and Aneta Stefanovska, “Dynamical Bayesian inference of time-evolving interactions: From a pair of coupled oscillators to networks of oscillators,” Physical Review E 86, 061126 (2012).
  * Stankovski  _et al._ (2017b) Tomislav Stankovski, Tiago Pereira, Peter VE McClintock, and Aneta Stefanovska, “Coupling functions: universal insights into dynamical interaction mechanisms,” Reviews of Modern Physics 89, 045001 (2017b).
  * Ota and Aoyagi (2014) Kaiichiro Ota and Toshio Aoyagi, “Direct extraction of phase dynamics from fluctuating rhythmic data based on a Bayesian approach,” arXiv preprint arXiv:1405.4126 (2014).
  * Onojima  _et al._ (2018) Takayuki Onojima, Takahiro Goto, Hiroaki Mizuhara, and Toshio Aoyagi, “A dynamical systems approach for estimating phase interactions between rhythms of different frequencies from experimental data,” PLoS computational biology 14, e1005928 (2018).
  * Suzuki  _et al._ (2018) Kento Suzuki, Toshio Aoyagi, and Katsunori Kitano, “Bayesian estimation of phase dynamics based on partially sampled spikes generated by realistic model neurons,” Frontiers in Computational Neuroscience 11, 116 (2018).
  * Arai  _et al._ (2021) Takahiro Arai, Yoji Kawamura, and Toshio Aoyagi, “Extracting phase coupling functions between networks of dynamical elements that exhibit collective oscillations: Direct extraction from time-series data,” arXiv preprint arXiv:2110.15598 (2021).
  * Brown  _et al._ (2004) Eric Brown, Jeff Moehlis, and Philip Holmes, “On the phase reduction and response dynamics of neural oscillator populations,” Neural Computation 16, 673–715 (2004).
  * Nakao  _et al._ (2005) Hiroya Nakao, Ken-suke Arai, Ken Nagai, Yasuhiro Tsubo, and Yoshiki Kuramoto, “Synchrony of limit-cycle oscillators induced by random external impulses,” Phys. Rev. E 72, 026220 (2005).
  * Arai and Nakao (2008) Kensuke Arai and Hiroya Nakao, “Phase coherence in an ensemble of uncoupled limit-cycle oscillators receiving common poisson impulses,” Physical Review E 77, 036218 (2008).
  * Wolf  _et al._ (1985) Alan Wolf, Jack B. Swift, Harry L. Swinney, and John A. Vastano, “Determining Lyapunov exponents from a time series,” Physica D: Nonlinear Phenomena 16, 285–317 (1985).
  * Knowles and Renka (2014) Ian Knowles and Robert J. Renka, “Methods for numerical differentiation of noisy data,” Electron. J. Differ. Equ 21, 235–246 (2014).
  * Hansen (1992) Per Christian Hansen, “Analysis of discrete ill-posed problems by means of the L-curve,” SIAM Review 34, 561–580 (1992).
  * van der Pol (1927) Balth van der Pol, “Forced oscillations in a circuit with non-linear resistance.(Reception with reactive triode),” The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 3, 65–80 (1927).
  * van der Pol and van der Mark (1927) Balth van der Pol and Jan van der Mark, “Frequency demultiplication,” Nature 120, 363–364 (1927).
  * van der Pol (1926) Balth van der Pol, “On relaxation-oscillations,” The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 2, 978–992 (1926).
  * Kato  _et al._ (2021) Yuzuru Kato, Jinjie Zhu, Wataru Kurebayashi, and Hiroya Nakao, “Asymptotic phase and amplitude for classical and semiclassical stochastic oscillators via Koopman operator theory,” Mathematics 9, 2188 (2021).
  * Zlotnik  _et al._ (2013) Anatoly Zlotnik, Yifei Chen, István Z. Kiss, Hisa-Aki Tanaka, and Jr-Shin Li, “Optimal waveform for fast entrainment of weakly forced nonlinear oscillators,” Phys. Rev. Lett. 111, 024102 (2013).