# Reconstruction of phase-amplitude dynamics from electrophysiological signals

Azamat Yeldesbay  University of Cologne
Institute of Zoology
Cologne, Germany  Cognitive Neuroscience
Institute of Neuroscience and Medicine (INM-3)
Research Centre Jülich
Jülich, Germany  Gemma Huguet  Universitat Politécnica de Catalunya
Departament Matemátiques
Barcelona, Spain  Centre de Recerca Matemática
Barcelona, Spain  Silvia Daun  Cognitive Neuroscience
Institute of Neuroscience and Medicine (INM-3)
Research Centre Jülich
Jülich, Germany  University of Cologne
Institute of Zoology
Cologne, Germany

###### Abstract

Signals from interacting brain regions display transient synchronization of
phases and amplitudes in different frequencies. Commonly, the interaction
between regions of the brain is quantitatively described by either analyzing
the correlations of amplitudes of the measured signals or by calculating
phase-synchronization measures. However, for a complete picture of the
interactions it is important to analyze the dynamics of both the amplitude and
the phase.

In this work, we present a new method for finding the coupling between brain
regions by reconstructing the phase-amplitude dynamics directly from the
measured electrophysiological signals. For this purpose, we use the recent
advances in the field of phase-amplitude reduction of oscillatory systems,
which allow the representation of an uncoupled oscillatory system as a phase-
amplitude oscillator in a unique form using transformations (parametrizations)
related to the eigenfunctions of the Koopman operator. By combining the
parametrization method and the Fourier-Laplace averaging method of finding the
eigenfunctions of the Koopman operator, we developed a novel method of
assessing the transformation functions from the signals of the interacting
oscillatory systems. The resulting reconstructed dynamical system is a network
of phase-amplitude oscillators with the interactions between them represented
as coupling functions in phase and amplitude coordinates.

Using synthetic signals generated from several models with known and unknown
theoretical phase-amplitude reduced forms, we demonstrate that our method is
capable of finding the proper unique dynamic form of these oscillatory systems
in the reduced phase-amplitude space.

Our method can be applied to describe any network of interacting oscillators
as a dynamical system using signals of the network elements. In particular, to
analyze the coupling between distant brain regions using high time resolution
signals, such as electroencephalography (EEG) or magnetoencephalography (MEG).
Further simulation and study of the reconstructed dynamical system then
enables the construction and investigation of a mathematical model of various
neural pathologies and disorders of the brain.

##### keywords

Oscillators, synchronization, experimental data, network, communication
between brain regions

##  1 Introduction

Performing a cognitive task involves large-scale integration of multiple brain
regions, which is accompanied by an emergence of transient synchronization
over multiple frequency bands [1, 2, 3]. The transient synchronization
observed during the stimulus of a cognitive task, e.g. the processing of a
movement, its preparation and execution stages, is functionally relevant for
the activation and coordination of a widespread network of oscillatory neural
populations across brain areas [4, 2]. A detailed study of the emergence
mechanisms of the transient synchronization is important to understand the
interaction between brain areas involved in such a task.

An assessment of the coupling between distant brain regions by means of the
synchronization of oscillatory activity observed in EEG/MEG signals can be
performed in different ways [5, 6]: by considering the relation between the
phase of the signals and calculating phase synchronization measures [7, 8, 9],
by analyzing the relationship between amplitudes with respect to an event
using event-related synchronization (ERS) and desynchronization (ERD) measures
[10], or by calculating phase-amplitude coupling (PAC) measures [11, 12, 13],
which consider the correlation between phase and amplitude of different
frequency bands. All these conventional measures provide only an estimation of
the interaction magnitude between the oscillatory sources, but do not provide
information about the direction of the coupling and can not explain the
mechanism of the transient synchronization. The latter is possible by
analyzing a dynamical system that describes the oscillatory network and
reproduces the observed transient synchronization.

The reconstruction of a dynamical system in a form of coupled phase
oscillators from observations was first presented in [14] and was developed
further in the following works [15, 16, 17, 18, 19, 20, 21, 22, 23]. These
methods are based on the phase reduction approach [24, 25, 26], which is
usually limited by weak coupling between oscillators, which, however, is
usually not the case for the connection between brain regions. For strong
coupling cases consideration of the amplitude of the oscillations is required.
In the recent decade several approaches of phase-amplitude reduction were
developed [27], namely methods based on the parametrization approach [28, 29,
30, 31, 32, 33, 34], the Koopman operator [35, 36, 37, 38, 39] and the
isostable coordinates approaches [40, 41, 42, 43, 44], with further
development as data-driven phase-amplitude reduction approaches using
isostable coordinates [45, 46, 47, 48, 49], using dynamic mode decomposition
(DMD) [50, 51, 52, 53, 54], and as an extension of dynamic causal modelling
[55].

In this work we aim at developing a new method to build a dynamical system
that will be able to describe the coupling between regions of the brain (the
sources of measured signals) in the framework of transient synchronization,
reflected in phase and amplitude of the measured signals (Fig. 1). We achieve
this by reconstructing the network of dynamical systems in a unique
(invariant) form using the amplitude and phase of the measured signals by
adapting the recent advances in the phase-amplitude reduction theory. In
particular, the phase-amplitude reduction performs a transformation of
observable variables (i.e. the phase and amplitude of the measured signals
$(\theta_{i},r_{i})$, see Fig. 1) into uniquely defined reduced variables (the
reduced amplitude and phase $(\varphi_{i},\sigma_{i})$, see Fig. 1). This
reduction allows us to find a unique dynamical description of the network and
the coupling terms between the regions. Thereby, we make several assumptions:
(1) during the performance of a cognitive task some intrinsic processes within
the cortical column of a brain region activate oscillatory activities in
different frequencies, (2) in the considered interval of time this oscillating
system has a limit cycle and the amplitude and phase are determinable around
this limit cycle, (3) the interaction between brain areas is reflected in the
transient synchronization of phase and amplitude of the measured signals.

Figure 1: From electrophysiological signals to a network model. We extract the
amplitude and phase of the brain signals of different sources obtained from
experimental measurements during the performance of a cognitive task. The
extracted (observable) amplitude and phase $(\theta_{i},r_{i})$ for each
region $i$ are transformed into the reduced amplitude and phase
$(\varphi_{i},\sigma_{i})$ by means of the phase-amplitude reduction. The
reconstructed network is a dynamical system, where the edges are given by
coupling functions. The parameters of the reconstructed system and of the
transformation from observable to reduced space are obtained from the
experimental data.

The outline of this article is the following. We start with the introduction
of main concepts of the phase-amplitude reduction and split the methods part
into two sub-sections. First, we a give brief overview of the theoretical
methods of phase-amplitude reduction of an uncoupled oscillator with known
vector filed (VF). There we introduce the parametrization and the Koopman
operator approaches. Further, based on these approaches we present a numerical
method of finding the transformation functions that perform phase-amplitude
reduction for a system with known VF. In the second part of the methods we
introduce a method of reconstruction of the phase-amplitude dynamics of the
observed system from measured signals, where we approximate the VF of the
system from measurements and find the coupling functions as the result of the
phase-amplitude reduction. In the results part we demonstrate the
applicability of the phase-amplitude reduction method to several models of
inter-connected oscillators with different amplitude and phase properties.
Finally, a discussion and an overview of further general application
possibilities are given in the last section.

##  2 Methods

###  2.1 Introduction to phase-amplitude reduction

In this section, we introduce the amplitude and phase reduction method and
review several concepts that will form the basis for the methods presented
further in this study.

####  2.1.1 Koopman operator and phase-amplitude reduction

Let us consider a dynamical system described by an autonomous system of ODEs

$\dot{\mathbf{x}}=\mathbf{F(\mathbf{x})},\mathbf{x}\in\mathcal{M}\subset\mathbb{R}^{n},$
(1)

which has a hyperbolic limit cycle $\Gamma$ of period $T$, parameterized by
$\varphi=\omega t$ with $\omega=2\pi/T$ as

$\begin{array}[]{rccl}\gamma:&[0,2\pi)&\rightarrow&\mathbb{R}^{n}\\\
&\varphi&\mapsto&\gamma(\varphi)\end{array}$

such that $\gamma(\varphi+2\pi)=\gamma(\varphi)$ and
$\mathbf{x}(t)=\gamma(\omega t)$ is a solution of Eq.(1).

We introduce the phase-amplitude reduction following [39]. Let us assume that
there exists a local diffeomorphism
$\mathcal{Z}(\mathbf{x}):\mathcal{M}\subset\mathbb{R}^{n}\rightarrow\mathbb{C}^{n}$
that transforms the state space coordinates $\mathbf{x}$ into new coordinates
$\mathbf{z}$, such that the dynamics in the new coordinates are represented as

$\dot{z}_{i}=\mu_{i}z_{i},\;i\in[1,n].$ (2)

Then, we redefine the coordinates in the following way. If
$\mu_{i}\in\mathbb{R}$, we define a new coordinate
$\sigma_{i}=z_{i}\in\mathbb{R}$ and denote $\lambda_{i}\triangleq\mu_{i}$. If
$\mu_{i}\notin\mathbb{R}$ we define $\sigma_{i}=|z_{i}|\in\mathbb{R}$ and
$\varphi_{i}=\angle z_{i}\in[0,2\pi)$, and also denote
$\lambda_{i}\triangleq\mathcal{R}\\{\mu_{i}\\}$ and
$\omega_{i}\triangleq\mathcal{I}\\{\mu_{i}\\}$. By removing repeated variables
due to complex conjugate pairs, we obtain an $n$-dimensional system of ODEs in
the amplitude and angle variables given by

$\begin{array}[]{rcl}\dot{\varphi}_{k}&=&\omega_{k},\\\
\dot{\sigma}_{l}&=&\lambda_{l}\sigma_{l},\end{array}$ (3)

where $k=1,\ldots,m$, $l=1,\ldots,n-m$, and $m$ is the number of complex
conjugate $\mu$’s. This transformation of the system Eq. (1) into the linear
system system Eq. (3) is called phase-amplitude reduction [39] and corresponds
to the parametrization method presented in [30] and [31]. To perform the
phase-amplitude reduction, we should find a transformation from the state
space coordinates $\mathbf{x}\in\mathbb{R}^{n}$ to a new phase and amplitude
coordinates
$\chi:(\varphi_{1},\ldots,\varphi_{m},\sigma_{1},\ldots,\sigma_{n-m})\in[0,2\pi)^{m}\times\mathbb{R}^{n-m}$,
where the dynamics is linear as given in Eq. (3). We will refer to this space
as the reduced space.

This transformation is closely related to the Koopman operator as shown in
[37, 38, 39]. The semi-group of Koopman operators associated with system Eq.
(1) is a linear one-parameter semi-group of operators $\\{U^{t}\\}_{t\geq
0}:\mathcal{F}\rightarrow\mathcal{F}$ acting on the space $\mathcal{F}$ of
scalar functions (observables)
$g:\mathcal{M}\subset\mathbb{R}^{n}\rightarrow\mathbb{R}$ as follows [39, 56]

$U^{t}g=g\circ s^{t},$ (4)

where $s^{t}\triangleq s(t,\cdot):\mathcal{M}\rightarrow\mathcal{M}$ is the
flow of system Eq. (1), i.e $s^{t}(\mathbf{x})=s(t,\mathbf{x})$ is the
solution of system Eq. (1) with initial condition $\mathbf{x}\in\mathcal{M}$.
In other words, if the state $\mathbf{x}$ after time $t$ transfers into
$s(t,\mathbf{x})$ due to the dynamics of system Eq. (1), then the observable
$g(\mathbf{x})$ of the system transfers into $g(s(t,\mathbf{x}))$ due to the
linear operator $U^{t}$. Thus, the Koopman operator lifts the dynamics from
the state space ($\mathbf{x}$) to the observable space ($g(\mathbf{x})$) such
that the dynamics in the reduced space is linear but infinite-dimensional.

Due to linearity we can define the eigenfunctions of the Koopman operator as

$U^{t}\phi_{\mu}=e^{\mu t}\phi_{\mu},\;\forall t\geq 0,$ (5)

where $\mu\in\mathbb{C}$ is the corresponding eigenvalue [39, 56]. If we
define the new variable $z(t)=\phi_{\mu}(s(t,\mathbf{x}))$ then, using Eq. (4)
and Eq. (5), we have

$\dot{z}(t)=\frac{d}{dt}\phi_{\mu}(s(t,\mathbf{x}))=\frac{d}{dt}U^{t}\phi_{\mu}(\mathbf{x})=\frac{d}{dt}e^{\mu
t}\phi_{\mu}(\mathbf{x})\\\ =\mu e^{\mu t}\phi_{\mu}(\mathbf{x})=\mu z(t).$
(6)

Thus, the eigenfunctions of the Koopman operator transform the dynamics of
system Eq. (1) into the linear dynamics given in Eq. (2), and, thereby, are
the transformation functions that perform the phase-amplitude reduction.

In general the Koopman operator associated with the system Eq. (1) has $n$
eigenfunctions, therefore $n$ variables in the reduced phase-amplitude
coordinate space $\chi$. One of the phase variables in Eq. (3), e.g.
$\varphi_{1}$, describes the periodic dynamics along the limit cycle $\Gamma$,
whereas other phase $\varphi_{k}$ and amplitude $\sigma_{l}$ variables
describe the transient dynamics. Moreover, on the limit cycle, we have
$\sigma_{l}=0$ for all $l$. In our study we assume that the dynamics of a
considered system can be approximated as planar, meaning that all amplitude
variables except one, e.g. $\sigma_{1}$, with
$\mathcal{R}\\{\lambda_{1}\\}>\mathcal{R}\\{\lambda_{i}\\}$, decay rapidly to
the limit cycle, thereby leaving only two variables - one phase and one
amplitude variable - that describe the state of the system. Further in the
text we use these variables without indices.

The dynamics in the reduced space then reads

$\left\\{\begin{array}[]{cc}\dot{\varphi}=\omega,\\\
\dot{\sigma}=\lambda\sigma,\end{array}\right.$ (7)

where $\omega=2\pi/T$ is the radial frequency of the oscillation and $\lambda$
is the real part of the eigenvalue with smallest modulus.

It needs to be mentioned that the definition of the coordinates
$(\varphi,\sigma)$ in the reduced space is not unique. Any angular variable
$\hat{\varphi}=\varphi+\vartheta$, where $\vartheta\in[0,2\pi)$ and any
amplitude variable $\hat{\sigma}=b\sigma$, where $b\in\mathbb{R}$ also satisfy
Eq. (7). This uncertainty can be solved if we fix $\varphi=0$ anywhere on the
limit cycle, and set a condition to the gradient of the transformation
function with respect to $\sigma$, as it will be mentioned later.

####  2.1.2 Parametrization method and the invariance equation

Here we present the parametrization method which was introduced in [30] and
used in [31, 57] as a phase-amplitude reduction method. Assume we have a
planar system ($\mathbf{x}\triangleq(x,y)$) – for higher dimensions see [34] –
and we know the transformation from the phase $\varphi$ and amplitude $\sigma$
variables into the state variables $\mathbf{x}$, i.e.
$\mathbf{x}=K(\varphi,\sigma)$. Using the chain rule we obtain

$\dot{\mathbf{x}}=\frac{\partial
K}{\partial\varphi}\dot{\varphi}+\frac{\partial
K}{\partial\sigma}\dot{\sigma},$ (8)

using Eq. (1) with $n=2$ and Eq. (7) we find

$\omega\frac{\partial K}{\partial\varphi}+\lambda\sigma\frac{\partial
K}{\partial\sigma}=\mathbf{F}(K(\varphi,\sigma)).$ (9)

This equation is called the invariance equation [31]. One can also derive the
inverse invariance equations as follows

$\displaystyle F^{(x)}(\mathbf{x})\frac{\partial\Phi(\mathbf{x})}{\partial
x}+F^{(y)}(\mathbf{x})\frac{\partial\Phi(\mathbf{x})}{\partial y}=\omega,$
(10) $\displaystyle
F^{(x)}(\mathbf{x})\frac{\partial\Sigma(\mathbf{x})}{\partial
x}+F^{(y)}(\mathbf{x})\frac{\partial\Sigma(\mathbf{x})}{\partial
y}=\lambda\Sigma(\mathbf{x}),$ (11)

where $F^{(x)}$, $F^{(y)}$ are the corresponding components of $\mathbf{F}$,
$\Phi(\mathbf{x})$ and $\Sigma(\mathbf{x})$ are the phase and amplitude
transformations, i.e. from state space into the reduced space.

Note that Eqs. (9),(10) and (11) can also be expressed in polar coordinates
($\mathbf{x}\triangleq(\theta,r)$). Then, Eqs. (10) and (11) read as

$\displaystyle
F^{(\theta)}(\mathbf{x})\frac{\partial\Phi(\mathbf{x})}{\partial\theta}+F^{(r)}(\mathbf{x})\frac{\partial\Phi(\mathbf{x})}{\partial
r}=\omega,$ (12) $\displaystyle
F^{(\theta)}(\mathbf{x})\frac{\partial\Sigma(\mathbf{x})}{\partial\theta}+F^{(r)}(\mathbf{x})\frac{\partial\Sigma(\mathbf{x})}{\partial
r}=\lambda\Sigma(\mathbf{x}),$ (13)

where $F^{(r)}(\mathbf{x})$ and $F^{(\theta)}(\mathbf{x})$ are the components
of $\mathbf{F}(\mathbf{x})$ in the polar coordinates.

Figure 2: Phase-amplitude transformations. The transformation that performs
the phase-amplitude reduction from the observable space $(\theta,r)$ to the
reduced space $(\varphi,\sigma)$ is the inverse transformation. The
transformation that changes the coordinates from the reduced space to the
observable space is the forward transformation.

To further distinguish the transformations, we will use the following
notation. The transformation $K(\varphi,\sigma)$, that maps the variables
$(\varphi,\sigma)$ from the reduced space to the variables
$\mathbf{x}=(\theta,r)$ of the state space (see Fig. 2) is called forward
transformation. We use superscripts to indicate its components, for instance
$K^{(x)}(\varphi,\sigma)$. The transformation
$(\Phi(\mathbf{x}),\Sigma(\mathbf{x}))$ that maps the state space variables
into the reduced space variables is called _inverse transformation_. The
latter is associated to the eigenfunctions $\phi(\mathbf{x})$ of the Koopman
operator defined up to some constant multiplier or shift as discussed above.

It was shown in [31, 32] that the forward transformation (the parametrization)
can be approximated around the limit cycle by Fourier-Taylor series

$K(\varphi,\sigma)=\sum^{\infty}_{n=0}K_{n}(\varphi)\sigma^{n},$ (14)

where $K_{n}(\varphi)$ are $2\pi$ periodic functions in $\varphi$, which are
approximated using Fourier series as follows

$K_{n}(\varphi)=\sum_{k=-\infty}^{\infty}q_{n,k}e^{\mathrm{i}k\varphi}.$ (15)

Using the invariance equation Eq. (9), it is easy to show that
$K_{0}(\varphi)$ is a parameterization of the limit cycle $\Gamma$.

Following [31], we insert Eq. (14) into Eq. (9) and expand the expression
using Taylor series around the limit cycle, i.e. $\sigma=0$. Then, by
collecting coefficients of $\sigma$ on both sides, we obtain a homogeneous
linear differential equation for $K_{1}(\varphi)$:

$\omega\frac{d}{d\varphi}K_{1}(\varphi)+\lambda
K_{1}(\varphi)=D\mathbf{F}(K_{0}(\varphi))\,K_{1}(\varphi),$ (16)

with $\omega=2\pi/T$. Thereby, $K_{1}$ is an eigenfunction of the linear
differential operator

$\mathcal{L}\triangleq\omega\frac{d}{d\varphi}-D\mathbf{F}(K_{0}(\varphi))$

with eigenvalue $-\lambda$. As it is shown in [30] and [31], we have that
$K_{1}(\varphi)$ is a periodic eigenfunction of $\mathcal{L}$ with eigenvalue
$-\lambda$ if and only if $K_{1}(0)$ is an eigenvector of the monodromy matrix
$\mathbf{M}_{\varphi=2\pi}$ with eigenvalue $e^{\lambda T}$. The matrix
$\mathbf{M}_{\varphi=2\pi}$ is a fundamental solution of the first variational
equation

$\frac{d}{d\varphi}\mathbf{M}_{\varphi}=\frac{T}{2\pi}D\mathbf{F}(K_{0}(\varphi))\mathbf{M}_{\varphi},$
(17)

evaluated at $\varphi=2\pi$, with initial condition
$\mathbf{M}_{\varphi=0}=\mathbf{I}$. Once $K_{1}(0)$ is chosen as an
eigenvector, $K_{1}(\varphi)$ is determined for all $\varphi\in[0,2\pi)$ by
solving Eq. (16), that is, $K_{1}(\varphi)=e^{-\lambda
T/(2\pi)\varphi}M_{\varphi}K_{1}(0)$.

By repeating the procedure of inserting Eq. (14) into Eq. (9) and collecting
terms with $\sigma^{n}$ one can derive non-homogeneous equations for
$K_{n}(\varphi)$, $n\geq 2$ which depend on $K_{j}(\varphi),j<n$, thereby
finding the higher terms of Eq. (14) as shown in [31]. Moreover, as it is
shown in [57, 34], using the expression for the parametrization
$K(\varphi,\sigma)$ (the forward transformation function) one can find the
expression for the inverse transformation functions $\Phi(\mathbf{x})$ and
$\Sigma(\mathbf{x})$.

###  2.2 Phase-amplitude reduction for an uncoupled oscillatory system

In this section we describe how to find the transformation that performs the
phase-amplitude reduction of an uncoupled system. For this we assume that we
have a good approximation of the vector field (VF) of the uncoupled system in
a certain (not necessarily small) region around the limit cycle in the
observable space $(\theta,r)$, e.g. the radius and the angle of a signal
measured from this uncoupled system, given by:

$\displaystyle\dot{\theta}$ $\displaystyle=F^{(\theta)}(\theta,r),$
$\displaystyle\dot{r}$ $\displaystyle=F^{(r)}(\theta,r),$

where $F^{(\theta)}(\theta,r)$ and $F^{(r)}(\theta,r)$ are the two components
of the vector field $\mathbf{F}(\theta,r)$ expressed in polar coordinates. How
a good approximation of the VF can be calculated will be elaborated in the
next section. Recall that knowing the vector field in polar coordinates is
equivalent to know it in the cartesian variables $(x,y)$ via standard
transformations.

Furthermore, we assume that the approximated system has a limit cycle, which
could be represented as $r=\tilde{\gamma}^{(r)}(\theta)$, and be approximated
using Fourier series

$\tilde{\gamma}^{(r)}(\theta)=\sum_{k=-N_{g}}^{N_{g}}g_{k}e^{\mathrm{i}k\theta}.$
(18)

The observable variables are the radius $r$ and angle $\theta$ of the measured
signal. Given the vector field $\mathbf{F}$ describing the dynamics, then the
state space (in polar coordinates) is equivalent to the observable space. To
distinguish the variables further in the text, we refer to $(\theta,r)$ as
observable variables and $(\varphi,\sigma)$ as reduced variables (see Fig. 2).

We search for forward and inverse transformation functions in Fourier-Taylor
series form as given in Eq. (14) and Eq. (15). Namely, for the observable
angle and radius we have:

$\theta=K^{(\theta)}(\varphi,\sigma)=\varphi+\sum_{n=0}^{\infty}\sum_{k=-\infty}^{\infty}q^{(\theta)}_{n,k}\sigma^{n}e^{\mathrm{i}k\varphi},$
(19)
$r=K^{(r)}(\varphi,\sigma)=\sum_{n=0}^{\infty}\sum_{k=-\infty}^{\infty}q^{(r)}_{n,k}\sigma^{n}e^{\mathrm{i}k\varphi},$
(20)

Note that for $n=0$, i.e., on the limit cycle, we have:

$\displaystyle\theta$ $\displaystyle=$
$\displaystyle\varphi+\sum_{k=-\infty}^{\infty}q^{(\theta)}_{0,k}e^{\mathrm{i}k\varphi}=\gamma^{(\theta)}(\varphi),$
$\displaystyle r$ $\displaystyle=$ $\displaystyle\gamma^{(r)}(\varphi),$

where the form of the phase transformation coincides with the forward
transformation of the phase reduction given in [16, 23].

For the reduced angle and amplitude we have:

$\varphi=\Phi(\theta,r)=\theta+\sum_{n=0}^{\infty}\sum_{k=-\infty}^{\infty}q^{(\varphi)}_{n,k}r^{n}e^{\mathrm{i}k\theta},$
(21)

and

$\sigma=\Sigma(\theta,r)=\sum_{n=0}^{\infty}\sum_{k=-\infty}^{\infty}q^{(\sigma)}_{n,k}r^{n}e^{\mathrm{i}k\theta},$
(22)

where $q^{(\theta)}_{n,k}$, $q^{(r)}_{n,k}$, $q^{(\varphi)}_{n,k}$,
$q^{(\sigma)}_{n,k}$ are the coefficients we need to find. Here again, if we
set $r=\tilde{\gamma}^{(r)}(\theta)$ from Eq. (18) into Eq. (21) the Fourier-
Taylor series of the phase transformation turns into Fourier series, which
coincides with the inverse transformation of the phase reduction given in
[23].

####  2.2.1 Fourier-Laplace averaging integrals

In the previous section 2.1.1 we showed that the transformation we search for
is the eigenfunction of the Koopman operator. Recall that we denote the limit
cycle of an uncoupled system as $\Gamma$ and the trajectory of system Eq. (1)
starting at $\mathbf{x}_{0}$ as $s(t,\mathbf{x}_{0})$. In [39] it is
demonstrated that an eigenfunction of the Koopman operator that corresponds to
the eigenvalue $\mu$ can be found using the time averaging integral over an
observable function $g(\mathbf{x})$ along the trajectory provided that the
averaging is not equal to zero everywhere and the observable function
$g(\mathbf{x})$ fulfills specific properties, i.e.

$g^{*}_{\mu}(\mathbf{x}_{0})=\lim\limits_{T\rightarrow\infty}\frac{1}{T}\int_{0}^{T}g(s(\tau,\mathbf{x}_{0}))e^{-\mu\tau}d\tau.$
(23)

Namely, if the eigenvalue $\mu$ is real, the integral (23) is called Laplace
average and the observable function should vanish on the limit cycle
$g(\mathbf{x})=0,\forall\mathbf{x}\in\Gamma$. Whereas, if the eigenvalue is
imaginary, $\mu=i\omega$, $\omega\in\mathbb{R}$, then the integral is called
Fourier average and it is equal to the eigenfunction that corresponds to
$i\omega$, i.e. $\phi_{i\omega}$, for almost all choices of the observable
function $g(\mathbf{x})$ (for details see [38, 39]).

The integral Eq. (23) can be used to find the inverse phase-amplitude
reduction transformations $\Phi(\mathbf{x})$ and $\Sigma(\mathbf{x})$. To
confirm this, we apply the Laplace averaging integral to an observable
function defined as
$g(s(t,\mathbf{x_{0}}))=r(t)-\tilde{\gamma}^{(r)}(\theta(t))$, which vanishes
on the limit cycle.

For an uncoupled planar system the solution of Eq. (7) in the reduced space is

$\displaystyle\varphi(t)$ $\displaystyle=$ $\displaystyle\omega
t+\varphi_{0},$ (24) $\displaystyle\sigma(t)$ $\displaystyle=$
$\displaystyle\sigma_{0}e^{\lambda t},$ (25)

where $\varphi_{0}$ and $\sigma_{0}$ are the initial conditions. Using this
solution and the approximation given in Eq. (20) we find the Laplace averaging
integral for the initial condition $\mathbf{x}_{0}$ in the observable space:

$g^{*}_{\lambda}(\mathbf{x}_{0})=\lim\limits_{T\rightarrow\infty}\frac{1}{T}\int_{0}^{T}\left[r(\tau)-\tilde{\gamma}^{(r)}(\theta(\tau))\right]e^{-\lambda\tau}d\tau\\\
=\lim\limits_{T\rightarrow\infty}\frac{1}{T}\int_{0}^{T}\left[K^{(r)}(\varphi(\tau),\sigma(\tau))-K^{(r)}(\varphi(\tau),0)\right]e^{-\lambda\tau}d\tau\\\
=\lim\limits_{T\rightarrow\infty}\frac{1}{T}\int_{0}^{T}\sum_{n=1}^{\infty}\sum_{k=-\infty}^{\infty}q^{(r)}_{n,k}\sigma_{0}^{n}e^{\mathrm{i}k\varphi_{0}}e^{(n-1)\lambda\tau+\mathrm{i}k\omega\tau}d\tau.$

Since we have a hyperbolic limit cycle the eigenvalue is negative, i.e.
$\lambda<0$. Therefore, taking the limit $T\rightarrow\infty$ all terms tend
to 0 except the term for $n=1$. Hence,

$g^{*}_{\lambda}(\mathbf{x}_{0})=\sigma_{0}\sum_{k=-\infty}^{\infty}q^{(r)}_{1,k}e^{\mathrm{i}k\varphi_{0}}\lim\limits_{T\rightarrow\infty}\frac{1}{T}\int_{0}^{T}e^{\mathrm{i}k\omega\tau}d\tau=\\\
=\sigma_{0}\left(q^{(r)}_{1,0}+\lim\limits_{T\rightarrow\infty}\frac{1}{T}\sum_{\begin{subarray}{c}k\in\mathbb{Z}\\\
k\neq
0\end{subarray}}q^{(r)}_{1,k}\frac{e^{\mathrm{i}k\varphi_{0}}}{\mathrm{i}k\omega}\left(e^{\mathrm{i}k\omega
T}-1\right)\right)=\sigma_{0}q^{(r)}_{1,0}.$

Therefore

$g^{*}_{\lambda}(\mathbf{x}_{0})=q^{(r)}_{1,0}\sigma_{0}=q^{(r)}_{1,0}\Sigma(\mathbf{x}_{0}),$
(26)

which indicates that the Laplace averaging integral gives us the inverse
transformation function $\Sigma$ scaled by a constant factor $q^{(r)}_{1,0}$.

In the same way we can show that the angle of a Fourier averaging integral
along the trajectory over an arbitrary observable function (for example, the
$x$-component of the state space or the measured signal) results in the
inverse transformation function for $\varphi$:

$g^{*}_{\mathrm{i}\omega}(\mathbf{x}_{0})=\lim\limits_{T\rightarrow\infty}\frac{1}{T}\int_{0}^{T}x(\tau)e^{-\mathrm{i}\omega\tau}d\tau\\\
=\lim\limits_{T\rightarrow\infty}\frac{1}{T}\int_{0}^{T}K^{(x)}(\varphi(\tau),\sigma(\tau))e^{-\mathrm{i}\omega\tau}d\tau\\\
=\lim\limits_{T\rightarrow\infty}\frac{1}{T}\int_{0}^{T}\sum_{n=0}^{\infty}\sum_{k=-\infty}^{\infty}q^{(x)}_{n,k}\sigma_{0}^{n}e^{\mathrm{i}k\varphi_{0}}e^{n\lambda\tau+\mathrm{i}(k-1)\omega\tau}d\tau\\\
=q^{(x)}_{0,1}e^{\mathrm{i}\varphi_{0}}+\lim\limits_{T\rightarrow\infty}\frac{1}{T}\sum_{\begin{subarray}{c}k\in\mathbb{Z}\\\
k\neq
1\end{subarray}}q_{0,k}^{(x)}e^{\mathrm{i}k\varphi_{0}}\frac{e^{\mathrm{i}(k-1)\omega
T}-1}{\mathrm{i}(k-1)\omega}\\\
+\lim\limits_{T\rightarrow\infty}\frac{1}{T}\sum_{n=1}^{\infty}\sum_{k\in\mathbb{Z}}q^{(x)}_{n,k}\sigma_{0}^{n}e^{\mathrm{i}k\varphi_{0}}\frac{e^{(n\lambda+\mathrm{i}(k-1)\omega)T}-1}{(n\lambda+\mathrm{i}(k-1)\omega)}\\\
=q^{(x)}_{0,1}e^{\mathrm{i}\varphi_{0}}.$ (27)

Thus,

$\angle g^{*}_{\mathrm{i}\omega}(\mathbf{x}_{0})=\angle
q^{(x)}_{0,1}e^{\mathrm{i}\varphi_{0}}=\varphi_{0}=\Phi(\mathbf{x}_{0}).$ (28)

####  2.2.2 Numerical calculation of the inverse and forward transformations

Numerical calculation of the Laplace averaging integral Eq. (23) is
challenging. First of all, we need the values of $\lambda$ and $\omega$ of the
reduced system Eq. (7). Second, an optimal time range $T$ is needed to
calculate the time averaging integral Eq. (23) (see. [39]). In this section,
we provide the details of these computations.

Knowing the VF $\mathbf{F}(\mathbf{x})$ describing the system dynamics, we
estimate the value of $\omega$ by finding the period of the limit cycle, and
the value of $\lambda$ by finding the eigenvalues of the monodromy matrix (Eq.
(17)). Namely, we simulate the system for a long enough period of time until
the radius relaxes to the limit cycle $\tilde{\gamma}^{(r)}(\theta)$. The
latter can be found by finding the coefficients $g_{k}$ of the Fourier series
in Eq. (18) by fitting the time courses $\theta(t)$ and $r(t)$ on the limit
cycle by means of the least squares error method (see. A).

Then, we find the derivatives of $\mathbf{F}$ with respect to $r$ and $\theta$
along the limit cycle and solve the linear variational equation (17) for one
period using the initial condition $\mathbf{M}_{\varphi=0}=\mathbf{I}$. The
eigenvalues of the resulting monodromy matrix $\mathbf{M}_{\varphi=2\pi}$ are
equal to $\hat{\lambda}_{k}=e^{\lambda_{k}T}=e^{\lambda_{k}2\pi/\omega}$. The
largest negative $\lambda_{k}$ is the searched value of $\lambda$.

The calculation of the Laplace averaging integral is restricted by numerical
precision [39]. Let us consider the time course of an observable

$\rho(t)=r(t)-\tilde{\gamma}^{(r)}(\theta(t)).$ (29)

The time course of $\rho(t)$ is shown in the upper left panel in Fig. 3. Due
to the limitation of the precision of the numerical integration, the
multiplication of two competing functions - exponentially growing $e^{-\lambda
t}$ and $\rho(t)$ that exponentially approaches the limit cycle - leads to
numerical errors of the averaging integral. This can be demonstrated by
plotting $|\rho(t)|$ in logarithmic scale (blue curve in right panel of Fig.
3), where numerical instabilities occur after reaching the numerical precision
threshold (horizontal dashed black line in right panel of Fig. 3).

However, one can find the initial values $\sigma_{0}$ and $\varphi_{0}$
without calculating the integral, but using a finite time and logarithmic
scale as discussed further.

The observable $\rho(t)$ can be expressed as

$\rho(t)=r(t)-\tilde{\gamma}^{(r)}(\theta(t))=K^{(r)}(\varphi(t),\sigma(t))-K^{(r)}_{0}(\varphi(t))\\\
=\sum_{n=1}^{\infty}\sum_{k=-\infty}^{\infty}q^{(r)}_{n,k}\sigma_{0}^{n}e^{\mathrm{i}k\varphi_{0}}e^{n\lambda
t+\mathrm{i}k\omega t}.$

If we factor out $\sigma_{0}e^{\lambda t}$ we have

$\rho(t)=\sigma_{0}e^{\lambda
t}\sum_{n=1}^{\infty}\sum_{k=-\infty}^{\infty}q^{(r)}_{n,k}\sigma_{0}^{n-1}e^{\mathrm{i}k\varphi_{0}}e^{(n-1)\lambda
t+\mathrm{i}k\omega t}\\\ =\sigma_{0}e^{\lambda t}R(t).$ (30)

Since $\lambda<0$, after a transient time, the terms involving $e^{\lambda t}$
in the expression for $R(t)$ decay down to zero and, therefore, only
oscillatory terms with $n=1$ remain, i.e.

$R(t)\sim\sum_{k=-\infty}^{\infty}q^{(r)}_{1,k}e^{\mathrm{i}k\varphi_{0}}e^{\mathrm{i}k\omega
t}\triangleq\hat{R}(t),$ (31)

the average of which over one period is

$\langle\hat{R}(t)\rangle_{T}=\left\langle\sum_{k=-\infty}^{\infty}q^{(r)}_{1,k}e^{\mathrm{i}k\varphi_{0}}e^{\mathrm{i}k\omega
t}\right\rangle_{T}=q^{(r)}_{1,0}.$

Hence, from Eqs. (30) and (31) we have

$\langle|\rho(t)|e^{-\lambda
t}\rangle_{T}=\langle|\sigma_{0}R(t)|\rangle_{T}\sim|q^{(r)}_{1,0}\sigma_{0}|.$
(32)

On the other hand we have

$\log|\rho(t)|\sim\lambda t+\log\sigma_{0}\hat{R}(t),$ (33)

that is, a line with slope $\lambda$ plus an oscillatory part around
$\log|q^{(r)}_{1,0}\sigma_{0}|$. This is shown in the right panel of Fig. 3,
where $|\rho(t)|$ is plotted in semi-logarithmic scale. Here we see that the
calculation of $|\rho(t)|$ becomes unstable after approximately the time point
$t_{2}$, which is defined using the precision threshold of the numerical
integration (horizontal dashed line). If we consider an interval of several
periods before reaching this threshold ($t\in[t_{1},t_{2}]$, where
$t_{2}=t_{1}+2T$), the time course of $|\rho(t)|$ in logarithmic scale
oscillates around the function $q^{(r)}_{1,0}\sigma_{0}e^{\lambda t}$, which
appears as a line in logarithmic scale (red dotted line).

Therefore, we average the last several periods of $|\rho(t)|e^{-\lambda t}$
before reaching the numerical precision threshold to find
$|q^{(r)}_{1,0}\sigma_{0}|$, as shown in the bottom panel of Fig. 3, where an
interval with the duration of two periods (red line) was taken. Notice that
the time course of $|\rho(t)|e^{-\lambda t}$ after some transient time is in
general oscillatory (see Eq. (31)) as shown in Fig. 13 in the Appendix.
Moreover, if the initial approximation of $\lambda$ using the solution of Eq.
(17) is not precise, one can find the correction to the value of $\lambda$ by
using the average slope of the curve $\log|\rho(t)|-\lambda t$, which is also
shown in Fig. 13 in the Appendix.

As we discussed previously, the reduced amplitude $\sigma$ can be defined up
to a constant multiplier. Thus, we can redefine the amplitude $\sigma_{0}$ as
$q^{(r)}_{1,0}\sigma_{0}$ and thereby fixing the scaling of $\sigma$. The sign
of $\sigma_{0}$ is defined by the sign of $\rho(0)$.

Figure 3: Illustration of how to find the amplitude variable $\sigma_{0}$
corresponding to an initial point $(\theta_{0},r_{0})$ in the observable space
for the canonical model close to the numerical instabilities. Upper left
panel: Time course of the observable $\rho$ introduced in Eq.(30) (solid blue
curve) that vanishes on the limit cycle $\Gamma$. Right panel: The absolute
value of $\rho$ in semi-logarithmic scale, where the expression
$\sigma_{0}e^{\lambda t}$ appears as a line (dotted red line) with constant
part equal to $\sigma_{0}$ (red star, the crossing with the vertical axes).
The horizontal dashed black line is the threshold of numerical precision, and
the right vertical dot-dashed black line marks the time point $t_{2}$ at this
threshold. The left vertical dot-dashed black line at $t_{1}$ is in two-
period-distance to the right one. The solid red line denotes the time course
of $|\rho(t)|,t\in[t_{1},t_{2}]$. Bottom left panel: The time course of
$|\rho(t)|e^{-\lambda t}$ (solid blue curve) approaches the horizontal line
$\langle|\hat{R}(t)\sigma_{0}|\rangle_{T}\sim|q^{(r)}_{1,0}\sigma_{0}|$
(dotted red line, red star) before the numerical instabilities. The solid red
line and the two vertical dot-dashed black lines are the same as in the right
panel. The vertical dotted lines indicate the oscillation periods.

In a similar way we can find the initial value $\varphi_{0}$ by finding the
angle of the average of $x(t)e^{-\mathrm{i}\omega
t}=r(t)\cos(\theta(t))e^{-\mathrm{i}\omega t}$ over several periods of
oscillations after a transient time and before reaching the numerical
precision threshold. This is equivalent to the calculation of the integral in
Eq. (27), where the terms involving $e^{n\lambda t}$ decay to 0 after some
transient time and the terms involving $e^{\mathrm{i}k\omega T}$ vanish due to
time averaging over several periods. We define $\varphi_{0}$ as
$\varphi_{0}=0$ when $\theta_{0}=0$.

Having computed several pairs of initial conditions in the observable space
$\\{\theta_{j,0},r_{j,0}\\}_{j\in\mathbb{N}}$ and the corresponding ones in
the reduced space $\\{\varphi_{j,0},\sigma_{j,0}\\}_{j\in\mathbb{N}}$ we can
then approximate the forward transformation functions in the forms given in
Eqs. (19) and (20) using the regularized least squares method (see A) with the
Fourier-Taylor basis functions

$\Psi^{(\chi)}=\left\\{\sigma^{n}e^{\mathrm{i}k\varphi},n=0,\ldots,N_{n},k=-N_{k},\ldots,N_{k}\right\\},$
(34)

where the superscript $(\chi)$ denotes the dependence on the reduced variables
$\varphi$ or $\sigma$.

In a similar way we can find the inverse transformations by approximating them
as in Eqs. (21) and (22) using pairs of points (initial conditions) in the
reduced space $\\{\varphi_{j,0},\sigma_{j,0}\\}_{j\in\mathbb{N}}$ and the
observable space $\\{\theta_{j,0},r_{j,0}\\}_{j\in\mathbb{N}}$ by means of the
regularized least squares method (see A) with the basis functions

$\Psi^{(x)}=\left\\{r^{n}e^{\mathrm{i}k\theta},n=0,\ldots,N_{n},k=-N_{k},\ldots,N_{k}\right\\},$
(35)

where again the superscript $(x)$ denotes the dependence on the observable
variables $\theta$ and $r$.

Finally, the whole pipeline of finding the transformation functions is the
following:

  1. 1.

Knowing the VF of a system or a good approximation of the VF on a given
region, numerically integrate it and find the limit cycle $\Gamma$, the
frequency $\omega$, and the approximation of the limit cycle Eq. (18) using
LSQ method (see. A).

  2. 2.

Using the limit cycle, solve the variational equation Eq. (17) and estimate
$\lambda$ by finding the eigenvalue of the monodromy matrix
$\mathbf{M}_{\varphi=2\pi}$.

  3. 3.

Numerically integrate the VF with different initial conditions on a grid
$\\{\theta_{j,0},r_{j,0}\\}_{j\in\mathbb{N}}$ for several periods and obtain
the time courses for $r(t)$ and $\theta(t)$.

  4. 4.

Compute $\rho(t)$ and find the correction to the value of $\lambda$ by finding
the average slope of $\log|\rho(t)|-\lambda t$ in the last two periods of the
simulation (or two periods before reaching the numerical precision threshold).
See C for details.

  5. 5.

Find $|\sigma_{0}|$ for a given $(\theta_{0},r_{0})$ by taking the average of
$|\rho(t)|e^{-\lambda t}$ over the last two periods (or two periods before
reaching the numerical precision threshold) and setting
$\mathrm{sign}(\sigma_{0})=\textrm{sign}(\rho(0))$.

  6. 6.

Find $\varphi_{0}$ for a given $(\theta_{0},r_{0})$ as the angle of the time
average of $r(t)\cos(\theta(t))$ over the last two periods of the time course
(or two periods before reaching the numerical precision threshold).

  7. 7.

Use the values of the initial conditions in the observable space
$\\{\theta_{j,0},r_{j,0}\\}_{j\in\mathbb{N}}$ and the reduced space
$\\{\varphi_{j,0},\sigma_{j,0}\\}_{j\in\mathbb{N}}$ to approximate the forward
transformations given in Eqs. (20) and (19) by means of the regularized least
squares method (see A) with the basis functions in Eq. (34).

  8. 8.

Similarly, approximate the inverse transformations in Eqs. (22) and (21) by
means of the regularized least squares method (see A) with the basis functions
in Eq. (35).

###  2.3 Amplitude and phase reduction for coupled oscillatory systems

In the previous section, we showed how to find the transformation functions
under the assumption that the VF of a system or a good approximation of it is
known. In this section, we present a method of finding the VF and the phase-
amplitude reduction for coupled systems using the signals measured from these
systems.

####  2.3.1 Approximation of the VF for coupled oscillatory systems

Let us assume that we have a set of $M$ trials of measured signals from $N$
different sources with high time resolution,
$\\{s_{i}(t)\\}_{p},\,i=1,\ldots,N,\,p=1,\ldots M$. Then, using methods to
extract the amplitude and angle of the signals, for example the Hilbert
transformation [58, 59], we obtain a set of trials with radius and angle time
courses $\\{\theta_{i}(t),r_{i}(t)\\}_{p}$, which we refer to as observable
amplitude and phase. Using finite differences we obtain the time derivatives
of the observable radius and phase
$\\{\dot{\theta}_{i}(t),\dot{r}_{i}(t)\\}_{p}$.

The dynamics of a region $i$ can be represented as

$\displaystyle\dot{x}_{i}=\mathbf{F}^{(x)}_{i}(\theta_{i},r_{i},\ldots,v_{j}(\theta_{j},r_{j}),\ldots),$
(36)

where $x$ is either $\theta$ or $r$, and $v_{j}$ is the input from other
regions coupled to region $i$, e.g. $v_{j}=r_{j}\cos(\theta_{j})$. As we have
shown in the previous sections, in order to find the transformation functions
we need an expression for the VF of the uncoupled system. The function
$\mathbf{F}_{i}$, in general, is non-linear with respect to its variables,
thus the dynamics of the uncoupled system can not be separated from the
dynamics induced by the coupling. However, if we approximate $\mathbf{F}_{i}$
using a Fourier-Taylor series, then all the terms in the series that depend
only on $\theta_{i}$ and $r_{i}$ correspond to the uncoupled dynamics, and the
rest to the coupling. Moreover, we assume that inputs from other regions do
not interact with each other and we can approximate the coupling part by only
considering pairwise interactions.

Thus, the approximated VF of region $i$ consists of its uncoupled dynamics and
the sum of pairwise couplings with other regions

$\dot{x}_{i}\approx
F^{(x)}_{i}(\theta_{i},r_{i})+\sum_{j}G^{(x)}_{i,j}(\theta_{i},r_{i},\theta_{j},r_{j}),$
(37)

where $x$ is again either $\theta$ or $r$. In Eq. (37) we refer to the
function representing the uncoupled dynamics $F^{(x)}_{i}(\theta_{i},r_{i})$
as uncoupled VF, which we approximate as

$F^{(x)}_{i}(\theta_{i},r_{i})=\sum_{n=0}^{N_{n}}\sum_{k=-N_{k}}^{N_{k}}q^{(x_{i})}_{n,k}r_{i}^{n}e^{\mathrm{i}k\theta_{i}},$
(38)

where $q^{(x_{i})}_{n,k}$ are the coefficients of the Fourier-Taylor series
for the oscillator $i$.

The pairwise coupling function
$G^{(x)}_{i,j}(\theta_{i},r_{i},\theta_{j},r_{j})$, we refer to as coupling
part of the VF, depends on the variables $\theta_{j}$ and $r_{j}$ of other
regions. We approximate it by a Fourier-Taylor series with the variables
$(\theta_{i},r_{i},\theta_{j},r_{j})$ as follows

$G^{(x)}_{i,j}(\theta_{i},r_{i},\theta_{j},r_{j})=\sum_{m=0}^{N_{m}}\sum_{l=0}^{m-l}\sum_{k_{i}=-N_{k_{i}}}^{N_{k_{i}}}\\\
\sum_{\begin{subarray}{c}k_{j}=-N_{k_{j}}\\\ k_{j}\neq
0\end{subarray}}^{N_{k_{j}}}c^{(x_{i,j})}_{m,l,k_{i},k_{j}}r_{i}^{l}r_{j}^{m-l}e^{\mathrm{i}(k_{i}\theta_{i}+k_{j}\theta_{j})},$
(39)

where $c^{(x_{i,j})}_{m,l,k_{i},k_{j}}$ are the coefficients which need to be
found, $x$ is either $\theta$ or $r$, and the subindices $i,j$ of the upper
indices $(x_{i,j})$ indicate the relation to the coupling $G^{(x)}_{i,j}$,
which is approximated up to the sum of the powers of $r_{i}$ and $r_{j}$ equal
to $N_{m}$. Note, the sum in Eq. (39) does not include terms with $k_{j}=0$.
The terms with $r_{j}^{0}$ and $k_{j}=0$ correspond to the uncoupled VF (Eq.
(38)), since here only terms with $r_{i}$ and $\theta_{i}$ remain. The terms
with $r_{j}^{m}$, $m\neq 0$ and $k_{j}=0$ are excluded to provide the
dependence on $v_{j}$, since the variables $r_{j}$ and $\theta_{j}$ do not
appear in the coupling part $G^{(x)}_{i,j}$ separately, but only in
combination via $v_{j}$.

In practice, the approximations in Eq. (37)-(38)-(39) are obtained by means of
the least squares method (see A) with a combination of the following basis
functions

$\Psi^{(x)}_{i}=\\{r_{i}^{n}e^{\mathrm{i}k\theta_{i}},n=0,\ldots,N_{n}\text{
and }k=-N_{k},\ldots,N_{k}\\},$ (40)

and

$\displaystyle\Psi^{(x)}_{ij}=\left\\{r_{i}^{m_{i}}r_{j}^{m_{j}}e^{\mathrm{i}(k_{i}\theta_{i}+k_{j}\theta_{j})}\right.,$
(41) $\displaystyle m_{i},m_{j}=0,\ldots,N_{m},\;m_{i}+m_{j}\leq N_{m},$
$\displaystyle k_{i}=-N_{k_{i}},\ldots,N_{k_{i}},$
$\displaystyle\left.k_{j}=-N_{k_{j}},\ldots,N_{k_{j}},\;k_{j}\neq 0\right\\},$

for every input region $j$. The superscript $(x)$ in Eqs. (40) and (41)
denotes the dependence on the state variables $\theta$ and $r$, and the lower
indices $i$ or $j$ denote of which region. Note, that we have different ranges
for the indices $n$ and $k$ in Eq. (40) and for the indices $m$, $k_{i}$ and
$k_{j}$ in Eq. (41).

Using the above equations, we can find an approximation of the uncoupled VF
$F^{(x)}_{i}$ and the coupling part of the VF $G^{(x)}_{i,j}$ for every region
and the incoming input to this region. The obtained functions give a good
approximation of the uncoupled VF and coupling in the domain that is covered
by the measured data. Furthermore, using the regularized least squares method
(see A) we escape over-fitting and provide a good approximation even in the
surrounding of the domain covered by the measured data, as it will be shown in
the results part.

####  2.3.2 Finding the coupling in the reduced space

The equations of the dynamics of a system with coupling in the reduced space
read as:

$\left\\{\begin{array}[]{cc}\dot{\varphi}_{i}=\omega_{i}+\sum_{j}g^{(\varphi)}_{i,j}(\varphi_{i},\sigma_{i},\varphi_{j},\sigma_{j}),\\\
\dot{\sigma}_{i}=\lambda_{i}\sigma_{i}+\sum_{j}g^{(\sigma)}_{i,j}(\varphi_{i},\sigma_{i},\varphi_{j},\sigma_{j}),\end{array}\right.$
(42)

where $g^{(\varphi)}_{i,j}$ and $g^{(\sigma)}_{i,j}$ are the pairwise coupling
functions between the oscillator $i$ and $j$.

Using the invariance equations Eqs. (12)-(13), and taking into account the
approximation of the VF Eqs. (LABEL:eq:VFappx)-(38)-(39) we obtain the
equations for $g^{(\varphi)}_{i,j}$ and $g^{(\sigma)}_{i,j}$:

$\displaystyle\sum_{j}g^{(\varphi)}_{i,j}=\frac{\partial\Phi_{i}}{\partial\theta_{i}}\sum_{j}G^{(\theta)}_{i,j}+\frac{\partial\Phi_{i}}{\partial
r_{i}}\sum_{j}G^{(r)}_{i,j},$ (43)
$\displaystyle\sum_{j}g^{(\sigma)}_{i,j}=\frac{\partial\Sigma_{i}}{\partial\theta_{i}}\sum_{j}G^{(\theta)}_{i,j}+\frac{\partial\Sigma_{i}}{\partial
r_{i}}\sum_{j}G^{(r)}_{i,j},$ (44)

where $\Phi_{i}=\Phi_{i}(\theta_{i},r_{i})$ and
$\Sigma_{i}=\Sigma_{i}(\theta_{i},r_{i})$ are the inverse transformation
functions approximated using the method presented in Section 2.2.2. Assuming
independence of the coupling terms between different oscillators, we express
the coupling in the reduced space separately for every pairwise coupling
$\\{i,j\\}$:

$\displaystyle
g^{(\varphi)}_{i,j}=\frac{\partial\Phi_{i}}{\partial\theta_{i}}G^{(\theta)}_{i,j}+\frac{\partial\Phi_{i}}{\partial
r_{i}}G^{(r)}_{i,j},$ (45) $\displaystyle
g^{(\sigma)}_{i,j}=\frac{\partial\Sigma_{i}}{\partial\theta_{i}}G^{(\theta)}_{i,j}+\frac{\partial\Sigma_{i}}{\partial
r_{i}}G^{(r)}_{i,j}.$ (46)

Note that the terms on the right hand side of Eqs. (45) and (46) depend on the
observable variables $(\theta,r)$. In the practical implementation, we define
a four dimensional grid $\\{\theta_{i},r_{i},\theta_{j},r_{j}\\}$ for every
$j\in[1,N]$ and find the values of $g^{(\varphi)}_{i,j}$ and
$g^{(\sigma)}_{i,j}$ with respect to this grid in the observable space. Then,
using the inverse transformations we find the corresponding four dimensional
grid in the reduced space
$\\{\varphi_{i},\sigma_{i},\varphi_{j},\sigma_{j}\\}$ and approximate the
coupling using Fourier-Taylor series as follows

$g^{(\chi)}_{i,j}(\varphi_{i},\sigma_{i},\varphi_{j},\sigma_{j})\approx\sum_{n_{i}=0}^{N_{n_{i}}}\sum_{n_{j}=0}^{N_{n_{j}}}\sigma_{i}^{n_{i}}\sigma_{j}^{n_{j}}\\\
\sum_{k_{i}=-N_{k_{i}}}^{N_{k_{i}}}\sum_{\begin{subarray}{c}k_{j}=-N_{k_{j}}\\\
k_{j}\neq
0\end{subarray}}^{N_{k_{j}}}p^{(\chi_{i,j})}_{n_{i},n_{j},k_{i},k_{j}}e^{\mathrm{i}(k_{i}\varphi_{i}+k_{j}\varphi_{j})},$
(47)

where $\chi$ is either $\varphi$ or $\sigma$, and again the subindices $i,j$
of the upper indices $(\chi_{i,j})$ denote the reference to the coupling
$g^{(\chi)}_{i,j}$.

Then, the coefficients $p^{(\chi_{i,j})}_{n_{i},n_{j},k_{i},k_{j}}$ together
with $\lambda_{i}$ and $\omega_{i}$ give us the whole set of parameters of the
amplitude and phase equations Eq. (42) in the reduced space.

The approximation Eq. (47) can be written in a form suitable for the
regularized least squares method (see A) using the basis functions

$\displaystyle\Psi^{(\chi)}_{{i,j}}=\left\\{\sigma_{i}^{n_{i}}\sigma_{j}^{n_{j}}e^{\mathrm{i}(k_{i}\varphi_{i}+k_{j}\varphi_{j})}\right.,$
(48) $\displaystyle
n_{i},n_{j}=0,\ldots,N_{m},k_{i}=-N_{k_{i}},\ldots,N_{k_{i}},$
$\displaystyle\left.k_{j}=-N_{k_{j}},\ldots,N_{k_{j}},\;k_{j}\neq 0\right\\}.$

##  3 Results

###  3.1 Description of the models tested

To test our method we used four models: (i) the radial isochron clock, (ii)
the canonical model for an oscillator, also known as Stuart-Landau oscillator,
(iii) the van der Pol oscillator and (iv) the Wilson-Cowan model (see B). The
first two models have circular limit cycles and known transformation functions
(see Eqs. (56) and (60), and Supplementary Materials for their derivations),
whereas the last two models have non-circular limit cycles and there are no
known analytical expressions for their transformation functions. Moreover, the
phase variable of the uncoupled radial isochron clock does not depend on the
amplitude (see Eq. (54)) and the input in the Wilson-Cowan model is non-linear
(see Eq. (63)). The properties of the models are summarized in Table 1.

Model name |  analytical transformation |  phase depends on amplitude |  circular limit cycle |  linear input
---|---|---|---|---
Radial isochron clock |  known |  no |  yes |  yes
Canonical model |  known |  yes |  yes |  yes
van der Pol oscillator |  unknown |  yes |  no |  yes
Wilson-Cowan |  unknown |  yes |  no |  no
Table 1: Table of the properties of the models studied.

We simulated these models as pairs of two coupled oscillators of the same kind
but with different parameter values (see B for the parameter values). We
generated 100 trials for every pair of oscillators with different initial
conditions. The initial conditions for the simulations were chosen randomly to
cover the region around the limit cycles. The trials lasted several periods of
the oscillations. In all simulations the first oscillator was driving the
second one.

###  3.2 Results of the VF reconstruction

In this Section, we present results of the reconstruction of the VF using
synthetically simulated signals (trials) from the different models presented
above.

From the trials we first find time courses of the radius $r(t)$ and the angle
$\theta(t)$. Then, we compute a numerical approximation of the derivatives
$\dot{r}(t)$ and $\dot{\theta}(t)$. Notice that, in practice, the signals
obtained from an experiment will be first filtered around a given frequency.
Thus, the time courses of the radius and angle, and their derivatives will be
smooth.

Using the approximated derivatives and the time courses of the variables we
approximate the VF with the regularized least squares method (see A) using the
basis functions in Eq. (40) and Eq. (41). This allows us to represent the VF
as a sum of the uncoupled VF and the pairwise coupling from other regions.

The results of the reconstruction of the uncoupled VF for all four different
models are given in the Supplementary Materials. Examples of the
reconstruction of the VF for two coupled radial isochron clocks and two
coupled van der Pol oscillators are given in Figs. 4 and 5, respectively. In
both figures, as goodness of approximation, the maximal deviation of the
reconstructed uncoupled VF from the theoretical one among all coordinates of
the data points is given on the top of every panel. Overall, one can see that
the approximations are close to the theoretical VFs.

Figure 4: An example of the uncoupled VF reconstruction for two non-
identically coupled radial isochron clocks. The left panels correspond to the
first, and the right ones to the second oscillator. The first oscillator
drives the second one. Red dots are the data points (trials). Approximated
uncoupled VFs are depicted as black mesh surfaces. The corresponding
theoretical VFs are plotted as green mesh surfaces. The numbers above each
panel correspond to the maximal deviation of the green and black surfaces at
the data points.  Figure 5: An example of the uncoupled VF reconstruction for
two non-identically coupled van der Pol oscillators. The notations are the
same as in Fig. 4.

Observe that in Fig. 4 the trials (data points in red) for the uncoupled
oscillator (left panels) gather around the surface of the uncoupled VF,
whereas for the coupled oscillator (right panels), the trials strongly deviate
due to the input received from the first oscillator. However, the
approximation of the VF manages to separate the uncoupled VF from the coupling
with a good accuracy, as it can be seen from the maximal deviation of the
theoretical (green mesh) and approximated (black mesh) surfaces at the
coordinates of data points, which is given above each panel.

Since the trials do not uniformly cover the space $(\theta,r)$, some areas on
the plane $(\theta,r)$ have a higher density of data points than others. This
uneven distribution leads to a slight deviation of the values of the
approximated VF from the theoretical one in the areas with less density or
absence of data points. Notably, due to the regularization (ridge regression)
used in the least squares method, the approximation does not strongly deviate
from the theoretical one even outside the region covered by the data points,
allowing for slight deviations in these uncovered areas to still maintain a
close fit to the theoretical model.

Even for a highly nonlinear VF function, such as for the van der Pol
oscillator (see Fig. 5), the VF can be well approximated outside of the
regions covered by the data points.

###  3.3 Reconstruction of the transformations

Having found the uncoupled VF of the oscillatory systems we use the
approximation method presented in Section 2.2.2 to find the individual inverse
and forward transformations of the oscillators. Namely, we build a grid of
initial conditions in the observable space. For every initial condition
$\\{\theta_{0,k},r_{0,k}\\}_{k\in\mathbb{N}}$ in the observable space we
simulate the system for several periods of oscillations and find the
corresponding initial conditions in the reduced space
$\\{\varphi_{0,k},\sigma_{0,k}\\}_{k\in\mathbb{N}}$ using the linear
approximation method discussed in the Methods Section. Then, having a set of
pairs of initial conditions in the observable space and in the reduced space,
the transformation functions were approximated using the regularized least
squares method (see A).

In Fig. 6 we show the comparison of the reconstructed inverse transformations
and the theoretical ones for the radial isochron clock and the canonical model
(see B.1 and B.2). Since for both models the theoretical expression of the
inverse transformation $\sigma=\Sigma(\theta,r)$ does not depend on the phase
$\theta$ (see Eqs. (56) and (60)), the comparison is performed only on the
observable radius $r$ and the reduced amplitude $\sigma$. The blue dots
correspond to the approximated initial values of the amplitude in the reduced
space $\sigma_{0}$.

In Fig. 7 the approximated inverse phase transformation is compared with the
theoretical one for the radial isochron clock (the left panel) and the
canonical model (the right panel). For the radial isochron clock the phase
transformation is $\varphi=\theta$ (Eq. (56)) and does not depend on the
radius (see the left panel in Fig. 7). For the canonical model the inverse
transformation depends logarithmically on the radius (Eq. (60)). In the right
panel of Fig. 7 the approximated initial values of $\varphi$ are compared with
the theoretical ones for different initial values of radius $r$ and different
fixed values of $\theta$. As one can see the approximations of the initial
amplitude and phase, and therefore the approximated amplitude and phase
inverse transformation functions, coincide with the theoretical values for the
radial isochron clock and the canonical model.

In Fig. 8, we show an example of the reconstructed inverse transformations for
the van der Pol oscillator (see B.3) and the Wilson-Cowan model (see B.4). In
both models, the inverse transformations depend non-linearly on both amplitude
and phase. The results of the reconstruction of the transformation functions
for all oscillator models, together with the derivation of the theoretical
expressions for the transformation functions are given in the Supplementary
Materials.

Figure 6: Comparison of the reconstructed inverse amplitude transformations
and theoretical ones for the radial isochron clock (left panel) and the
canonical model (right panel). Blue dots are the approximated initial values
of the amplitude in the reduced space $\sigma$ for the simulated trials with
the initial values of the radius $r$ in the observable space. Red curves are
the corresponding theoretical values of the inverse amplitude transformations.

Figure 7: Comparison of the reconstructed inverse phase transformations and
the theoretical ones for the radial isochron clock (left panel) and the
canonical model (right panel). Left panel: Blue dots are the approximated
values of the inverse phase transformation (the initial values of the phase in
the reduced space $\varphi$) for the radial isochron clock with respect to
different initial phases $\theta$ and with a fixed initial value of the radius
$r=1.2$ in the observable space. Red curve represents the corresponding
theoretical values of the inverse phase transformation for the radial isochron
clock. Right panel: Blue dots are the approximated values of the inverse phase
transformations (the initial values of the phase in the reduced space
$\varphi$) for the canonical model with respect to the initial values of the
radius $r$ for different fixed values of the initial phase $\theta$. Red lines
represent the corresponding theoretical values of the inverse phase
transformation for different fixed values of $\theta$.

Figure 8: An example of the reconstructed inverse transformation functions
$\Sigma(\theta,r)$ and $\Phi(\theta,r)$ for (A) the van der Pol oscillators
and (B) the Wilson-Cowan models. The left column of each panel corresponds to
the amplitude transformation $\sigma_{0}=\Sigma(\theta_{0},r_{0})$ and the
right column corresponds to the phase transformation
$\varphi_{0}=\Phi(\theta_{0},r_{0})$ from the observable to the reduced space.
The upper panels correspond to the first oscillator and the lower panels to
the second oscillator. The red dots are the initial values of the amplitude
and phase in the observable $(\theta_{0},r_{0})$ and reduced ($\varphi_{0}$,
$\sigma_{0}$) spaces. The black surfaces are the approximated inverse
transformations using the regularized least squares method.

###  3.4 Reconstruction of the coupling part

In this Section, we present the reconstruction results of the coupling part of
the VF for the four different models considered. As mentioned in the Methods
Section, we first find the coupling in the observable space using the
approximated coupling part of the VF and the gradients of the inverse
transformations. Therefore, we generate random points in the four dimensional
space $(\theta_{i},r_{i},\theta_{j},r_{j})$ for every pair of coupling $ij$.
The range of random points are chosen to be inside 0.8 of the maximum and
1/0.8 of the minimum radius of the data points. This is done to ensure good
approximation results, since the VF and the transformation functions have a
better approximation close to the limit cycle.

After estimating the approximation of the coupling $G^{(x)}_{i,j}$ in the
observable space as Eq. (39) with the basis functions Eq. (41), we find the
corresponding coordinates of the random points in the reduced space
$(\varphi_{i},\sigma_{i},\varphi_{j},\sigma_{j})$ using the inverse
transformation. We then approximate the coupling $g^{(\chi)}_{i,j}$ in the
reduced space as Eq. (47) with the basis functions Eq. (48). The results of
the approximation are the real-valued Fourier-Taylor coefficients
$p^{(\chi_{i,j})}_{n_{i},n_{j},k_{i},k_{j}}$ that define the coupling in the
amplitude and phase variables (see Eqs. (42)) for every pair of regions $ij$.

We present the approximated Fourier-Taylor coefficients
$p^{(\chi_{i,j})}_{n_{i},n_{j},k_{i},k_{j}}$ as a heat-map (see Fig. 9). The
coefficients are sorted in panels with respect to the amplitude terms of the
corresponding basis functions, i.e. with respect to the combinations of the
powers of $\sigma_{i}$ and $\sigma_{j}$, which is indicated by a label on the
top of every panel. The vertical axis of all panels corresponds to the
components of the basis functions that depend on the own amplitude and phase
(i.e. $\Psi_{i}$), whereas the horizontal axis corresponds to the ones that
depend on the amplitude and phase of an input region ($\Psi_{j}$). Then, an
intersection on the heat map is the value of a Fourier-Taylor coefficient that
corresponds to the multiplication of the vertical and horizontal basis
functions. This way of representation of the coefficients allows us to display
the approximated Fourier-Taylor coefficients as a matrix, which is build as a
product of basis functions as in Eq. (LABEL:eq:basisfuncX).

Figure 9: Explanation of the meaning of the Fourier-Taylor coefficients panels
for a particular example. See text for more details.

Examples of the reconstructed coupling coefficients for the radial isochron
clock and for the canonical model are presented in Figs. 10 and 11. For these
models, an analytic expression of the transformation functions, as well as of
the coupling functions can be derived (see B and Supplementary Materials). The
reconstructed coupling coefficients for these models are non-zero only for the
terms that have a specific combination of the cosine and sine functions (Figs.
10 and 11). For example, for the radial isochron clock the coupling
coefficients are non-zero for the terms involving
$\sin(\varphi_{i})\cos(\varphi_{j})$ (see Fig 10), which agrees with the
theoretical expression of the coupling functions Eq. (57). Moreover, the signs
of the coefficients are also in accordance with the theoretical ones, e.g. for
the amplitude coupling we have a positive sign for the terms with
$\sigma^{0}_{i}\sigma_{j}^{0}=1$, $\sigma_{j}$, $\sigma_{j}^{2}$ and
$\sigma_{i}^{2}$, and a negative sign for the terms with $\sigma_{i}$,
$\sigma_{i}\sigma_{j}$. Similar agreement between the theoretical expression
of the coupling coefficients (Eq. (61)) and the reconstructed coupling
coefficients can be observed for the canonical model (Fig. 11). The coupling
coefficients for the other models are given in the Supplementary Materials as
well.

Figure 10: An example of the computed Fourier-Taylor coefficients of the
reconstructed coupling in the amplitude equation (top panels) and phase
equation (bottom panels) for two uni-directionally coupled radial isochron
clocks (see B.1).

Figure 11: An example of the computed Fourier-Taylor coefficients of the
reconstructed coupling in the amplitude equation (top panels) and phase
equation (bottom panels) for two uni-directionally coupled canonical models
(see B.2).

Figure 12: An example of the computed Fourier-Taylor coefficients of the
reconstructed coupling in the amplitude equation (top panels) and phase
equation (bottom panels) for two uni-directionally coupled van der Pol
oscillators (see B.3).

For the van der Pol oscillator the reconstructed coupling coefficients are
presented in Fig. 12. In this case, we do not have an analytic expression for
the transformation functions, and neither for the coupling part in the reduced
space. However, it is notable that the non-zero coefficients of the coupling
part appear in a grid form, which indicates that the basis functions can be
sparsified. Moreover, the grid form suggests that the coupling part can be
approximated (factorized) as a multiplication of two functions: the first one,
which consists of the basis functions that depend only on the oscillator’s own
variables (i.e. on $\sigma_{i}$ and $\varphi_{i}$, the basis functions on the
vertical axes) and the second function, which has basis functions that depend
on the input variables (i.e. on $\sigma_{j}$ and $\varphi_{j}$, the basis
functions on the horizontal axes). Thus, the former function corresponds
either to the amplitude or the phase response functions of the oscillator, and
the latter function to the input. The representation of the coupling part as a
multiplication of the response functions and the input is possible due to the
linear coupling of the original system for the van der Pol oscillators (see
Eq. (62)). Therefore, the coefficients of the coupling part for the Wilson-
Cowan models do not have a clear grid form (see Figs.
LABEL:smfig:coupl_ampl_wcmod and LABEL:smfig:coupl_phase_wcmod), since the
system non-linearly depends on the input (see Eq. (63)).

For all models the calculated coefficients of the coupling part of the first
oscillator (without input) were negligibly small, which demonstrates the
ability of the presented method to detect the directionality of the coupling.
The results of the estimated coupling coefficients for all models and
oscillators are given in the Supplementary Materials.

##  4 Discussion

We have presented a new method for the reconstruction of the phase-amplitude
dynamics from high time resolution signals measured from oscillatory systems.
This method is intended to be applied to electrophysiological signals from the
brain, such as Magneto- or Electroencephalography (M/EEG),
Electrocorticography (ECoG), and intracranial EEG (iEEG), to build a dynamical
system that reproduces the transient synchronization properties between brain
regions within and between different frequency bands observed in the
experiments.

The presented method combines recent advances in the field of phase-amplitude
reduction. Namely, the parametrization method [28, 29, 30, 31, 32, 33, 34] and
the time averaging approach for finding the transformation functions [38, 40,
39], based on the Koopman operator approach [35, 36, 37]. The approximation of
the parametrization functions from [31] allowed us to explain the connection
between the eigenfunctions of the Koopman operator and the transformation
functions. Moreover, using the approximation of the parametrization, we
derived a method of finding the result of the time averaging integral suitable
for numerical estimation (Eq. (33) and Fig. 3).

The results of the phase-amplitude reconstruction of the different models show
a good VF approximation both for the uncoupled system as well as for coupled
oscillatory systems. The approximated transformation functions and the
theoretical ones show a robust agreement for the coupled radial isochron
clocks and the coupled canonical models, for which they can be derived
analytically. Moreover, the approximated coupling coefficients agree with the
theoretical values. Application of our method to the two other models, namely
the van der Pol oscillators and the Wilson-Cowan models, with non-linear
interdependence of the phase and amplitude in the transformation functions
also yielded convincing results.

The reconstruction of the phase-amplitude dynamics in the reduced space aimed
at building an observable independent form of the system dynamics. As it was
shown in previous studies [16, 23], it is necessary to take into account the
non-circularity of the limit cycle, or equivalently, the asymmetry of the
phase distribution when reconstructing the reduced phase dynamics from
measured signals. Neglecting the transformation between observable and
theoretical phases could cause spurious coupling. The same principles apply
for phase-amplitude reduction dynamics reconstructed from the amplitude and
phase of measured signals. Indeed, the system in Eqs. (42) implies that in the
absence of coupling ($g^{(\varphi)}_{i,j}=0$ and $g^{(\sigma)}_{i,j}=0$) the
phase grows linearly and the amplitude decays exponentially. If we dismiss the
non-circularity of the limit cycle, e.g. by taking the experimentally measured
phase and defining the amplitude variable as the radial distance to the limit
cycle, the deviation of the observed amplitude decay from the exponential one
and the variation of the rotation speed of the observed phase along the limit
cycle and in the surrounding region will account for the coupling. Thus, the
reconstructed transformation functions serve the purpose of describing these
dependencies of the amplitude and the phase on the domain surrounding the
limit cycle.

As it was mentioned, the transformation functions are unique up to a constant
factor for amplitude and a constant shift for phase. In the presented method,
the constant factor and the shift are taken into account at the numerical
approximation of the transformation functions, thereby providing the
reconstruction of a unique and observable independent dynamical system in the
reduced space.

In comparison to many data-driven methods presented in recent studies (e.g.
[46, 54, 49]) we use in this study only polar coordinates and Fourier-Taylor
approximations. The crucial part of the presented method is finding the VF
with the coupling using basis functions of Fourier-Taylor series. The VF
approximation method resembles the ones presented in [16] for phase coupling
and in the sparse indentification of nonlinear dynamics (SINDY) method [60,
61], however, without sparsification of the VF coefficients matrix. The
representation of the VF as the sum of the basis functions has another
advantage. Albeit we do not make any assumptions about the VF of the uncoupled
system and assume that the coupling is pairwise, the approximation as Fourier-
Taylor basis functions allows us to consider the coupling part of the VF as an
additive (linear) input. This simplifies the calculation of the coupling
functions in the reduced space.

Another difference of our method is the approximation of the coupling part as
a sum of pairwise couplings, each one approximated using two dimensional
Fourier-Taylor series that depend on amplitude and phase of the driving and
driven systems in a general form. This is requested by the specificity of the
high time resolution of brain signals. Unlike the data-driven phase-amplitude
reduction methods based on dynamic mode decomposition (DMD) with input and
control [50, 51, 52, 53], SINDY [60, 61, 54], first and higher order
approximation of the isochron equations [45, 46, 47, 49] and the method of
estimation of the response functions [48], in our case we do not have access
to the input or driving signals coming into the region of interest. In
contrary, we can measure only the resulting electrical/electromagnetic signals
emitted by different regions of the brain. Therefore, we assume that input to
a region of the brain and the interaction of the neuronal populations within
the brain region produce the oscillations in different frequencies observed in
the measurements. Thus, the reconstruction method presented here is intended
to be used as a tool to describe the transient synchronization between these
oscillations within and between distant brain regions in the form of a
dynamical system. Our method thereby extends conventional methods, like phase
synchronization measures in the phase domain, correlation analyses and event-
related (de-)synchronization (ERD/ERS) in the amplitude domain, or Phase-
Amplitude Coupling (PAC) and similar methods in the phase and amplitude
domain, by providing not only information about the directionality of the
coupling, but also about the intrinsic mechanisms of the dynamical
interaction.

The presented reconstruction method has several requirements for application
that we next discuss in detail. (1) A good VF approximation in the domain
outside of the coverage of the data points is necessary and was possible due
to the application of the regularized least squares method with generalized
cross-validation (GCV) (see A). The regularization restricts the value of the
coefficients by means of the $L^{2}$ measure, thereby escaping an over-fitting
of the data points. However, for a successful reconstruction of the VF with
coupling a good coverage of the observable amplitude and phase domain is
required. From a practical point of view this requires a data set (trials)
that starts from possibly different amplitude and phase values. (2) Our method
can be applied to signals originated from any network of oscillatory systems,
if one could define a limit cycle for the considered oscillatory system. Thus,
clear oscillatory patterns should be observed in the measurements to be able
to interpret the results. (3) Our method is intended to be applied to the
oscillatory signals that are filtered from the actual measurements around a
given frequency. Thus, the noise in the signals should be filtered out as much
as possible in order to obtain a good approximation. (4) For the extraction of
the phase and amplitude information from the signal one should use adapted
methods that reduce the dependence of the extracted information from the
future and the past of the signal as discussed in [58, 59, 48].

In conclusion, we present in this work a method that allows us to derive a
dynamical system from high time resolution signals originated from a network
of interacting oscillators, in particular, from interacting distant brain
regions. The analysis of the resulting dynamical system can provide a useful
tool to construct mathematical models and investigate neural pathologies and
disorders of the brain, such as in silico experiments with targeted
manipulations of activity within brain regions or couplings between them. This
can contribute to a deeper understanding of the neurobiological underpinnings
of cognitive processes and how network synchronization is impaired by
(non-)pathological aging.

##  Appendix A Regularized least squares method

Here, we present the regularized least squares method in the ridge regression
formulation [62, 63, 64].

Let us assume that we have the column $N$-vectors
$\mathbf{x}=\\{x_{i}\\}^{T}$, $\mathbf{y}=\\{y_{i}\\}^{T}$, $i=1,\ldots,N$ and
basis functions $\psi_{j}:\mathbb{R}\rightarrow\mathbb{R}$, $j=1,\ldots,M$.
Let $\mathbf{\Psi}$ be the $N\times M$ matrix of basis functions with elements
$\psi_{j}(x_{i})$. We want to approximate $\mathbf{y}$ using the basis
functions as

$\mathbf{y}\approx\mathbf{\Psi}\mathbf{q},$ (49)

where $\mathbf{q}=\\{q_{j}\\}^{T}$ is the column $M$-vector of coefficients
that can be obtained by solving the minimization problem

$\mathbf{q}=\arg\min_{\mathbf{q}\in\mathbb{R}^{M}}\left\\{\|\mathbf{\Psi}\mathbf{q}-\mathbf{y}\|^{2}+\kappa\|\mathbf{q}\|^{2}\right\\},$
(50)

where $\kappa$ is the regularization parameter. Then, the solution of Eq. (50)
is

$\mathbf{q}=\left(\kappa\mathbf{I}+\mathbf{\Psi}^{T}\mathbf{\Psi}\right)^{-1}\mathbf{\Psi}^{T}\mathbf{y}.$
(51)

We find the optimal regularization parameter using the generalized cross-
validation method (GCV) [62, 63, 64]. The optimal $\kappa$ is found by
minimizing

$GCV(\kappa)=\frac{\|\mathbf{\Psi}\mathbf{q}-\mathbf{y}\|^{2}}{\tau^{2}},$
(52)

where $\tau$ is the effective number of the degree of freedom:

$\tau=N-\sum_{l=1}^{r}\frac{\varrho_{l}^{2}}{\varrho_{l}^{2}+\kappa},$ (53)

where $\varrho_{l}$ are the singular values of $\mathbf{\Psi}$, or
$\varrho_{l}^{2}$ are the eigenvalues of $\mathbf{\Psi}^{T}\mathbf{\Psi}$, and
$r$ is the rank of the matrices.

##  Appendix B Models

In all simulations models of two uni-directionally coupled oscillators were
used: the second oscillator receives an input from the first oscillator.

###  B.1 Radial isochron clock

The system of equations for the coupled radial isochron clock in polar
coordinates reads as follows [65, 66]

$\displaystyle\frac{d\theta_{i}}{dt}$
$\displaystyle=1+\epsilon_{ij}\frac{r_{j}}{r_{i}}\cos(\theta_{i})\sin(\theta_{j}),$
(54) $\displaystyle\frac{dr_{i}}{dt}$
$\displaystyle=a_{i}r_{i}(1-r_{i}^{2})+\epsilon_{ij}\sin(\theta_{i})r_{j}\sin(\theta_{j}),$

where $a_{i}>0$ is a parameter controlling the decay rate to the limit cycle
and $\epsilon_{ij}$ is a parameter controlling the coupling strength.

The theoretical expressions for the forward transformation functions are

$\displaystyle K^{(\theta)}_{i}(\varphi_{i},\sigma_{i})=$
$\displaystyle\varphi_{i},$ (55) $\displaystyle
K^{(r)}_{i}(\varphi_{i},\sigma_{i})=$
$\displaystyle\sqrt{\frac{1}{1-2a_{i}\sigma_{i}}},$

and for the inverse transformation functions are

$\displaystyle\Phi_{i}(\theta_{i},r_{i})=$ $\displaystyle\theta_{i},$ (56)
$\displaystyle\Sigma_{i}(\theta_{i},r_{i})=$
$\displaystyle\frac{1}{2a_{i}}\frac{r_{i}^{2}-1}{r_{i}^{2}}.$

The system of equations for the coupled radial isochron clock in the reduced
space approximated using Fourier-Taylor series on both amplitude variables
$\sigma_{i}$ and $\sigma_{j}$ is

$\frac{d\varphi_{i}}{dt}\approx
1+\epsilon_{ij}\cos(\varphi_{i})\sin(\varphi_{j})\Bigl{(}1-a_{i}\sigma_{i}+a_{j}\sigma_{j}-\frac{1}{2}a_{i}^{2}\sigma_{i}^{2}\\\
-a_{i}a_{j}\sigma_{i}\sigma_{j}+\frac{3}{2}a_{j}^{2}\sigma_{j}^{2}+\ldots\Bigr{)},\\\
\frac{d\sigma_{i}}{dt}\approx-2a_{i}\sigma_{i}+\epsilon_{ij}\frac{\sin(\varphi_{i})\sin(\varphi_{j})}{a_{i}}\Bigl{(}1-3a_{i}\sigma_{i}+a_{j}\sigma_{j}\\\
+\frac{3}{2}a_{i}^{2}\sigma_{i}^{2}-3a_{i}a_{j}\sigma_{i}\sigma_{j}+\frac{3}{2}a_{j}^{2}\sigma_{j}^{2}+\ldots\Bigr{)},$
(57)

where we have provided the terms up to order 2 in the amplitude variables
$\sigma$. For the derivation of these equations we refer the reader to the
Supplementary Materials.

###  B.2 Canonical model

The system of equations for the canonical model for an oscillator, also known
as Stuart-Landau oscillator [67, 68, 25], in polar coordinates reads as
follows

$\displaystyle\dot{\theta_{i}}=$ $\displaystyle
1+\alpha_{i}a_{i}r_{i}^{2}-\epsilon_{ij}\frac{r_{j}}{r_{i}}\sin(\theta_{i})\cos(\theta_{j}),$
(58) $\displaystyle\dot{r_{i}}=$
$\displaystyle\alpha_{i}r_{i}(1-r_{i}^{2})+\epsilon_{ij}r_{j}\cos(\theta_{i})\cos(\theta_{j}),$

where oscillator $i$ receives input from oscillator $j$. The parameter values
used in the simulations are: $a_{1}=1.2$, $a_{2}=1.0$, $\alpha_{1}=1.5$,
$\alpha_{2}=2.0$, $\epsilon_{12}=0.0$, $\epsilon_{21}=0.3$.

The theoretical expressions for the forward transformation functions are

$\displaystyle K^{(\theta)}_{i}(\varphi_{i},\sigma_{i})=$
$\displaystyle\varphi_{i}+\frac{a_{i}}{2}\ln(1-2\alpha_{i}\sigma_{i}),$ (59)
$\displaystyle K^{(r)}_{i}(\varphi_{i},\sigma_{i})=$
$\displaystyle\sqrt{\frac{1}{1-2\alpha_{i}\sigma_{i}}},$

and for the inverse transformation functions are

$\displaystyle\Phi_{i}(\theta_{i},r_{i})=$
$\displaystyle\theta_{i}+a_{i}\ln(r_{i}),$ (60)
$\displaystyle\Sigma_{i}(\theta_{i},r_{i})=$
$\displaystyle\frac{1}{2\alpha_{i}}\left(1-\frac{1}{r_{i}^{2}}\right).$

The system of equations for the canonical model in the reduced space using
Fourier-Taylor series on both amplitude variables $\sigma_{i}$ and
$\sigma_{j}$ is

$\displaystyle\frac{d\sigma_{i}}{dt}=$
$\displaystyle-2\alpha_{i}\sigma_{i}+\epsilon_{ij}\Bigl{(}(\cos(\varphi_{i})\cos(\varphi_{j}))/\alpha_{i}$
$\displaystyle+(\sin(\varphi_{i})\cos(\varphi_{j})a_{i}-3\cos(\varphi_{i})\cos(\varphi_{j}))\sigma_{i}$
$\displaystyle+(\cos(\varphi_{i})\sin(\varphi_{j})a_{j}+\cos(\varphi_{i})\cos(\varphi_{j}))\alpha_{j}\sigma_{j}/\alpha_{i}$
$\displaystyle-((a_{i}^{2}-3)\cos(\varphi_{i})\cos(\varphi_{j})+4a_{i}\sin(\varphi_{i})\cos(\varphi_{j}))\alpha_{i}\sigma_{i}^{2}/2$
$\displaystyle+(\sin(\varphi_{i})\sin(\varphi_{j})a_{i}a_{j}+\sin(\varphi_{i})\cos(\varphi_{j})a_{i}$
$\displaystyle-3\cos(\varphi_{i})\sin(\varphi_{j})a_{j}-3\cos(\varphi_{i})\cos(\varphi_{j}))\alpha_{j}\sigma_{j}\sigma_{i}$
$\displaystyle-((a_{j}^{2}-3)\cos(\varphi_{i})\cos(\varphi_{j})-4a_{j}\cos(\varphi_{i})\sin(\varphi_{j}))\alpha_{j}^{2}\sigma_{j}^{2}/(2\alpha_{i})$
$\displaystyle+\ldots\Bigr{)}.$ $\displaystyle\frac{d\varphi_{i}}{dt}=$
$\displaystyle
1+\epsilon_{ij}\Bigl{(}-\sin(\varphi_{i})\cos(\varphi_{j})+\cos(\varphi_{i})\cos(\varphi_{j})a_{i}$
$\displaystyle+(a_{i}^{2}+1)\sin(\varphi_{i})\cos(\varphi_{j})\alpha_{i}\sigma_{i}$
$\displaystyle+\bigl{[}a_{i}a_{j}\cos(\varphi_{i})\sin(\varphi_{j})+a_{i}\cos(\varphi_{i})\cos(\varphi_{j})$
$\displaystyle-
a_{j}\sin(\varphi_{i})\sin(\varphi_{j})-\sin(\varphi_{i})\cos(\varphi_{j})\bigr{]}\alpha_{j}\sigma_{j}$
$\displaystyle+\bigl{[}-\cos(\varphi_{i})\cos(\varphi_{j})a_{i}^{3}+\sin(\varphi_{i})\cos(\varphi_{j})a_{i}^{2}$
$\displaystyle-\cos(\varphi_{i})\cos(\varphi_{j})a_{i}+\sin(\varphi_{i})\cos(\varphi_{j})\bigr{]}\alpha_{i}^{2}\sigma_{i}^{2}/2$
$\displaystyle+\bigl{[}a_{i}^{2}a_{j}\sin(\varphi_{i})\sin(\varphi_{j})+a_{i}^{2}\sin(\varphi_{i})\cos(\varphi_{j})$
$\displaystyle+a_{j}\sin(\varphi_{i})\sin(\varphi_{j})+\sin(\varphi_{i})\cos(\varphi_{j})\bigr{]}\alpha_{i}\alpha_{j}\sigma_{j}\sigma_{i}$
$\displaystyle-\bigl{[}(a_{j}^{2}-3)a_{i}\cos(\varphi_{i})\cos(\varphi_{j})+(3-a_{j}^{2})\sin(\varphi_{i})\cos(\varphi_{j})$
$\displaystyle-4a_{i}a_{j}\cos(\varphi_{i})\sin(\varphi_{j})+4a_{j}\sin(\varphi_{i})\sin(\varphi_{j})\bigr{]}\alpha_{j}^{2}\sigma_{j}^{2}/2+\ldots\Bigr{)},$
(61)

where we have provided the terms up to order 2 in the amplitude variables
$\sigma$. For the derivation of the equations we refer the reader to the
Supplementary Materials.

###  B.3 Van der Pol oscillator

The system of equations for the van der Pol oscillator in Cartesian
coordinates [69, 70, 66] reads as follows

$\displaystyle\frac{dx_{i}}{dt}=$ $\displaystyle y_{i},$ (62)
$\displaystyle\frac{dy_{i}}{dt}=$
$\displaystyle\mu_{i}(1-x_{i}^{2})y_{i}-x_{i}+\epsilon_{ij}x_{j},$

where an oscillator $i$ receives an input from an oscillator $j$. The
parameter values used in the simulations: $\mu_{1}=0.3$, $\mu_{2}=0.5$,
$\epsilon_{12}=0.0$, $\epsilon_{21}=0.1$.

###  B.4 Wilson-Cowan model

The system of equations for the Wilson-Cowan model [71, 24] in Cartesian
coordinates reads

$\displaystyle\frac{dx_{i}}{dt}=$ $\displaystyle-
x_{i}+S(\rho_{x,i}+a_{i}x_{i}-b_{i}y_{i}),$ (63)
$\displaystyle\frac{dy_{i}}{dt}=$ $\displaystyle-
y_{i}+S(\rho_{y,i}+c_{i}x_{i}-d_{i}y_{i}+\epsilon_{ij}(x_{j}-x_{0,j})),$

where oscillator $i$ receives input from oscillator $j$, $x_{0,j}$ is the
$x$-component of the unstable equilibrium point for the WC model corresponding
to uncoupled oscillator $j$, that can be found by solving the equation:

$\displaystyle 0=$ $\displaystyle-
x_{0,j}+S(\rho_{x,j}+a_{j}x_{0,j}-b_{j}y_{0,j}),$ (64) $\displaystyle 0=$
$\displaystyle-y_{0,j}+S(\rho_{y,j}+c_{j}x_{0,j}-d_{j}y_{0,j}).$ (65)

The parameter values used in the simulations:
$a_{1}=a_{2}=b_{1}=b_{2}=c_{1}=c_{2}=10$, $d_{1}=d_{2}=-2$,
$\rho_{x,1}=\rho_{x,2}=0.0$, $\rho_{y,1}=-6.75$, $\rho_{y,2}=-7.0$,
$\epsilon_{12}=0.0$, $\epsilon_{21}=1.0$.

##  Appendix C Finding initial conditions in the reduced space

In general, at the limit $T\rightarrow\infty$ the expression for $R(t)$ in Eq.
(31) is oscillatory as shown in Fig. 13 for the van der Pol oscillator. The
system is integrated for five periods, enough time for the trajectory to relax
to the limit cycle (upper panel). The red line in the middle panel of Fig. 13
marks the time course of $\log(|\rho(t)|)-\lambda t$, $t\in[t_{1},t_{2}]$,
where $t_{2}$ is the time point at the threshold of numerical precision and
$t_{1}=t_{2}-2T$. If the initial approximation of $\lambda$ found by using the
monodromy matrix and the solution of Eq. (17) is not precise the oscillations
of $\log(|\rho(t)|)-\lambda t$ at the limit have a non-zero slope, as shown in
the middle panel in Fig. 13 (here for demonstration purposes the initial
approximation of $\lambda$ was changed). The correction $\Delta\lambda$ to the
value of $\lambda$ can be estimated as the slope of the line, that passes
through the points $\log(|\rho(t_{1})|)-\lambda t_{1}$ and
$\log(|\rho(t_{2})|)-\lambda t_{2}$. After correction of the value of
$\lambda$ the time course of $|\rho(t)|e^{-(\lambda+\Delta\lambda)t}$ tends to
an oscillation with zero slope and mean equal to $|q^{(r)}_{1,0}\sigma_{0}|$.

Figure 13: An example of finding the correction to the value of $\lambda$ and
the initial condition $\sigma_{0}$ for the van der Pol oscillator. Upper
panel: the time course of an observable
$\rho(t)=r(t)-\gamma^{(r)}_{0}(\theta(t))$ (solid blue line) that vanishes on
the limit cycle $\Gamma$. Middle panel: the expression
$\log(|\rho(t)|)-\lambda t$ (solid blue line) relaxes to a periodic
oscillation. The slope of the oscillation $\Delta\lambda$ between the black
dashed line and the horizontal black dotted line gives us the correction to
$\lambda$. The two vertical dot-dashed lines are in two-period-distance to
each other ($t_{2}=t_{1}+2T$). The red line marks the time course within
$t\in[t_{1},t_{2}]$. The black star and black circle are the first and the
last point of this interval, respectively. Vertical dotted lines denote the
periods of the oscillation. Bottom panel: The blue line is the time course of
$|\rho(t)|e^{-(\lambda+\Delta\lambda)t}$. The vertical dotted lines and the
two vertical dot-dashed lines ($t_{1}$ and $t_{2}$) are the same as in the
middle panel. The time course within the interval $t\in[t_{1},t_{2}]$ is
marked by the red line, the mean of which is equal to
$|q^{(r)}_{1,0}\sigma_{0}|$ (red star). The horizontal red dotted line is
$\langle|\sigma_{0}\hat{R}(t)|\rangle_{T}$.

## References

  * [1] F. Varela, J. P. Lachaux, E. Rodriguez, and J. Martinerie, “The brainweb: phase synchronization and large-scale integration.,” Nature reviews. Neuroscience, vol. 2, no. 4, pp. 229–239, 2001.
  * [2] G. Buzsáki, D. Andreas, A. Draguhn, G. Buzsaki, and A. Draguhn, “Neuronal Oscillations in Cortical Networks,” Science, vol. 304, p. 1926, jun 2004.
  * [3] F. van Ede, A. J. Quinn, M. W. Woolrich, and A. C. Nobre, “Neural Oscillations: Sustained Rhythms or Transient Burst-Events?,” Trends in Neurosciences, vol. 41, pp. 415–417, jul 2018.
  * [4] a. K. Engel, P. Fries, and W. Singer, “Dynamic predictions: oscillations and synchrony in top-down processing.,” Nature reviews. Neuroscience, vol. 2, no. 10, pp. 704–16, 2001.
  * [5] R. E. Greenblatt, M. E. Pflieger, and a. E. Ossadtchi, “Connectivity measures applied to human brain electrophysiological data,” Journal of Neuroscience Methods, vol. 207, no. 1, pp. 1–16, 2012.
  * [6] F. Lopes da Silva, “EEG and MEG: Relevance to neuroscience,” Neuron, vol. 80, no. 5, pp. 1112–1128, 2013.
  * [7] J.-p. Lachaux, E. Rodriguez, J. Martinerie, and F. J. Varela, “Measuring phase synchrony in brain signals,” Human brain mapping, vol. 8, pp. 194–208, 1999.
  * [8] P. Sauseng, W. Klimesch, W. R. Gruber, and N. Birbaumer, “Cross-frequency phase synchronization: A brain mechanism of memory matching and attention,” NeuroImage, vol. 40, no. 1, pp. 308–317, 2008.
  * [9] P. Sauseng and W. Klimesch, “What does phase information of oscillatory brain activity tell us about cognitive processes?,” Neuroscience and Biobehavioral Reviews, vol. 32, no. 5, pp. 1001–1013, 2008.
  * [10] G. Pfurtscheller and F. H. Lopes Da Silva, “Event-related EEG/MEG synchronization and desynchronization: Basic principles,” Clinical Neurophysiology, vol. 110, no. 11, pp. 1842–1857, 1999.
  * [11] A. B. L. Tort, R. Komorowski, H. Eichenbaum, and N. Kopell, “Measuring Phase-Amplitude Coupling Between Neuronal Oscillations of Different Frequencies,” Journal of Neurophysiology, vol. 104, no. 2, pp. 1195–1210, 2010.
  * [12] B. Voytek, M. D’Esposito, N. Crone, and R. T. Knight, “A method for event-related phase/amplitude coupling,” NeuroImage, vol. 64, pp. 416–424, jan 2013.
  * [13] M. J. Hülsemann, E. Naumann, and B. Rasch, “Quantification of Phase-Amplitude Coupling in Neuronal Oscillations: Comparison of Phase-Locking Value, Mean Vector Length, Modulation Index, and Generalized-Linear-Modeling-Cross-Frequency-Coupling,” Frontiers in Neuroscience, vol. 13, pp. 1–15, jun 2019.
  * [14] M. G. Rosenblum and a. S. Pikovsky, “Detecting direction of coupling in interacting oscillators.,” Physical review. E, Statistical, nonlinear, and soft matter physics, vol. 64, no. 4 Pt 2, p. 045202, 2001.
  * [15] B. Kralemann, L. Cimponeriu, M. Rosenblum, A. Pikovsky, and R. Mrowka, “Uncovering interaction of coupled oscillators from data,” Physical Review E, vol. 76, pp. 1–4, nov 2007.
  * [16] B. Kralemann, L. Cimponeriu, M. Rosenblum, A. Pikovsky, and R. Mrowka, “Phase dynamics of coupled oscillators reconstructed from data,” Physical Review E, vol. 77, pp. 1–16, jun 2008.
  * [17] B. Kralemann, A. Pikovsky, and M. Rosenblum, “Reconstructing phase dynamics of oscillator networks,” Chaos, vol. 21, no. 2, pp. 1–10, 2011.
  * [18] B. Kralemann, A. Pikovsky, and M. Rosenblum, “Detecting triplet locking by triplet synchronization indices,” Physical Review E - Statistical, Nonlinear, and Soft Matter Physics, vol. 87, no. 5, pp. 1–6, 2013.
  * [19] B. Kralemann, A. Pikovsky, and M. Rosenblum, “Reconstructing effective phase connectivity of oscillator networks from observations,” New Journal of Physics, vol. 16, p. 085013, aug 2014.
  * [20] W. D. Penny, V. Litvak, L. Fuentemilla, E. Duzel, and K. Friston, “Dynamic Causal Models for phase coupling,” Journal of Neuroscience Methods, vol. 183, no. 1, pp. 19–30, 2009.
  * [21] T. Stankovski, T. Pereira, P. V. E. McClintock, and A. Stefanovska, “Coupling functions: Universal insights into dynamical interaction mechanisms,” Reviews of Modern Physics, vol. 89, p. 045001, nov 2017.
  * [22] T. Stankovski, T. Pereira, P. V. McClintock, and A. Stefanovska, “Coupling functions: Dynamical interaction mechanisms in the physical, biological and social sciences,” Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 377, no. 2160, 2019.
  * [23] A. Yeldesbay, G. R. Fink, and S. Daun, “Reconstruction of effective connectivity in the case of asymmetric phase distributions,” Journal of Neuroscience Methods, vol. 317, pp. 94–107, apr 2019.
  * [24] F. C. Hoppensteadt and E. M. Izhikevich, Weakly Connected Neural Networks, vol. 126 of Applied Mathematical Sciences.  New York, NY: Springer New York, 1997.
  * [25] A. Pikovsky, M. Rosenblum, and J. Kurths, Synchronization: A Universal Concept in Nonlinear Sciences.  Cambridge nonlinear science series, Cambridge University Press, 2001.
  * [26] A. Pikovsky and M. Rosenblum, “Dynamics of globally coupled oscillators: Progress and perspectives.,” Chaos (Woodbury, N.Y.), vol. 25, no. 9, p. 097616, 2015.
  * [27] B. Monga, D. Wilson, T. Matchen, and J. Moehlis, “Phase reduction and phase-based optimal control for biological systems: a tutorial,” Biological Cybernetics, vol. 113, no. 1-2, pp. 11–46, 2019.
  * [28] X. Cabre, E. Fontich, and R. de la Llave, “The parameterization method for invariant manifolds I: Manifolds associated to non-resonant subspaces,” Indiana University Mathematics Journal, vol. 52, no. 2, pp. 283–328, 2003\.
  * [29] X. Cabre, E. Fontich, and R. de la Llave, “The parameterization method for invariant manifolds II: regularity with respect to parameters,” Indiana University Mathematics Journal, vol. 52, no. 2, pp. 329–360, 2003.
  * [30] X. Cabré, E. Fontich, and R. de la Llave, “The parameterization method for invariant manifolds III: Overview and applications,” Journal of Differential Equations, vol. 218, no. 2, pp. 444–515, 2005.
  * [31] A. Guillamon and G. Huguet, “A Computational and Geometric Approach to Phase Resetting Curves and Surfaces,” SIAM Journal on Applied Dynamical Systems, vol. 8, no. 3, pp. 1005–1042, 2009.
  * [32] G. Huguet and R. De La Llave, “Computation of limit cycles and their isochrons: Fast algorithms and their convergence,” SIAM Journal on Applied Dynamical Systems, vol. 12, no. 4, pp. 1763–1802, 2013.
  * [33] A. Pérez-Cervera, P. Ashwin, G. Huguet, T. M. Seara, and J. Rankin, “The uncoupling limit of identical Hopf bifurcations with an application to perceptual bistability,” The Journal of Mathematical Neuroscience, vol. 9, p. 7, dec 2019.
  * [34] A. Pérez-Cervera, T. M-Seara, and G. Huguet, “Global phase-amplitude description of oscillatory dynamics via the parameterization method,” Chaos, vol. 30, no. 8, 2020.
  * [35] I. Mezić and A. Banaszuk, “Comparison of systems with complex behavior,” Physica D: Nonlinear Phenomena, vol. 197, no. 1-2, pp. 101–133, 2004.
  * [36] I. Mezic, “Spectral properties of dynamical systems, model reduction and decompositions,” Nonlinear Dynamics, vol. 41, pp. 309–325, 2005.
  * [37] M. Budišić, R. Mohr, and I. Mezić, “Applied Koopmanism,” Chaos, vol. 22, no. 4, 2012.
  * [38] A. Mauroy and I. Mezić, “On the use of Fourier averages to compute the global isochrons of (quasi)periodic dynamics,” Chaos, vol. 22, no. 3, 2012\.
  * [39] A. Mauroy and I. Mezić, “Global computation of phase-amplitude reduction for limit-cycle dynamics,” Chaos, vol. 28, no. 7, 2018.
  * [40] A. Mauroy, I. Mezić, and J. Moehlis, “Isostables, isochrons, and Koopman spectrum for the action-angle representation of stable fixed point dynamics,” Physica D: Nonlinear Phenomena, vol. 261, pp. 19–30, oct 2013\.
  * [41] D. Wilson and J. Moehlis, “Isostable reduction of periodic orbits,” Physical Review E, vol. 94, p. 052213, nov 2016.
  * [42] D. Wilson and B. Ermentrout, “Greater accuracy and broadened applicability of phase reduction using isostable coordinates,” Journal of Mathematical Biology, vol. 76, no. 1-2, pp. 37–66, 2018.
  * [43] D. Wilson, “Phase-amplitude reduction far beyond the weakly perturbed paradigm,” Physical Review E, vol. 101, no. 2, pp. 1–21, 2020.
  * [44] D. Wilson, “Analysis of input-induced oscillations using the isostable coordinate framework,” Chaos, vol. 31, no. 2, pp. 1–15, 2021.
  * [45] D. Wilson, “A data-driven phase and isostable reduced modeling framework for oscillatory dynamical systems,” Chaos, vol. 30, no. 1, 2020.
  * [46] D. Wilson, “Data-driven inference of high-accuracy isostable-based dynamical models in response to external inputs,” Chaos: An Interdisciplinary Journal of Nonlinear Science, vol. 31, pp. 1–33, jun 2021.
  * [47] T. Ahmed, A. Sadovnik, and D. Wilson, “Data-driven inference of low-order isostable-coordinate-based dynamical models using neural networks,” Nonlinear Dynamics, 2022.
  * [48] R. Cestnik, E. T. K. Mau, and M. Rosenblum, “Inferring oscillator’s phase and amplitude response from a scalar signal exploiting test stimulation,” New Journal of Physics, vol. 24, p. 123012, dec 2022.
  * [49] D. Wilson, “A direct method approach for data-driven inference of high accuracy adaptive phase-isostable reduced order models,” Physica D: Nonlinear Phenomena, vol. 446, p. 133675, 2023.
  * [50] J. H. Tu, C. W. Rowley, D. M. Luchtenburg, S. L. Brunton, and J. N. Kutz, “On dynamic mode decomposition: Theory and applications,” Journal of Computational Dynamics, vol. 1, no. 2, pp. 391–421, 2014.
  * [51] M. O. Williams, I. G. Kevrekidis, and C. W. Rowley, “A Data–Driven Approximation of the Koopman Operator: Extending Dynamic Mode Decomposition,” Journal of Nonlinear Science, vol. 25, pp. 1307–1346, dec 2015.
  * [52] J. L. Proctor, S. L. Brunton, and J. N. Kutz, “Dynamic mode decomposition with control,” SIAM Journal on Applied Dynamical Systems, vol. 15, no. 1, pp. 142–161, 2016.
  * [53] J. L. Proctor, S. L. Brunton, and J. Nathan Kutz, “Generalizing koopman theory to allow for inputs and control,” SIAM Journal on Applied Dynamical Systems, vol. 17, no. 1, pp. 909–930, 2018.
  * [54] E. Kaiser, J. N. Kutz, and S. L. Brunton, “Data-driven discovery of Koopman eigenfunctions for control,” Machine Learning: Science and Technology, vol. 2, p. 035023, sep 2021.
  * [55] E. D. Fagerholm, R. J. Moran, I. R. Violante, R. Leech, and K. J. Friston, “Dynamic causal modelling of phase-amplitude interactions,” NeuroImage, vol. 208, no. December 2019, 2020.
  * [56] A. Mauroy, Y. Susuki, and I. Mezić, Introduction to the Koopman Operator in Dynamical Systems and Control Theory, pp. 3–33.  Cham: Springer International Publishing, 2020.
  * [57] O. Castejón, A. Guillamon, and G. Huguet, “Phase-amplitude response functions for transient-state stimuli.,” Journal of mathematical neuroscience, vol. 3, no. 1, p. 13, 2013.
  * [58] M. Rosenblum, A. Pikovsky, A. A. Kühn, and J. L. Busch, “Real-time estimation of phase and amplitude with application to neural data,” Scientific Reports, vol. 11, no. 1, pp. 1–11, 2021.
  * [59] J. L. Busch, L. K. Feldmann, A. A. Kühn, and M. Rosenblum, “Real-time phase and amplitude estimation of neurophysiological signals exploiting a non-resonant oscillator,” Experimental Neurology, vol. 347, no. August 2021, p. 113869, 2022.
  * [60] S. L. Brunton, J. L. Proctor, and J. N. Kutz, “Discovering governing equations from data by sparse identification of nonlinear dynamical systems,” Proceedings of the National Academy of Sciences of the United States of America, vol. 113, no. 15, pp. 3932–3937, 2016.
  * [61] S. Brunton and J. Kutz, Data-Driven Science and Engineering: Machine Learning, Dynamical Systems, and Control.  Cambridge University Press, 2019.
  * [62] G. H. Golub, M. Heath, and G. Wahba, “Generalized cross-validation as a method for choosing a good ridge parameter,” Technometrics, vol. 21, no. 2, pp. 215–223, 1979.
  * [63] D. M. Bates and G. Wahba, “A truncated singular value decomposition and other methods for generalized cross-validation,” Technical Report, vol. 715, 1983.
  * [64] C. Fenu, L. Reichel, G. Rodriguez, and H. Sadok, “GCV for Tikhonov regularization by partial SVD,” BIT Numerical Mathematics, vol. 57, pp. 1019–1039, dec 2017.
  * [65] A. Winfree, The Geometry of Biological Time.  Biomathematics (Berlin), Springer Verlag, 1980.
  * [66] S. Strogatz, Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry and Engineering.  Studies in nonlinearity, Westview, 2000.
  * [67] L. D. Landau, “On the Problem of Turbulence,” Compt. Rend. Acad. Sci. URSS, vol. 44, 1944.
  * [68] J. T. Stuart, “On the non-linear mechanics of hydrodynamic stability,” Journal of Fluid Mechanics, vol. 4, no. 1, p. 1–21, 1958.
  * [69] B. van der Pol, “A theory of the amplitude of free and forced triode vibrations,” Radio Review, vol. 1, pp. 701–710, 754–762, 1920.
  * [70] N. N. Bogoliubov and Y. A. Mitropolski, Asymptotic Methods in the Theory of Non-Linear Oscillations.  London: Gordon and Breach, 1961.
  * [71] H. R. Wilson and J. D. Cowan, “Excitatory and Inhibitory Interactions in Localized Populations of Model Neurons,” Biophysical Journal, vol. 12, no. 1, pp. 1–24, 1972.