# Inferring oscillator’s phase and amplitude response from a scalar signal
exploiting test stimulation

Rok Cestnik  Erik T. K. Mau  Michael Rosenblum  Department of Physics and
Astronomy, University of Potsdam, Karl-Liebknecht-Str. 24/25, D-14476 Potsdam-
Golm, Germany

(August 28, 2024)

###### Abstract

The phase sensitivity curve or phase response curve (PRC) quantifies the
oscillator’s reaction to stimulation at a specific phase and is a primary
characteristic of a self-sustained oscillatory unit. Knowledge of this curve
yields a phase dynamics description of the oscillator for arbitrary weak
forcing. Similar, though much less studied characteristic, is the amplitude
response that can be defined either using an ad hoc approach to amplitude
estimation or via the isostable variables. Here, we discuss the problem of the
phase and amplitude response inference from observations using test
stimulation. Although PRC determination for noise-free neuronal-like
oscillators perturbed by narrow pulses is a well-known task, the general case
remains a challenging problem. Even more challenging is the inference of the
amplitude response. This characteristic is crucial, e.g., for controlling the
amplitude of the collective mode in a network of interacting units; this task
is relevant for neuroscience. Here, we compare the performance of different
techniques suitable for inferring the phase and amplitude response,
particularly with application to macroscopic oscillators. We suggest
improvements to these techniques, e.g., demonstrating how to obtain the PRC in
case of stimuli of arbitrary shape. Our main result is a novel technique
denoted by IPID-1, based on the direct reconstruction of the Winfree equation
and the analogous first-order equation for isostable dynamics. The technique
works for signals with or without well-pronounced marker events and pulses of
arbitrary shape; in particular, we consider charge-balanced pulses typical in
neuroscience applications. Moreover, this technique is superior for noisy and
high-dimensional systems. Additionally, we describe an error measure that can
be computed solely from data and complements any inference technique.

phase response, amplitude response, phase – isostable reduction, inference

##  I Introduction

Analysis and control of real-world oscillatory dynamics require inference of
oscillator’s parameters from observations. In particular, one can explore the
system by applying a specifically designed perturbation and measuring the
reaction. This paper discusses experiments extracting a self-sustained unit’s
phase and amplitude response. We consider oscillators with weakly-stable limit
cycles and stimulation with pulses of arbitrary shape.

In the first approximation in perturbation’s strength, the oscillator’s phase
response is quantified by the phase sensitivity function or phase response
curve (PRC) Winfree (1980); Mackey and Glass (1988); Rinzel and Ermentrout
(1998); Canavier (2006), a crucial characteristic of a self-sustained
oscillator. Knowing the PRC one describes the oscillator’s phase dynamics via
the Winfree equation

$\dot{\varphi}=\omega+Z(\varphi)p(t)\;,$ (1)

where $\varphi$ is the phase, $\omega$ is the natural frequency, $Z(\varphi)$
is the PRC, and $p(t)$ is the external perturbation. Equation (1) predicts,
e.g., the system’s response to an arbitrary stimulus or synchronization of the
oscillator by an external force. A description of oscillatory systems in terms
of PRCs is useful in various fields, e.g., in computational neuroscience
Achuthan and Canavier (2007).

For known oscillator’s equations, one can obtain PRC by solving the adjoint
problem Ermentrout (2002). Estimation of PRC in an experiment is a less
trivial task. The standard approach is to apply weak, narrow pulses at
different phases and measure the phase shifts caused by stimulation Mackey and
Glass (1988); Canavier (2006). Indeed, if the pulse at phase $\varphi$ is the
Dirac’s delta function, then the induced phase shift equals precisely
$Z(\varphi)$. It is easy to implement this approach for neuronal oscillators,
where the distance between the spikes gives the variation of the period and,
hence, the phase shift. It is not that easy to exploit this technique to
investigate oscillators without well-pronounced marker events that can be
assigned a specific phase value. Furthermore, for any oscillator, the applied
pulses should reasonably approximate the Dirac’s delta pulses, which is often
impossible because stimulation of live tissue shall fulfill specific criteria,
namely to be charge-balanced.

A related, much more demanding problem is the inference of the amplitude
response of an oscillator. This problem naturally arises, e.g., in model
studies Tass (1999); *Tass-00; *Tass-01; *Tass_2001; *Tass-02; *Rosenblum-
Pikovsky-04; *Rosenblum-Pikovsky-04a; *Popovych-Hauptmann-Tass-05; *Tukhlina-
Rosenblum-Pikovsky-Kurths-07; *Hauptmann-Tass-09; *wilson2011; *Popovych-
Tass-12; *Lin_2013; *Zhou_2017; *Wilson-Moehlis-16; *Popovych_et_al-17;
*Krylov-Dylov-Rosenblum-20; Montaseri  _et al._ (2013); Holt  _et al._ (2016);
Rosenblum (2020); Duchet  _et al._ (2020) of a clinical technique known as
deep brain stimulation (DBS) Benabid  _et al._ (1991); *Benabid_et_al-09;
*Kuehn-Volkmann-17. DBS aims to modulate a brain rhythm by stimuli delivered
through implanted microelectrodes. Thus, a question arises: stimulation at
which phase results in the most significant change of the amplitude Holt and
Netoff (2014); Holt  _et al._ (2016); Duchet  _et al._ (2020).

This paper critically assesses previously developed techniques and proposes
two new approaches for estimating an oscillator’s phase and amplitude response
from observations, exploiting test pulses of arbitrary shapes. The first novel
technique relies on the instantaneous amplitude and phase of a signal and,
hence, is in the spirit of traditional time series analysis. The second
technique is model-based; it exploits recent advances in dynamical systems’
description and reconstructs the equations for phase and, for the first time,
for the isostable variable. We exploit several test systems to probe our
techniques’ performance and compare them with those known from the literature.
In particular, we concentrate on systems exhibiting signals without well-
pronounced marker events, i.e., non-spiky signals, where the standard
technique is not efficient.

The paper is organized as follows. In the rest of this section, we discuss the
state of the art and problem formulation in detail. In Section II we present
our test models. Section III.4 presents the signal analysis-based techniques
and their comparison, while Section III.5 introduces the reconstruction of the
phase – isostables dynamics and demonstrates the performance of this approach.
Section IV summarizes and discusses the results. Appendix section introduces
the phase – isostables representation and presents the details of the test
systems and inference techniques.

###  I.1 State of the art and problems to be solved

Here we briefly review the currently known approaches to estimate the phase
and amplitude response in experiments.

####  I.1.1 Estimating phase response in experiments

PRC estimation for neuron-like oscillators with well-defined spikes is well
known, e.g., see Refs Mackey and Glass (1988); Canavier (2006). Indeed,
without loss of generality, we can assign zero phase value (or $2\pi$) to the
spikes. Let the unperturbed period (inter-spike interval) be $T_{0}$. Suppose
a stimulus perturbs the oscillator at time $t$ after a spike, $0\leq t<T_{0}$;
the phase at the stimulus’s onset is $\varphi=2\pi t/T_{0}$. The phase shift
caused by the stimulus is $\Delta\varphi=2\pi(T-T_{0})/T_{0}$, where $T$ is
the perturbed inter-spike interval. The dependence
$\Delta\varphi=\Delta\varphi(\varphi)$ obtained by multiple stimulations at
different phases yields the estimation of the PRC, $Z(\varphi)$.

However, this simple and reliable technique quantifies the response to pulses
of the particular shape used for stimulation. This means, if the pulse is
narrow and represents a good approximation to the Dirac’s delta pulse, then
the obtained response is indeed an estimation of the PRC $Z(\varphi)$.
However, if the pulse’s shape is essentially different, then the response
function also differs from $Z(\varphi)$. In particular, we are interested in
the charge-balanced stimuli of length $\delta$ that fulfill the condition

$\int_{0}^{\delta}{\cal P}(t)\,\mathrm{d}t=0\;,\qquad{\cal P}(t)=0\text{ for
}t\notin[0,\delta]\;.$ (2)

For neuroscience applications, such stimuli are required to avoid charge
accumulation in live tissue. We denote the response to non-Dirac stimuli as
empirical PRC $Z_{\cal P}(\varphi)$. Thus, the first problem is to relate
$Z(\varphi)$ and $Z_{\cal P}(\varphi)$ for given ${\cal P}(t)$. In the
following, we infer the oscillator’s response by perturbing it with the pulse
train

$p(t)=\sum_{k}{\cal P}(t-t_{k})\;;$

we specify the choice of stimulation times $t_{k}$ in test examples below. The
phase shift caused by a stimulus certainly also depends on the stimulus’
amplitude and hence has to be normalized. For unimodal stimuli, it is natural
to normalize by the integral $f=\int_{0}^{\delta}{\cal P}\,\mathrm{d}t$,
commonly referred to as the stimulus’ ”action” 111We highlight that in the
context of electrical stimulation, e.g., in neuroscience applications, $p(t)$
is a current or a voltage and the introduced action has the physical meaning
of electrical charge.. For bimodal charge-balanced pulses the normalization is
ambiguous; we choose $f=\frac{1}{2}\int_{0}^{\delta}|{\cal P}|\,\mathrm{d}t$.
Thus, $Z_{\cal P}=\Delta\varphi/f$.

The success of the standard technique is due to the presence of slow and fast
motion in neuronal oscillators. A spike corresponds to the epoch of fast
motion, where the isochron density in the phase space is low. Thus, spike
detection via threshold-crossing is weakly dependent on the threshold and,
hence, the precision of the phase estimation is high. Consider now a system
without slow-fast motion that exhibits nearly sinusoidal oscillation. Although
there are no well-defined markers, the points of a constant phase value can
also be determined by threshold-crossing events. However, the error in phase
estimation should be higher than in the case of neuron-like oscillators,
especially in the presence of noise. Therefore, the natural idea is to avoid
determining the marker events and estimate the instantaneous phase. The first
step in this direction has been done by Holt et al. Holt and Netoff (2014);
Holt  _et al._ (2016). They suggested fitting sine waves to the several
oscillation periods before and after each stimulus. The phase shift caused by
the stimulus is then easily obtained from these two sines. Duchet et al.
Duchet  _et al._ (2020) approached this problem using the Hilbert Transform:
they estimated the instantaneous phase before and after the stimulus and in
this way computed the PRC. 222Instantaneous phase and amplitude of a real-
valued signal $s(t)$ are the argument and absolute value of the complex-valued
analytical signal $s(t)+\mathrm{i}s_{H}(t)$, where $s_{H}$ is the Hilbert
Transform of $s$, see, e.g., pikovsky2001; King-09; Feldman (2011).
Unfortunately, neither of the techniques published in Holt and Netoff (2014);
Holt  _et al._ (2016); Duchet  _et al._ (2020) has been tested on an
oscillatory model with known PRC; therefore, the performance of these
techniques is unclear. Below, we test these techniques and propose a
different, more precise approach.

####  I.1.2 Estimating amplitude response in experiments

The problem of amplitude response estimation appears, e.g., in the context of
DBS, where the goal is to modulate a brain rhythm by weak pulses. Such
modulation is possible only if the amplitude’s perturbation decays slowly so
that the effects of several pulses accumulate. For a dynamical model, it means
that the brain oscillator’s limit cycle is weakly stable, i.e., its largest
Floquet multiplier is close to one (Floquet exponent close to 0). A popular
model of brain rhythm generation is a large neuronal population exhibiting
collective mode oscillation due to the synchronization of the population’s
elements. If the amplitude of the collective mode is not too large, i.e., the
system is close to the synchronization transition point, then the limit cycle
corresponding to the collective oscillation is weakly stable even if the
cycles of individual neurons are strongly stable. This consideration motivates
our analysis of the amplitude response using the models with a weakly stable
cycle.

The first problem with the amplitude response is its definition. While the
phase of a limit-cycle oscillator is unambiguously (up to an additive
constant) defined via isochrons, there is no universal definition for the
amplitude variable. A pragmatic approach Duchet  _et al._ (2020) is to
estimate the amplitude of an observed signal using the Hilbert Transform.
Comparing the amplitude before and after the stimulus, for different phases of
stimulation, Duchet et al. Duchet  _et al._ (2020) in this way introduced the
amplitude response curve, ARC.

Another approach to characterization of the amplitude variations exploits the
notion of isostables Wilson and Moehlis (2016b); Wilson and Ermentrout
(2018a); Wilson (2019). In this framework, the Winfree Eq. (1) is complemented
by an equation for the isostable variable $\psi$ that can be interpreted as a
deviation from the limit cycle:

$\dot{\psi}=\kappa\psi+I(\varphi)p(t)\;.$ (3)

The real-valued parameter $\kappa$ (Floquet exponent) determines the stability
of the cycle and is negative for the stable one. The function $I(\varphi)$
quantifies the phase dependence of the stimulation’s effect. We denote
$I(\varphi)$ as the isostable response curve, IRC. For an introduction of the
phase-isostable representation, see Appendix A. We remind that Eqs. (1,3)
represent first-order approximation in the perturbation’s strength.

This paper elaborates on the phase and amplitude response inference employing
test stimulation. To evaluate the performance of different techniques, we test
them on systems for which this response is known. First, we analyze the
existing signal analysis techniques and propose a new, more efficient one.
Next, we present the method for inference of the phase-isostable Eqs. (1,3)
for the first time.

##  II The test models

###  II.1 Two-dimensional test models

Our first test model is the Stuart-Landau (SL) oscillator. The reasons for
using this system are threefold. First, the parameters of this system
explicitly govern the limit cycle’s stability and shape of isochrons. Second,
this system’s phase and amplitude response can be obtained analytically.
Third, this system represents the normal form of the Hopf bifurcation and,
therefore, serves as the model equation for the collective mode of an
oscillator population close to the synchronization transition point. The
system’s equations are:

$\displaystyle\dot{x}$ $\displaystyle=\mu x-\eta y-(x^{2}+y^{2})(x-\alpha
y)+\sigma\xi_{x}(t)+\cos\beta\cdot p(t)\;,$ (4) $\displaystyle\dot{y}$
$\displaystyle=\mu y+\eta x-(x^{2}+y^{2})(y+\alpha
x)+\sigma\xi_{y}(t)+\sin\beta\cdot p(t)\;,$

where the parameter $\mu$ governs the limit cycle’s stability as it is
directly related to its Floquet exponent by $\kappa=-2\mu$ and $\alpha$ is the
non-isochronicity parameter. The frequency of the limit-cycle oscillation is
given by $\omega:=\eta-\alpha\mu$. $\sigma$ is the noise strength, and
$\xi_{x,y}$ are two realizations of the Gaussian $\delta$-correlated white
noise with zero average and unit variance. $p(t)$ is the external perturbation
specified in the test examples below, and parameter $\beta$ determines the
perturbation’s direction with respect to $x$ and $y$ variables. The solution
of the noise-free and unperturbed SL system is a sine with the amplitude
$\sqrt{\mu}$. For definiteness, throughout this paper, we assume that the
perturbation is state-independent in the chosen Cartesian coordinates (4).
Then, the SL system (4) has PRC (see, e.g., Rosenblum and Pikovsky (2019))

$Z(\varphi)=-\frac{1}{\sqrt{\mu}}\left(\sin(\varphi-\beta)+\alpha\cos(\varphi-\beta)\right)\;$
(5)

and IRC

$I(\varphi)=\frac{2}{\sqrt{\mu}}\cos(\varphi-\beta)\;.$ (6)

In the following, we need the empirical amplitude response
$A(\varphi)=R_{a}/R_{b}$, where $R_{b,a}$ are the values of the amplitude
variable $R$ immediately before and after the applied pulse. For a narrow
monophasic pulse this expression reads:

$A(\varphi)=1+\frac{f\cos(\varphi-\beta)}{\sqrt{\mu}}\;,$ (7)

where $f$ is the pulse’s action. For derivation of IRC and ARC for the SL
model and their interrelation, see Appendix A.1.

We modify the SL oscillator to obtain the second test model that provides a
non-harmonic solution but is still analytically tractable. We denote this
system as the modified SL oscillator (mSL). The system’s equations that we
derive in Appendix A.2 include the frequency of the limit cycle’s oscillation
$\omega$, the Floquet exponent $\kappa$, and non-isochronicity $\alpha$ as
parameters. With an account of the perturbations, these equations read:

$\displaystyle\dot{x}$
$\displaystyle=\omega\left[xC(x,y)-y\right]+\frac{\kappa}{2}\left[D(x,y)-1\right]\left[x+\alpha\left(xC(x,y)-y\right)\right]+\sigma\xi_{x}(t)+\cos\beta\cdot
p(t)\;,$ (8) $\displaystyle\dot{y}$
$\displaystyle=\omega\left[yC(x,y)+x\right]+\frac{\kappa}{2}\left[D(x,y)-1\right]\left[y+\alpha\left(yC(x,y)+x\right)\right]+\sigma\xi_{y}(t)+\sin\beta\cdot
p(t)\;,$

where the functions $C$ and $D$ are:

$\displaystyle C(x,y)$ $\displaystyle=-\frac{2xy}{(r+2)x^{2}+ry^{2}}\;,$ (9)
$\displaystyle D(x,y)$
$\displaystyle=\frac{(x^{2}+y^{2})^{2}}{(r+2)x^{2}+ry^{2}}\;,$ (10)

and $r$ is a positive parameter determining the shape of the limit cycle in
state space. System (8) has IRC

$\displaystyle
I(\varphi)=2\frac{(r+1)\cos(\varphi-\beta)+\cos(3\varphi-\beta)}{(r+2\cos^{2}(\varphi))^{\frac{3}{2}}}$
(11)

and PRC

$\displaystyle
Z(\varphi)=-\frac{\sin(\varphi-\beta)}{\sqrt{r+2\cos^{2}(\varphi)}}-\alpha\frac{(r+1)\cos(\varphi-\beta)+\cos(3\varphi-\beta)}{(r+2\cos^{2}(\varphi))^{\frac{3}{2}}}\;.$
(12)

The system’s dynamics are illustrated in Fig. 1.

Figure 1:  Panel (a) shows the phase space of the modified SL model (11) with
the stable limit cycle as the solid black line, the isochrons as solid brown
lines and the isostables as dashed grey lines. The blue arrow indicates the
direction in which stimulation enters the system, i.e., here $\beta=0$. An
epoch of the time evolution of the $x$\- and $y$-component on the limit cycle
is illustrated in panel (b) as dashed and solid line, respectively. In panel
(c), the PRC $Z$ (greenish, triangles down) and IRC $I$ (orange, triangles up)
for forcing in $x$-direction are shown. Parameter values are $\omega=1$,
$\kappa=-0.1$, $\alpha=0.0$, $r=0.75$.

###  II.2 High-dimensional test model

Our third test system is the simplest model of macroscopic neuronal
oscillation. For this goal, we take a system of $N$ globally coupled
Bonhoeffer–van der Pol oscillators:

$\displaystyle\dot{x}_{k}$
$\displaystyle=x_{k}-x_{k}^{3}/3-y_{k}+I_{k}+\varepsilon X+p(t)\;,$ (13)
$\displaystyle\dot{y}_{k}$ $\displaystyle=0.1(x_{k}-0.8y_{k}+0.7)\;,$

where $k$ is the oscillator index, $k=1,\ldots,N$, and the term $\varepsilon
X$, where $X=N^{-1}\sum_{k}x_{k}$, describes the mean-field coupling. In the
following, we take $N=1000$. The oscillators’ frequencies are determined by
the parameter $I_{k}$ that is Gaussian-distributed with the mean $0.6$ and
standard deviation $0.1$. The coefficient $\varepsilon$ explicitly describes
the interaction between the elements of the ensemble. We choose
$\varepsilon=0.023$; for this parameter’s value the system exhibits collective
chaos, so that the mean-field oscillation is amplitude-modulated. Figure 2
presents the two-dimensional trajectory in coordinates $X$ and
$Y=N^{-1}\sum_{k}y_{k}$ in (a) and $X(t)$ in (b). With this system we imitate
natural variability of real-world signals.

Figure 2: Collective oscillation on an ensemble of globally coupled units Eqs.
(13). (a) The trajectory in mean-field coordinates (see text). (b) Time
dependence of the mean-field oscillation exhibits amplitude modulation,
typical, e.g., for band-pass filtered brain activity, cf. Fig. 1 in Ref.
Duchet  _et al._ (2020).

##  III Methods and their performance

We formulate the inference problem as follows. Suppose we perturb the system
by applying some known stimulation $p(t)$ and measure the system’s scalar
output $s(t)$. We process $s(t)$ to infer the PRC $Z(\varphi)$ and IRC
$I(\varphi)$ or ARC $A(\varphi)$ as specified below. Generally, we are free to
construct $p(t)$, e.g., as a sequence of pulses. However, some constraints
exist in specific settings, e.g., the pulses must be charge-balanced in
neuroscience applications. In this section, we first introduce the error
measures. Then we suggest a simple approach relating the phase response to
stimuli of arbitrary shape back to the infinitesimal PRC. Next we review the
standard PRC inference technique, followed by testing the performance of other
techniques in use. And finally, we describe and test our approach to the
problem.

###  III.1 Error measures

We define two error measures with which we quantify the goodness of inference.
If the true response curve is known from the theory, we evaluate the goodness
of the inference by computing the normalized $L_{2}$ distance between the
inferred and true curves. We denote these errors with $L_{Z,A,I}$, where the
subscript indicates which curves we are comparing (either PRC, ARC, or IRC).
Thus, for error of, e.g., the PRC recovery, we obtain:

$L_{Z}=\sqrt{{\int\limits_{0}^{2\pi}(Z(\varphi)-Z^{\text{rec}}(\varphi))^{2}\,\mathrm{d}\varphi}\;\big{/}\;{\int\limits_{0}^{2\pi}(Z-\bar{Z})^{2}(\varphi)\,\mathrm{d}\varphi}}\;.$
(14)

where $\bar{Z}$ is the average value of the PRC
$\bar{Z}=(2\pi)^{-1}\int\limits_{0}^{2\pi}Z(\varphi)\,\mathrm{d}\varphi$, and
$Z^{\text{rec}}$ the inferred (recovered) PRC. The error value is
$L_{Z,A,I}=0$ if the two curves coincide, it is of order one if the two curves
are of the same order of magnitude but not similar, and can also take higher
values if one curve is significantly larger on average.

The other measure quantifies how well the inferred curve represents the
dynamical model. Since the dynamics of phase and amplitude are distinct, see
Eqs. (1) and (3), respectively, the definition of the corresponding error
measures differs as well. We first introduce the error measure for the phase
response and later in Sec III.5.2 when introducing the method we also specify
an analog for the isostable variable error. Note that the inferred phase
response curve $Z(\varphi)$ in conjunction with the natural frequency $\omega$
can be used in Eq. (1) to reproduce $\varphi(t)$ for a given realization of
stimulation $p(t)$ by means of numerical integration. In the first step, we
threshold the observed signal to determine time events $\tau_{i}$ that
correspond to the same phase; for details, see Appendix C. Without loss of
generality, we set this phase to zero. Next, starting at $\tau_{i}$ we
reproduce the phase evolution up to the time $\tau_{i+1}$. In the ideal noise-
free case, if the stimulation is weak and the inferred curve is exact, the
reproduced phase $\varphi(\tau_{i+1})=\Phi_{i}$ equals $2\pi$ (or,
equivalently, zero since we consider the wrapped phase). In practice, this
reproduced phase $\Phi_{i}$ deviates from $2\pi$, and it is precisely this
deviation that we use to quantify the inference’s quality. We define the error
of PRC reconstruction $E_{Z}$ as the standard deviation of $\Phi_{i}$ from
$2\pi$:

$E_{Z}=\langle(\Phi_{i}-2\pi)^{2}\rangle^{1/2}\;,$ (15)

where $\langle\cdot\rangle$ denotes averaging over the index $i$. The value of
$E_{Z}$ should be compared to a measure of the signal’s irregularity. We
evaluate the latter by a measure that is proportional to the standard
deviation of the inter-event intervals $T_{i}=\tau_{i+1}-\tau_{i}$
(instantaneous periods):

$E_{Z0}=\frac{2\pi}{\langle T_{i}\rangle}\langle(T_{i}-\langle
T_{i}\rangle)^{2}\rangle^{1/2}\;.$ (16)

A natural measure of the inference’s quality is the ratio ${E_{Z}}/{E_{Z0}}$
since it roughly describes how much of the signal’s irregularity is explained
by the inferred phase model. If one considers a PRC that is identically zero,
this ratio equals one.

We emphasize an essential advantage of the error measure ${E_{Z}}/{E_{Z0}}$.
The computation of the $L_{2}$-based measure Eq. (15) requires knowledge of
the ground truth and, therefore, helps only in testing the techniques on model
data from limit-cycle oscillators. In contrast, the error measure
${E_{Z}}/{E_{Z0}}$ is obtained solely from data and, hence, is a helpful tool
for experimental data analysis.

###  III.2 Relation between empirical $Z_{\mathcal{P}}$ and infinitesimal $Z$
PRCs

From the infinitesimal PRC $Z(\varphi)$ one can compute the empirical PRC
$Z_{\mathcal{P}}$ as a response to an arbitrary stimulus $\mathcal{P}$. One
simply has to evaluate the effective phase shift for the chosen stimulation
$\mathcal{P}$ at every phase by integrating the Winfree phase equation (1).
For example, let us focus on evaluating the phase response for stimuli
$\mathcal{P}$ at a particular phase $\varphi^{*}$. One first needs to compute
the phase in time for the duration of the pulse:

$\varphi(t)=\varphi^{*}+\int\limits_{0}^{t}(\omega+Z(\varphi(t^{\prime}))\mathcal{P}(t^{\prime}))\,\mathrm{d}t^{\prime}\;.$
(17)

Then the empirical PRC at the chosen phase corresponds to the phase shift that
is induced by $\mathcal{P}$ over the duration of the stimulation:

$Z_{\mathcal{P}}(\varphi^{*})=\frac{1}{f}\int\limits_{0}^{\delta}Z(\varphi(t))\mathcal{P}(t)\,\mathrm{d}t\;,$
(18)

where $f$ is the action of the pulse.

The inverse problem of how to determine the infinitesimal response
$Z(\varphi)$ when given the empirical one $Z_{\mathcal{P}}$ and stimulation
shape $\mathcal{P}$, is in general much harder. However, for pulses with small
action $f$, the empirical response $Z_{\mathcal{P}}(\varphi)$ is well
approximated as the convolution of the infinitesimal curve $Z(\varphi)$ and
the pulse shape $\mathcal{P}(t)$. In Fourier space, convolution corresponds to
a multiplication, meaning the inverse operation is represented by division.
One can, therefore, easily deconvolve an empirical response to the
infinitesimal by dividing the Fourier representations of $Z_{\mathcal{P}}$ and
$\mathcal{P}$; we present the technique in more detail in Appendix B. There we
also present an algorithm that does not assume the smallness of $f$. The
transformation $Z_{\cal P}\to Z$ is straightforward and can be combined with
any technique that infers an empirical response.

###  III.3 Standard approach (phase response only)

We start with a brief review of a practical implementation of the standard
approach for the case where well-defined marker events in the observed signal
are absent. To trace the variation of the period, we determine states with
equal phases employing threshold-crossing. 333Alternatively, one can detect
the signal’s maxima; this approach corresponds to threshold-crossing for the
derivative. For example, we apply a stimulus and look for the events with zero
crossings from below. However, if the limit cycle is not strongly stable,
observing such an event immediately after the stimulus’s application at
$t_{s}$ does not suffice. Indeed, we have to wait until the system returns to
the limit-cycle. Thus, we take the instant $\tau_{n}$ of the $n^{\text{th}}$
threshold-crossing event and compute

$Z_{\cal P}(\varphi)=2\pi\frac{\tau_{n}-\tau_{0}-nT_{0}}{fT_{0}}\;,$ (19)

where $t_{0}$ is the threshold-crossing event preceding the stimulus, and the
phase of the stimulus’ application is $\varphi=2\pi(t_{s}-\tau_{0})/T_{0}$
where $T_{0}$ is the natural period. The choice of $n$ depends on the
stability of the limit cycle, quantified with the Floquet exponent $\kappa$;
as a rule of thumb, we suggest $n\kappa T_{0}\gg 1$. However, in practice,
$\kappa$ is unknown, and $n$ shall be chosen by trial and error. To obtain the
function $Z_{\cal P}(\varphi)$, we repeat the perturbation $\mathcal{P}$ for
different $\varphi$. To this end, we either choose random intervals between
stimuli or stimulate periodically with the period incommensurate with $T_{0}$.
In any case, the stimulation interval shall be large enough to ensure
observation of $n$ threshold-crossing events after each applied pulse.

We illustrate the inference via the standard technique in Fig. 3 where we
consider a charge-balanced stimulation and infer the response for the two test
models presented earlier, see Eqs. (4,8). For the SL system (4) the parameter
values are: $\omega=1$, $\kappa=-0.1$, $\alpha=-0.3$ and for the mSL system
(8) they are $\omega=1$, $\kappa=-0.1$, $\alpha=0$, $r=0.75$. There is no
noise $\sigma=0$. For the observable we take the $x$ variable and $n=3$
threshold-crossing events after stimulation were considered for the phase
shift determination. The charge-balanced pulse used has a short interval of
positive stimulation (duration 0.2), followed by an interval without
stimulation (duration 0.4) and then a longer interval of negative stimulation
(duration 1.0). The amplitude of stimulation depends on the oscillator, the
action $f=\frac{1}{2}\int|\mathcal{P}|\,\mathrm{d}t$ is 0.01 for SL and 0.07
for mSL, see Appendix section D for details on perturbation generation. We
show the following curves: (i) theoretical PRC $Z(\varphi)$ according to Eqs.
(5) and (12), (ii) theoretical response to charge-balanced pulses $Z_{\cal
P}$, (iii) empirical PRC $Z_{\cal P}$ obtained via Eq. (19) (“direct
inference”) 444The empirical curve is an $8^{\text{th}}$-order Fourier fit of
experimental points., and (iv) the result of the deconvolution of the
empirical PRC $Z_{\cal P}$ as described in Appendix B (“deconvolved”).

We emphasize that there is no analog to the standard technique for inferring
the amplitude response.

Figure 3: Inference of the phase response by the standard technique, for
charge-balanced stimulation and for two test oscillators, (a) Stuart-Landau
and (b) its modification Eq. (8). Orange (full) and gray (dashed) bold
transparent curves depict the theoretical curves $Z(\varphi)$ and $Z_{\cal
P}(\varphi)$. Red (full) and blue (dashed) thin curves show the results of the
direct inference and of the deconvolution (see text for explanations). The
purple line demonstrates the pulse shape as a function of
$t=\frac{\varphi}{2\pi}T_{0}$ (as it would appear next to an unperturbed
signal). The pulse is vertically scaled to fit the plot. We see a very good
correspondence between the theoretical and inferred curves.

###  III.4 Inferring response measuring instantaneous phase and amplitude

The idea of the approach is straightforward. Suppose we perturb the oscillator
by a sequence of pulses at instants $t_{k}$ and measure the system’s output
$s(t)$. Computing the instantaneous phase and amplitude of $s(t)$ before and
after each pulse, we obtain the phase shift and amplitude variation as
functions of the stimulation’s phase. To obtain a reasonable estimation of
these functions, we have to apply the pulses at different phases – we ensure
this by either using stimulation with a period incommensurate with that of the
unperturbed oscillator or by choosing random $t_{k}$.

Let the instantaneous phase and amplitude immediately before and after a
finite-width stimulus be $\varphi_{a,b,}$ and $a_{a,b}$, respectively, where
$a$ and $b$ stand for “after” and “before”. From these quantities, we obtain
the empirical PRC $Z_{\cal
P}(\varphi_{b})=(\varphi_{a}-\varphi_{b}-\omega\delta)/f$, where $\delta$ is
the pulse’s width, $\omega$ is the frequency of the unperturbed oscillation,
and $f$ is the normalization factor. Next, following Duchet et al. Duchet  _et
al._ (2020), we introduce the empirical amplitude response curve (ARC) as
$A(\varphi)=a_{a}/a_{b}$, to be distinguished from IRC. Below, we explore the
performance of three different techniques.

####  III.4.1 Fitting the signal by a sine (phase response only)

Holt et al. Holt and Netoff (2014); Holt  _et al._ (2016) suggested fitting a
harmonic to several periods of the signal $s(t)$ before the stimulus and
another harmonic to an interval of the same length $T_{\text{fit}}$ after the
stimulus. Then, they exploited the Fourier Transform of both fitted functions
to find the phase shift evoked by the stimulus. Obviously, this technique does
not provide information on the amplitude response. Indeed, the perturbation in
the amplitude normally decays rather quickly and cannot be captured by a fit
over several oscillation periods. We illustrate the performance of this
technique on the test model (4) in Fig. 4, computing the phase shift in the
time domain. Suppose the pulse of length $\delta$ occurs at $t_{s}$. Let the
fitted functions be $s(t)\approx\
a_{1}\cos[\omega(t-(t_{s}-T_{\text{fit}}))+\eta_{1}]$ for
$t_{s}-T_{\text{fit}}<t<t_{s}$ and $s(t)\approx\
a_{2}\cos[\omega(t-(t_{s}-T_{\text{fit}}))+\eta_{2}]$ for
$t_{s}+\delta<t<t_{s}+\delta+T_{\text{fit}}$. The phase of the second cosine
at $t=t_{s}+\delta$ is $\varphi=\omega(T_{\text{fit}}+\delta)+\eta_{2}$. If
there were no stimulus, the phase of the first cosine at this point would be
$\varphi=\omega(\delta+T_{\text{fit}})+\eta_{1}$. Hence, the phase shift is
$\Delta\varphi=\eta_{2}-\eta_{1}$. 555Note that the oscillation frequency
$\omega$ shall be either determined by fit or set to be equal to the central
frequency of the bandpass filter, see Holt and Netoff (2014); Holt  _et al._
(2016).

Naturally, this technique applies to the signals close to a harmonic
oscillation. Therefore, we test it on the SL oscillator only. We use the same
parameters as used to illustrate the standard technique, see Fig. 3. For
simplicity, we use monophasic rectangular pulses with amplitude $0.1$ and
width $\delta=0.03$, therefore having action $f=0.003$. Such a stimulus is a
good approximation for the delta-pulse, hence, we can expect a good
correspondence of the inferred PRC with the theoretical one. However, the
results summarized in Fig. 4 demonstrate a limited precision of the approach.
As another drawback of the technique, we mention that, since fitting requires
several periods before and after the stimulus, the stimulation pulses shall be
relatively rare and, hence, the time interval required for PRC estimation is
rather long.

Figure 4: PRC obtained via fitting a cosine prior and after the stimulus, for
different lengths of the fitting interval (the interval is given in periods of
the unperturbed system). The solid line depicts the theoretical PRC
$Z(\varphi)$, see Eq. (5).

####  III.4.2 Estimating phase and amplitude using the Hilbert Transform

A typical way to obtain the narrow-band signal’s instantaneous phase and
amplitude is to use the Hilbert Transform (HT). So, Duchet et al. Duchet  _et
al._ (2020) exploited HT to compute the phase and amplitude response of tremor
oscillation to deep brain stimulation of the thalamus. They took the
instantaneous phase and amplitude immediately before and after the stimulus
for this computation 666In their experiments, Duchet et al. Duchet  _et al._
(2020) used a complex stimulus consisting of 25 bursts of high-frequency
pulses, of the total length of 5 s..

However, being a nonlocal operation, the HT is not suitable for measuring
response to a pulse, e.g., the HT shows the systems’ reaction even before the
stimulus begins. We illustrate this property of the HT by perturbing the SL
oscillator with a rectangular pulse. For a better visibility, we choose a wide
pulse (amplitude 0.1, $\delta=0.6$). The parameters of the SL system are the
same as before. Figure 5 illustrates the results, see also Fig. 3.5 in Feldman
(2011). We see that the instantaneous Hilbert-based amplitude immediately
before the pulse and immediately after it strongly deviates from the actual
value. This deviation is of the same order of magnitude as the amplitude’s
change due to the pulse. We observe similar behavior for the instantaneous
phases.

Figure 5: HT-based amplitude and phase of the perturbed SL oscillator. (a)
Black dashed-dotted curve shows the time dependence of the amplitude variable
$R=\sqrt{x^{2}+y^{2}}$ of the SL system, and the bold red curve shows the
instantaneous amplitude $a_{H}$ computed via the Hilbert Transform. The thin
blue curve depicts the observation $x(t)$. (b) Here, the black dashed-dotted
curve shows the evolution of the oscillator’s phase $\varphi$, while the red
curve illustrates the HT-based phase $\varphi_{H}$ (for both phases, we
subtract the average growth with the frequency $\omega$). The vertical lines
in all panels indicate the beginning and end of the rectangular pulse. In the
absence of stimulation, the HT amplitude $a_{H}(t)$ reproduces $R(t)$, but the
discrepancy is essential before and after the pulse. For the phases, the
effect is even more pronounced.

In summary, we shall interpret the HT-based phase and amplitude with caution
in this context.

We suggest improving the performance of the HT-based inference in the
following way. We neglect the instantaneous phases in the intervals
$(t_{s}-\delta_{\text{off}},t_{s})$ and
$(t_{s}+\delta,t_{s}+\delta+\delta_{\text{off}})$, where $t_{s}$ is the
instant of the stimulus’s onset and $\delta_{\text{off}}$ is the offset time.
Next, we obtain the phase $\varphi_{b}$ by extrapolating the linear fit of
$\varphi(t)$ over the interval
$(t_{s}-\delta_{\text{off}}-\delta_{\text{fit}},t_{s}-\delta_{\text{off}})$ to
the instant $t_{s}$; the parameter $\delta_{\text{fit}}$ is the length of the
fitting interval. Note that fitting requires the unwrapped phase. Similarly,
extrapolating the fit over the interval
$(t_{s}+\delta+\delta_{\text{off}},t_{s}+\delta+\delta_{\text{off}}+\delta_{\text{fit}})$
to $t_{s}+\delta$ yields $\varphi_{a}$. Moreover, linear fit provides the
frequency $\omega$. Hence, we compute $Z_{\cal P}(\varphi_{b})$. The amplitude
response (ARC) we compute as
$A(\varphi_{b})=a_{H}(t_{s}+\delta+\delta_{\text{off}})/a_{H}(t_{s}-\delta_{\text{off}})$.
We emphasize that the algorithm has two parameters, $\delta_{\text{off}}$ and
$\delta_{\text{fit}}$.

We demonstrate the performance of the described algorithm perturbing the SL
system (4) by rectangular pulses, taking the inter-pulse interval
incommensurate with the natural period $T_{0}$ – in this way, we ensure that
stimulation occurs at different phases. We set the observable $s(t)=x(t)$ and
first choose $\beta=0$. Other parameters are the same as in the illustration
of the standard and sine-fitting techniques. We illustrate the results in Fig.
6a,b. Here, we show the theoretical curves for $Z(\varphi)$, $A(\varphi)$
according to Eqs. (6,7) along with inferred characteristics $Z_{H}(\varphi)$,
$A_{H}(\varphi)$. To compute the latter, we choose
$\delta_{\text{fit}}=50\delta$ and $\delta_{\text{off}}=2\delta$ and perform
an 8th-order Fourier fit of obtained points. Numerical tests show that
although the results are not very sensitive to the choice of
$\delta_{\text{fit}}$, the choice of $\delta_{\text{off}}$ is crucial. We
demonstrate this in Fig. 6c,d by showing the dependence of the error of
inference on $\delta_{\text{off}}$. We compute this error according to Eq.
(14). We see, that a proper choice of the offset $\delta_{\text{off}}$
essentially affects the inference; the reasonable results in Fig. 6a,b are due
to the optimal value $\delta_{\text{off}}/\delta=2$. Unfortunately, we do not
see any practical way to choose $\delta_{\text{off}}$ when the true curve is
unknown 777For PRC estimation, a possible approach would be computing the
error measure $E_{Z}$ according to Eq. (15) as a function of
$\delta_{\text{off}}$ and searching for the minimum. However, this approach
does not apply to ARC estimation since Hilbert technique does not yield the
Floquet exponent. While we can use the same $\delta_{\text{off}}$ value to
infer both curves, there is no guarantee that the optimal value is the same. .
Moreover, the results strongly depend on how the stimulation enters the
system’s equations. To demonstrate this, we plot in Fig. 6c,d the
corresponding curves for $\beta=\pi/4$ and $\beta=\pi/2$. We see that
generally, the error of inference is not small. Summarizing this example, we
say that the HT-based phase and amplitude response inference generally yields
imprecise results. Though one can use this approach to obtain some empirical
measures of the response, see Duchet  _et al._ (2020), these inference’s
results may be only loosely related to the theory.

Figure 6: HT-based phase (a) and amplitude (b) responses’ inference for the SL
oscillator. Dashed black and solid red curves show the theoretical and
inferred curves, respectively, for $\beta=0$. Parameter $\beta$ describes how
the stimulation enters Eqs. (4). Panels (c,d) show the inference errors as a
function of the offset parameter $\delta_{\text{off}}$ for $\beta=0$ (black
circles), $\beta=\pi/4$ (red crosses), and $\beta=\pi/2$ (blue stars). An
optimal choice of $\delta_{\text{off}}$ ensures successful inference shown in
(a,b). However, the optimization requires knowledge of the investigated
system, making the HT-based technique of limited use in practical
applications.

Additionally, we mention that HT serves merely as a signal embedding
technique, yielding the angle variable or protophase; see Kralemann  _et al._
(2007, 2008) for a discussion. A protophase generally depends on the
embedding; so, e.g., the Hilbert-based protophase does not necessarily
coincide with the angle variable in the $x,y$ plane. Although protophases and
the true (asymptotic) phase provide the same average frequencies, they
generally differ microscopically, i.e., on a time scale smaller than the
period. This difference is due to the non-uniform rotational velocity of the
protophase.

####  III.4.3 Estimating phase and amplitude using a virtual auxiliary
oscillator

Finally, we adapted a technique Rosenblum  _et al._ (2021) for real-time
estimation of phase and amplitude for our purpose. This technique exploits two
virtual linear oscillators - one for phase and one for amplitude determination
- to yield a causal estimation. Namely, we use the signal $s(t)$ as an input
to a damped oscillator $\ddot{x}+\alpha_{a,\varphi}\dot{x}+\eta^{2}x=s(t).$ We
choose the oscillators’ frequencies $\eta$ to be much larger than the
characteristic frequency $\nu$ of $s(t)$ so that the systems are far from the
resonance and the response weakly depends on $\nu$. Next, we take the damping
parameters $\alpha_{\varphi}$ and $\alpha_{a}$ for the phase and amplitude
measurement to ensure a simple relation between the phase
$\varphi(t)=\arctan(-\dot{x}/\nu x)$ and amplitude
$\sqrt{x^{2}+(\dot{x}/\nu)^{2}}$ of the forced oscillation and those of the
investigated signal. The implementation is simple and boils down to numerical
integration of the linear oscillator’s equation driven by a signal given at
discrete time points. For details, see Appendix E.

We expect this causal approach to yield a precise estimation of the amplitude
or phase before the stimulus. Immediately after the stimulus, the estimation
is poor due to transients. Indeed, the approach implies that the oscillation
with the oscillator’s frequency $\eta$ decays, and only the oscillation with
the frequency of the input $s(t)$ remains. We suggest the following solution
to this problem. We compute the amplitude twice, first for the original signal
and then for the time-inverted one. Namely, we flip the time series $s_{k}$ so
that the time “runs” backward, i.e., from the $N$-point time series $s_{k}$,
$k=1,2,\ldots,N$, we construct the flipped series $\hat{s}_{1}=s_{N}$,
$\hat{s}_{2}=s_{N-1}$, …, $\hat{s}_{N}=s_{1}$. Next, we compute the amplitude
of $\hat{s}_{k}$ via the same algorithm and flip it in time to obtain the
flipped amplitude denoted by $a_{f}(t)$. $a_{f}(t)$ is a precise estimation
for the time interval after the stimulus while the transients corrupt the
estimation before the stimulus. However, since the phase changes its sign with
the time inversion, we must reflect the obtained phase,
$\varphi_{f}(t)\to-\varphi_{f}(t)$. We illustrate this algorithm in Fig. 7;
parameters of the virtual oscillator are $\eta=5$, $\alpha_{a}=6$,
$\alpha_{\varphi}=0.2$.

Figure 7: Amplitude and phase of the perturbed SL oscillator, obtained via the
linear non-resonant oscillator technique. (a) Black dashed-dotted curve shows
the time dependence of the amplitude variable $R=\sqrt{x^{2}+y^{2}}$ of the SL
system, and the red dashed and solid blue curves show the instantaneous
amplitude $a_{\text{nro}}$, where the subscript stays for non-resonant
oscillator. Note that dashed and solid lines correspond to time series $s_{k}$
and $\hat{s}_{k}$, respectively. (b) Here, the black dashed-dotted curve shows
the evolution of the oscillator’s phase $\varphi$, while the red dashed and
blue solid curves illustrate the estimated phase $\varphi_{\text{nro}}$,
computed from $s_{k}$ and $\hat{s}_{k}$, respectively (for both phases, we
subtract the average growth with the frequency $\omega$). The vertical lines
in all panels indicate the beginning and end of the rectangular pulse. We see
that the amplitude, computed in forward time yields a good estimate of the
amplitude immediately before the pulse. Similarly, the amplitude computed in
backward time reasonably traces the amplitude after the pulse. For the phases,
the estimation is less successful (at least for this particular pulse).

A remark is in order. The linear oscillator used to compute the phase has a
smaller damping parameter than the oscillator for the amplitude estimation.
Correspondingly the transients in the phase measurement are essentially longer
than in the case of the amplitude measurement, which makes the phase
measurement less precise. We present the inference’s results for the SL model
in Fig. 8. Additionally, just like any approach relying on estimating small
differences of the angle variable, it generally suffers significant errors
because it does not include information on local isochrons.

Figure 8: Phase (a) and amplitude (b) responses for the SL oscillator,
inferred using auxiliary non-resonant linear oscillators. Dashed black and
solid red curves show the theoretical and inferred curves, respectively. We
show the results for $\beta=\pi/2$, where the HT-based technique fails. We
recall that parameter $\beta$ describes how the stimulation enters Eqs. (4).

###  III.5 Inferring response curves by reconstructing the first order phase –
isostable dynamics: the IPID-1 technique

Here, we describe a method for inferring the first-order phase-isostable
dynamics from observations by fitting the model Eqs. (1,3). We name the method
IPID-1 standing for “Inferring Phase–Isostable Dynamics of order 1”. We carry
out the procedure in two steps. First, we infer the PRC and instantaneous
phase by adapting the algorithm introduced in Ref Cestnik and Rosenblum (2018)
for the case of pulse stimulation. Next, we use the inferred phase to
reconstruct the isostable dynamics. As in the previous sections, we assume
that a scalar signal $s(t)$ and the perturbation $p(t)$ are known.

####  III.5.1 Inference of the phase response

The key step in the inference is determining time instants $\tau_{i}$
corresponding to the same asymptotic phase $\varphi$; these instants must be
extracted from the observed signal $s(t)$. Sometimes the choice of such events
is obvious, e.g., in the case of a spiky signal where spikes indicate the same
phase. In general one considers threshold-crossing events,
$s(\tau_{i})=s_{\text{thr}}$, choosing one crossing per period, e.g., with an
additional condition $\frac{d}{dt}s(\tau_{i})>0$. The choice of the threshold
$s_{\text{thr}}$ is important. The proper thresholding should closely match
the crossing of a local isochron, see Appendix C for details. Since isochrons
make a full rotation over a limit cycle, at least two such thresholds exist
for any scalar signal - we find appropriate thresholds with a direct search,
seeing which threshold yields the best fit in terms of the error (15).

The core of the inference is fitting the Winfree equation (1) integrated over
single periods, determined as the time interval between two threshold-crossing
events, $[\tau_{i},\tau_{i+1}]$:

$2\pi=\omega(\tau_{i+1}-\tau_{i})+\int\limits_{\tau_{i}}^{\tau_{i+1}}Z(\varphi(t))p(t)\,\mathrm{d}t\;.$
(20)

The left-hand side equals $2\pi$ due to the definition of events $\tau_{i}$
having the same phase. On the right-hand side, we approximate the PRC by a
finite Fourier series of order $N_{F}$:
$Z(\varphi)=\sum_{n=0}^{N_{F}}\left[z_{n}^{\cos}\cos(n\varphi)+z_{n}^{\sin}\sin(n\varphi)\right]$.
By interchanging the order of integration and summation, we obtain from Eqs.
(20) a linear system for the unknown Fourier coefficients $z_{n}$ and
frequency $\omega$. Since we have as many equations (20) as there are periods
of the observed signal, for a sufficiently long data set we have more
equations than unknowns. Thus, we can solve the linear system, e.g., by least-
squares minimization. The coefficients $\int\cos(n\varphi)p(t)\,\mathrm{d}t$,
$\int\sin(n\varphi)p(t)\,\mathrm{d}t$, however, cannot be evaluated yet
because the phase $\varphi(t)$ is not known to us a priori. We overcome this
with an iterative procedure; by first approximating the phase to obtain an
approximate solution, and then exploiting this solution to improve the phase
estimate 888Each phase approximation $\varphi^{(m)}$ is used to compute the
integrals that represent coefficients of the linear system:
$\int\cos(n\varphi^{(m)})p(t)\,\mathrm{d}t$,
$\int\sin(n\varphi^{(m)})p(t)\,\mathrm{d}t$. Then when solving the linear
system and obtaining a better PRC approximation $Z^{(m+1)}(\varphi)$, it is
used to recompute the next approximation of the phase $\varphi^{(m+1)}$ by
integrating the Winfree Eq. (1).

Approximating the phase initially as linearly growing between the events,
$\varphi^{(0)}(t)=2\pi\frac{t-\tau_{i}}{\tau_{i+1}-\tau_{i}}$, we obtain the
first-approximation solution of Eq. (20), namely $\omega^{(1)}$ and
$Z^{(1)}(\varphi)$. (The superscripts denote the iteration number.) Then, we
obtain the next approximation of the phase by integrating the Winfree Eq. (1)
between events. In general, the $m^{\text{th}}$-order approximation for the
phase is obtained as:

$\varphi^{(m)}(t)=\int_{\tau_{i}}^{t}\left[\omega^{(m)}+Z^{(m)}(\varphi^{(m)}(t^{\prime}))\right]\,\mathrm{d}t^{\prime}\;.$
(21)

Estimated in this way, the phase at the end of a period generally differs from
$2\pi$: $\Phi_{i}^{(m)}=\lim\limits_{t\uparrow\tau_{i+1}}\varphi^{(m)}(t)\neq
2\pi$. Thus, we additionally re-scale the phase with the factor
$\Phi_{i}^{(m)}$ to ensure $\varphi^{(m)}(\tau_{i+1})=2\pi$. As a result, the
approximations gradually improve through iterations.

Estimated phases at the end of periods $\Phi_{i}$ indicate how well our
inference fits the observations, see error measure (15). For each iteration
$(m)$ we can compute the error and monitor the convergence. For further
details of the technique, we refer to Ref. Cestnik and Rosenblum (2018).

####  III.5.2 Inference of the amplitude response

This section extends our approach to cover the isostable dynamics
reconstruction. We still need the signal $s(t)$ and perturbation $p(t)$, and
since we first infer the PRC we also have the instantaneous phase
$\varphi(t)$.

Like the PRC technique, the isostable inference relies on time events of equal
phase $\tau_{i}$ and, additionally, on estimating the isostable amplitude at
those events, $\psi_{i}\equiv\psi(\tau_{i})$. The time events are
straightforwardly obtained from the instantaneous phase $\varphi(t)$ 999Note
that the events $\tau_{i}$ used for isostable inference generally differ from
those used for the PRC reconstruction, though we keep the same notation.
Practically, to determine $\tau_{i}$ we use linear interpolation of the
monotonically growing function $\varphi(t)$ given in a discrete set of
points., while the isostable variable needs to be estimated from the observed
signal $s(t)$. The IRC function describes a linear response and is operable in
the close vicinity of the limit cycle where isochrons can be approximated by
straight lines. Thus in the first-order approximation, the isostable variable
linearly depends on signal $s(t)$:

$\psi_{i}=c\ (s(\tau_{i})-s_{0})+\mathcal{O}((s(\tau_{i})-s_{0})^{2})\;,$ (22)

where factors $c$ and $s_{0}$ are generally phase dependent. However, since
all crossing events $\tau_{i}$ correspond to the same phase, $c$ and $s_{0}$
are constant in this context. Additionally, the isostable variable is
inherently determined up to a constant factor. Therefore, without loss of
generality, we can set $c\equiv 1$. Note that approximation via Eq. (22) is
also justified for high-dimensional systems given if one Floquet exponent is
significantly smaller in absolute value than the rest, i.e., there is a slow
manifold (see discussion in Appendix A).

We now consider the isostable dynamics (3) integrated over periods determined
by the events of the same phase $[\tau_{i},\tau_{i+1}]$:

$\psi_{i+1}-\psi_{i}=\kappa\int_{\tau_{i}}^{\tau_{i+1}}\psi(t)\,\mathrm{d}t+\int_{\tau_{i}}^{\tau_{i+1}}I(\varphi)p(t)\,\mathrm{d}t\;.$
(23)

This system of equations is similar to Eq. (20) in that it can be approximated
with a linear system by expanding the unknown function in a Fourier series,
$I(\varphi)=\sum_{n=0}^{N_{F}}\left[u_{n}^{\cos}\cos(n\varphi)+u_{n}^{\sin}\sin(n\varphi)\right]$,
and interchanging the order of integration and summation. Furthermore, since
we know the phase as a function of time, the $i$-dependent factors of Fourier
coefficients can be computed directly. Using approximation (22) we evaluate
the left-hand side as $\psi_{i+1}-\psi_{i}=s(\tau_{i+1})-s(\tau_{i})$ because
the still unknown constant $s_{0}$ cancels in the difference. The integral of
the isostable amplitude is challenging since it includes time-dependent
amplitude, which we do not have.

There are two issues with evaluating the amplitude integral
$\int_{\tau_{i}}^{\tau_{i+1}}\psi(t)\,\mathrm{d}t$: (i) we do not know the
constant $s_{0}$ a priori, and (ii) we want to use approximation (22) only in
discrete events $t_{i}$ (otherwise, we would have to consider $c$ and $s_{0}$
as phase-dependent). We tackle the first issue by splitting the integral in
two:

$\kappa\int_{\tau_{i}}^{\tau_{i+1}}\psi(t)\,\mathrm{d}t=-\kappa
s_{0}(\tau_{i+1}-\tau_{i})+\kappa\int_{\tau_{i}}^{\tau_{i+1}}(\psi(t)+s_{0})\,\mathrm{d}t\;.$
(24)

In this representation $\kappa s_{0}$ becomes just another unknown variable
that will be determined while solving the linear system $\eqref{eq:iso_int}$,
which together with determining $\kappa$ allows expressing $s_{0}$ as the
ratio $s_{0}=\frac{\kappa s_{0}}{\kappa}$. The remaining integral involves the
continuous quantity $\psi(t)+s_{0}$, which coincides with our observable
$s(t)$ at events $\tau_{i}$, but generally this quantity differs from the
signal for $t\neq\tau_{i}$: $s(t)\neq\psi(t)+s_{0}$. The property
$s(\tau_{i})=\psi(\tau_{i})+s_{0}$ helps us resolve the second issue of only
using discrete values $s(\tau_{i})$ to estimate the integral of the continuous
amplitude. Just as in the phase response method, we approach this iteratively.
First, we linearly interpolate $\psi(t)+s_{0}$ between known events so that we
can evaluate the last term:
$\int_{\tau_{i}}^{\tau_{i+1}}(\psi(t)+s_{0})\,\mathrm{d}t\approx[s(\tau_{i})+s(\tau_{i+1})]/2$.
Then we solve the linear system approximating Eq. (23) and obtain the first-
approximation solution $s_{0}^{(1)}$, $\kappa^{(1)}$, $I^{(1)}(\varphi)$.
Next, we use this solution to improve the estimate of the continuous amplitude
by integrating the underlying dynamical equation (3). In general, from the
approximate solution of order $m$, we obtain the subsequent amplitude
estimation via integration 101010Since the forcing term is just a function of
time, the isostable dynamics can be integrated via variation of constant, in
which case the time integral can be computed explicitly:
$\psi^{(m)}(t)=(s(\tau_{i})-s_{0}^{(m)})e^{\kappa^{(m)}(t-\tau_{i})}+\int_{\tau_{i}}^{t}I^{(m)}(\varphi(t^{\prime}))p(t^{\prime})e^{\kappa^{(m)}(t-t^{\prime})}\,\mathrm{d}t^{\prime}$:

$\psi^{(m)}(t)=s(\tau_{i})-s_{0}^{(m)}+\int_{\tau_{i}}^{t}[\kappa^{(m)}\psi^{(m)}(t^{\prime})+I^{(m)}(\varphi)p(t^{\prime})]\,\mathrm{d}t^{\prime}\;.$
(25)

Estimated in this way, the amplitude $\psi^{(m)}$ starts the interval
$[\tau_{i},\tau_{i+1}]$ as a direct estimation from the signal:
$s(\tau_{i})-s_{0}^{(m)}$, but further along the interval it can deviate due
to the approximation in integrands $\kappa^{(m)},I^{(m)}$. Thus, at the end of
the interval this quantity does not exactly correspond to
$s(\tau_{i+1})-s_{0}^{(m)}$, which is the value that we take to start the next
interval $[\tau_{i+1},\tau_{i+2}]$. The estimated timeseries of the isostable
amplitude obtained in this way is thus discontinuous at the events $\tau_{i}$.
If the inference were perfect this discontinuity would disappear, which means
we can use the magnitude of the discontinuity to quantify the quality of the
model as follows. Let us denote the amplitude estimated at the end of a period
as $\Psi_{i}^{(m)}=\lim\limits_{t\uparrow\tau_{i+1}}\psi^{(m)}(t)\neq
s(\tau_{i+1})-s_{0}^{(m)}$. We define the error of the fit as the standard
deviation of the difference:

$E_{I}^{(m)}=\langle(\Psi_{i}^{(m)}-(s(\tau_{i+1})-s_{0}^{(m)}))^{2}\rangle^{1/2}\;,$
(26)

Just like the PRC error was compared to the irregularity of inter-event
intervals, see Eq. (16), this error should be compared to the irregularity of
the isostable variable at the events:

$E_{I0}=\langle(s(\tau_{i})-\langle s(\tau_{i}\rangle)^{2}\rangle^{1/2}\;.$
(27)

The entire procedure of inferring the isostable amplitude response is
condensed into step-by-step instructions in the Appendix F.

####  III.5.3 Performance of the IPID-1 phase – isostable reconstruction

We test the performance of the introduced IPID-1 method on the two example
oscillators (4) and (8). We infer the PRC $Z(\varphi)$ and IRC $I(\varphi)$
with the approaches explained in Sections III.5.1 and III.5.2 respectively.
The results can be seen in Fig. 9. Parameters are the same as used in Sec.
III.3 and Fig. 3. We stimulate with a biphasic charge-balanced pulse with a
period of positive stimulation lasting 0.2, period of no stimulation lasting
0.4, and a period of negative stimulation lasting 1.0 (pulse shape can be seen
compared to the signal in Fig. 9b,h. We used relatively long time series (1500
periods) and there is no noise or other unknown inputs. Additionally, we
determined events of equal phase with the optimal threshold 111111The optimal
threshold corresponds to thresholding with a Poincaré section that is
tangential to a local isochron, see Fig. 15. We find it by performing a direct
search over a reasonable range of threshold values and choosing the one that
corresponds to the lowest error value (15), as explained in the method paper
Cestnik and Rosenblum (2018). by performing a direct search over possible
threshold values, minimizing the error (15). As a result, the inferred curves
(red) accurately reflect the true ones (thick gray). The slight deviations are
mostly due to using a finite Fourier representation ($N_{F}=10$). As a
byproduct of the inference, we also obtain the asymptotic phase and isostable
amplitude as functions of time, and plot them against their true counterparts
in Fig. 9e,f,k,l. The obtained isostable time series is also a good
representation of the signal envelope if shifted to match the signal maxima,
as we show in a later Figure. 13.

Figure 9: Results of PRC and IRC inference with the IPID-1 method from Sec.
III.5 for the two test oscillators: Stuart-Landau (4) (top row) and its
generalization Eq. (8) (bottom row). From left to right: (a,g) trajectory in
state space, (b,h) an epoch of the observed signal and pulsatile forcing (blue
and green), (c,i) phase response curve, (d,j) isostable response curve, (e,k)
asymptotic phase comparison, (f,l) isostable amplitude comparison. We put the
labels for the vertical axes on the top of the plots to save space.

###  III.6 Inference in the presence of noise

Real-world oscillators are inevitably noisy, and therefore, we test the
performance of introduced techniques in the presence of dynamical noise. For
this purpose, we simulate the same two systems (4) and (8) with different
strengths $\sigma$ of Gaussian white noise. We start with the SL oscillator
and test the sine-fitting, Hilbert-based, and auxiliary oscillator techniques;
again we use the monophasic pulses. In Fig. 10 we illustrate the performance
of the latter method. (Since the sine-fitting and Hilbert-based techniques are
less efficient for noise-free systems, we do not illustrate them here.) We
show the inference errors as a function of the noise intensity $\sigma$, for
three different values of $\beta$, cf. Eq. (4). As expected, for this method
the PRC inference is more sensitive to noise than the ARC reconstruction.

Figure 10: Error (14) of PRC (a) and ARC (b) inference using an auxiliary
oscillator as explained in Sec. III.4.3. The error is measured for different
noise strengths $\sigma$.

Next, we proceed with testing the most promising IPID-1 technique, based on
the phase – isostable dynamics reconstruction. We keep the forcing strength
constant and then for each $\sigma$ generate several trajectory realizations
and for each one we compute the $L$ error measure using Eq. (14). These errors
are then depicted with standard box plots in Fig. 11a,b. For comparison, we
also show the results for the standard technique (PRC only). We conclude that
the PRC inference through the technique introduced in Section III.5.1 is most
stable to noise; in particular, it outperforms the standard approach. The
direct comparison of the amplitude response reconstruction is not that easy
since the auxiliary oscillator approach yields the empirical amplitude and,
correspondingly, the ARC, while the isostable reconstruction provides the IRC.
However, the results presented in Fig. 10b and Fig. 11b indicate approximately
equal performance. We underline that stronger forcing would imply lower errors
in the strong noise regime since the forcing-to-noise ratio would increase. On
the other hand, the errors in the weak noise regime would increase since the
linear approximation works worse for strongly forced systems. We also mention
that the standard technique (red) has stronger limitations on the
observations, namely, the pulses have to be rare (one pulse per several
periods), which is why in the case of no noise, the standard technique has a
smaller error. For our technique we considered a more realistic case of more
than one pulse per period on average. It means our technique requires much
shorter time series, see Appendix section D for details.

Figure 11: Performance of the standard technique III.3 (red, crossed) and
introduced IPID-1 methods III.5.1 and III.5.2 (green, plain) in relation to
dynamical Gaussian white noise strength $\sigma$ (the red crossed and green
plain box plots are slightly shifted horizontally to the left and to the
right, respectively, to avoid overlap). Two test oscillators were used,
Stuart-Landau (4) (a,b) and the modification (8) (c,d). Panels (a,c) and (b,d)
show the error measures according to Eq. (14) for the phase response
$Z(\varphi)$ and the amplitude response $I(\varphi)$, respectively. For each
value of noise $\sigma$ several trajectory realizations were simulated, each
realization yielding one error measurement. The $y$-axis is logarithmic, and
the measurements for the same value of noise are represented with standard box
plots. Parameters of the oscillators and the perturbing pulses are the same as
used before in Sec. III.3 and Figs. 3 and 9.

Finally, we test our phase – isostable reconstruction technique on the
complicated case of the modified SL oscillator (8). The results in Fig. 11c,d
confirm the advantage of this technique, namely, it is not constricted to
sinusoidal signals.

###  III.7 Response of a high-dimensional system

As a final test we apply the studied and newly-developed techniques to a
signal generated by a globally coupled oscillatory ensemble, see Eqs. (13).
The individual units are governed by the Bonhoeffer–van der Pol equations.
This system can be treated as a simple model of neuronal rhythmical activity.

We consider a large number of oscillators and choose the parameters such that
the system exhibits weak collective chaos, see test oscillators Sec. II for
more details on the system and Fig. 2 for a depiction of the mean-field orbit.
We assume that we observe collective dynamics. This example is a hard test. In
this case, not only is the ground truth phase model unknown to us, but since
the system is chaotic, we know that there exists no phase description that can
exactly describe the local deviations Schwabedal  _et al._ (2012). We,
therefore, estimate the goodness of inference solely based on how well the
inferred model reproduces observations; namely, we consider the ratio of
errors (15) and (16): $E_{Z}/E_{Z0}$. Our IPID-1 method performs best and
yields an error ratio of less than half, while other methods perform rather
poorly. Figure 12 shows all the corresponding phase response curves and the
values of errors in a bar plot. We generated the empirical PRCs using the mean
field in the $x$ variable as our observable, $s=X$, and exploiting biphasic
pulses as before. Next, we deconvolved all curves to obtain the theoretical
phase response curves. We performed the inference on a different observable
$s=X+2Y$ of the same system and obtained similar results.

Figure 12: Application of all phase response methods on a high dimensional
example of a chaotic oscillatory ensemble (13). (a) Inferred PRCs and (b) the
ratio of errors (15) and (16), representing the goodness of the inferred model
(ratio 0 would correspond to a perfect model, and ratio one is as good as
considering no response).

Only three techniques can infer the amplitude response: the Hilbert-based,
auxiliary oscillator-based, and IPID-1. The first two methods yield the
effective ARC. In contrast, our technique estimates the infinitesimal
isostable curve $I(\varphi)$. Hence, for comparison, we recompute $I(\varphi)$
into the ARC in the following way: first, using the infinitesimal curve
$I(\varphi)$ and pulse shape $\mathcal{P}(t)$ we evaluate the isostable shift
$\Delta\psi$. Then we relate it to the ARC. For a particular phase
$\varphi^{*}$, we integrate the first-order dynamics (1,3) with the considered
pulse shape $\mathcal{P}(t)$ to obtain phase as a function of time,
$\varphi(t)$, via Eq.(17), and then compute the effective isostable shift as:

$\Delta\psi(\varphi^{*})=\int\limits_{0}^{\delta}(\kappa\psi(t)+(\varphi(t))\mathcal{P}(t))\,\mathrm{d}t\;.$
(28)

Next, we have to relate the isostable variable $\psi$ to amplitude $a$. As
already mentioned, locally, the isostable linearly depends on the distance
from the limit cycle, see Eq. (22) : $\psi\approx
c(\varphi)(R-R_{0}(\varphi))$. Here $R$ represents the distance from the
origin, while $c(\varphi)$ and $R_{0}(\varphi)$ are phase-dependent and hold
the information of the isostable structure and parametrization of the limit
cycle, respectively. If amplitude is simply defined as $a=R$, then
$a=R_{0}(\varphi)+\psi/c(\varphi)$. Recalling the definition of ARC as a ratio
of amplitudes before and after a pulse, we write:

$A(\varphi)=\frac{a_{a}}{a_{b}}=1+\frac{\Delta\psi(\varphi)}{c(\varphi)R_{0}(\varphi)}\;.$
(29)

We cannot estimate $c(\varphi)$ nor $R_{0}(\varphi)$ from a scalar signal -
therefore, here we approximate them by constants. We choose $c=1$ since this
value was used in the inference, and for $R_{0}$ we use the maximal value of
the signal: $R_{0}=s_{\text{max}}$. 121212Ignoring the phase dependence of
$c(\varphi)$ is similar to ignoring the information of local isochrons in the
phase response, as most techniques do, e.g., the Hilbert-based one. We stress
that difficulties of recomputing IRC into ARC are due to an ad hoc definition
of the amplitude. Indeed, defined as the distance from the origin, it depends
on the projection space (depends on both the observable and embedding
technique). In contrast, $I(\varphi)$ is an invariant characteristic of a
limit cycle.

We compare the curves in Fig. 13a. Similarly to the characterization of the
phase response, we estimate the goodness of the IRC inference quantifying how
well the dynamical model fits the observation. Namely, we compute the ratio of
errors (26) and (27). We obtain $E_{I}/E_{I0}=0.28$, which indicates a decent
inference (we remind that this measure is 0 for a perfect fit and of order one
for a completely wrong curve). Note that two other techniques do not yield
dynamical equations; therefore, no errors were computed. As a byproduct, our
technique IPID-1 also yields an estimation of $\psi$ as a function of time. In
Fig. 13b,c, we demonstrate that the isostable variable provides a good
envelope if shifted to match the signal’s maxima (or minima). We depict such
envelope along with the Hilbert amplitude as the commonly used alternative.
Notice how Hilbert amplitude significantly fluctuates on the timescale of one
period, while the isostable envelope changes mainly in times of stimulation
and otherwise slowly follows the signal.

Figure 13: Inferred amplitude response (a) and the signal envelope (b,c) for
the high-dimensional chaotic ensemble (13). Three methods were used for
inferring the amplitude response: Hilbert-based, virtual oscillator-based, and
IPID-1. Since IPID-1 yields the theoretical IRC, we recompute it into ARC, see
text. In (b,c) we show two signal’s envelopes: the Hilbert amplitude (green
dash-dotted) and shifted isostable amplitude $\psi$ (red full line). Plot (b)
shows the behavior on a short time scale (few periods), and plot (c) shows it
on a long timescale (dozens of periods).

##  IV Discussion and conclusions

This paper addressed the problem of inferring an oscillator’s phase and
amplitude response from observations using test stimulation. In addition to
critically testing several techniques presented in the literature, we also
developed a new method for fitting first-order phase-isostable dynamics to
observations, denoted as the IPID-1 method. We discuss these two tasks of
phase and amplitude inference separately.

###  IV.1 Phase response

PRC determination for noise-free neuron-like oscillators has been known for
decades. We here concentrated on noisy signals without well-defined marker
events, such as spikes. As our first result, we mention the technique for
recomputing the empirical PRC obtained in response to stimuli of arbitrary,
but known shape, into the theoretical PRC, $Z_{\cal P}(\varphi)\to
Z(\varphi)$. This technique can be used along with any approach for PRC
inference, in particular, it can enhance the applicability of the standard
approach. This addition is especially beneficial for analyzing biological
systems where some restrictions on the stimulus’s shape may apply.

Next, we suggested a simple approach to inference’s error estimation. Without
knowing the ground truth and using only the observations, we quantify the
inferred PRC’s usefulness. Namely, we quantify how much of the signal’s
variability can be explained by the inferred PRC. Again, this approach can
complement any PRC estimation technique. We strongly recommend using this
error estimation tool in experiments. Indeed, blind application of an
inference technique always provides some response curve, but in the case of
noisy or chaotic systems, the obtained curve may not describe the underlying
dynamics.

We critically considered previously developed techniques that were introduced
without being tested on examples where the ground truth is known. We performed
tests and highlighted the potential weaknesses and drawbacks. In particular,
we suggested improving the inference approach based on the popular analytical
signal technique exploiting the HT by ignoring the data points immediately
before and after stimuli. Furthermore, we compared the performance of the
known and newly developed techniques in the noise-free and noisy cases using
specially designed test data with the known ground truth. We tested the
dependence of the results on the noise level and the observable used for data
analysis. In particular, we have shown that the HT-based technique is
generally unreliable. As a result, we have demonstrated the essential
advantage of our IPID-1 approach based on the direct reconstruction of the
Winfree equation: it is not restricted to sine-like signals, it immediately
provides the theoretical PRC, and is robust against noise. We confirmed this
conclusion by applying different inference techniques to a high-dimensional
system representing the simplest model of brain rhythm generation. Estimating
the error, we have demonstrated that IPID-1 is the only approach to yield a
good result for this challenging test.

We conclude the discussion of the phase response inference with the following
remark. While testing the algorithms on the model systems, we ensured that the
stimulation was weak. In practice, one must perform stimulation with stimuli
of different strengths and check whether the inferred PRC depends on this
strength. No dependence means validity of the linear approximation and, hence,
of the PRC description. Otherwise, the revealed curve is not the theoretical
(infinitesimal) PRC but quantifies the effect of a finite-strength
perturbation and is amplitude-dependent.

###  IV.2 Amplitude response

The main problem with the amplitude response’s inference is the definition of
the amplitude. While phase can be uniquely defined for any point in the basin
of attraction of a limit cycle, the definition of the amplitude is ambiguous.
An operational approach introduces the (empirical) amplitude variable as an
envelope to the observable and computes the envelope employing, e.g., the
Hilbert Transform. However, this definition obviously depends on the
observable. Moreover, for non-harmonic signals, HT and other techniques
generally do not provide a “good” envelope since they show variation for
perfectly periodic signals.

Another approach involves equating the amplitude with the slowest decaying
isostable variable; eigenvariable of the linearized system, which as such
provides a universal description independent from the choice of observable.
The IPID-1 approach developed here reconstructs the first-order isostable
dynamics directly from a scalar time series. The advantage of this technique
is that inferred Eq. (3) can be further exploited to predict the effect of
stimuli of a different shape. Furthermore, like the PRC inference, the
technique provides an error measure that can be computed solely from data. We
also mention the limitations of the procedure. It seems that inferring the
isostable response is generally harder than inferring the phase response, and
this is reflected in the performance of the IPID-1 technique by, e.g.,
comparing the inference where the ground truth is known. We also stress that
Eq. (3) suffices for describing the dynamics only if the oscillator is two-
dimensional or generally high-dimensional, but perturbations’ decay in one
direction is much slower than in other directions. The reduction of limit-
cycle oscillators to phase-isostable dynamics is rather recent Wilson and
Moehlis (2016b) and the physical interpretation of the isostable structure is
still a matter of discussion. Another problem is to relate the isostable
variables to envelopes of observed signals. Generally, inference and
quantification of the amplitude response remains a subject of further studies.

###  IV.3 Relevance for oscillatory dynamics control

Knowledge of the phase and amplitude response allows efficient control of an
oscillatory system. Indeed, probing the system by pulses of a known shape, we
infer the Winfree Eq. (1) and then can exploit this equation to optimize the
stimulus’ shape, e.g., aiming to maximize the response.

In many cases, the goal of the control is to suppress or enhance the
oscillation. The most illustrative example is deep brain stimulation (DBS),
aiming to quench the Parkinsonian tremor. The design of efficient control
schemes requires the determination of vulnerable phases where the stimulation
is most efficient. A possible approach relies on an adaptive control scheme
Montaseri  _et al._ (2013); Rosenblum (2020). An alternative idea is to use
first a test stimulation to infer the response properties and then exploit the
corresponding curves to determine the proper stimulation phases Duchet  _et
al._ (2020). Maxima and minima of the IRC or ARC provide the optimal
stimulation phases for correspondingly enhancing and suppressing the
oscillation. We emphasize that PRC alone does not yield the required
information. Although a relationship between PRC and ARC has been demonstrated
for some examples Duchet  _et al._ (2020), generally, this relationship does
not hold, see Appendix A.3 for details. Another idea is stimulating around
zeros of $Z(\varphi)$. If the phase response is zero, the stimulus acts along
the isochron and modifies only the amplitude. This idea might be a good
starting point, but it does not guarantee that this phase is optimal \- the
final amplitude variation depends on isostables density. We believe that
reconstruction of the isostable equation from data provides a means to design
an optimal stimulation. Just as in the case of phase response, using the
reconstructed Eq. (3) we obtain the amplitude responses to arbitrary stimuli
and can exploit this equation to optimize the stimulus’ shape.

###### Acknowledgements.

R.C. and E.T.K.M. acknowledge financial support from Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation), Project-IDs PI
220/21-1 and 424778381 – TRR 295, respectively.

##  Appendix A Phase – isostable coordinates

Isostable coordinates extend the phase description of a limit cycle by adding
variables $\psi_{i}$ that describe the decay or growth of deviations from the
limit cycle in the isochronal hyperplane. In the vicinity of a limit cycle,
any $N$-dimensional system can be represented in the isostable coordinate
system

$\displaystyle\dot{\varphi}$ $\displaystyle=\omega\;,$ (30)
$\displaystyle\dot{\psi_{i}}$ $\displaystyle=\kappa_{i}\psi_{i}\;,\qquad
i=1,\dots,N-1\;,$ (31)

where $\omega$ is the frequency and $\kappa_{i}$ are Floquet exponents of the
limit cycle. For a system with a stable limit cycle, all Floquet exponents
have a negative real part, $\Re(\kappa_{i})<0$. However, in many cases, it is
sufficient to consider only one isostable variable $\psi_{\max}$ corresponding
to the largest real part of its Floquet exponent. This two-dimensional
approximation is justified if there is a clear separation of time scales and
all neglected $\psi_{i}$ decay considerably faster than $\psi_{\max}$.
Restricting our consideration to only one isostable coordinate $\psi_{\max}$
complementing $\varphi$, we also assume the corresponding Floquet exponent to
be real. It means that we do not address the case of the slowest mode
described by a complex-conjugate pair of Floquet exponents. Thus
$\kappa:=\max_{i}\Re(\kappa_{i})$ and below we omit the subscript ${\max}$ in
the notation $\psi_{\max}$. For a detailed discussion, in particular, on how
to interpret a general $N$-dimensional system of phase-isostable coordinates
with pairs of complex conjugate Floquet exponents, see Wilson and Moehlis
(2016b); Wilson and Ermentrout (2018a); Wilson (2019).

As one can check from the autonomous dynamical Eqs. (31), all variables
$\psi_{i}$ are invariant to scaling with some non-zero factor $c$ – the
dynamical equations remain unchanged after a change of variables
$\tilde{\psi}:=c\psi$. As soon as external perturbations $p(t)$ are introduced
to the system via

$\displaystyle\dot{\varphi}$
$\displaystyle=\omega+G_{\varphi}(\varphi,\psi)p(t)\,,$ (32)
$\displaystyle\dot{\psi}$
$\displaystyle=\kappa\psi+G_{\psi}(\varphi,\psi)p(t)\,,$ (33)

the scaling of isostable variable leads to the adjustment in the stimulation
function

$\displaystyle
G_{\tilde{\psi}}(\varphi,\tilde{\psi})=cG_{\psi}(\varphi,\frac{\tilde{\psi}}{c})\;.$
(34)

Thus, the IRC (defined as the isostable amplitude response evaluated at the
LC) is also scaled as a result of the scaling of $\psi$:

$\displaystyle
I_{\tilde{\psi}}(\varphi)=G_{\tilde{\psi}}(\varphi,0)=cI_{\psi}(\varphi)\,.$
(35)

The PRC is invariant under rescaling of the isostable $\psi$. In the following
derivations, we keep the scaling factor $c$ to show where it appears in the
equations.

###  A.1 Isostable reduction for the SL system

The Stuart-Landau system, Eq. (4), is the analytically solvable normal form of
a Hopf bifurcation: For $\mu>0$, a radial limit cycle at $R_{0}=\sqrt{\mu}$
establishes with a basin of attraction being the entire phase space except for
the fixed point at the origin. Thus, each point in the phase space, here given
in polar coordinates $R$ and $\theta$, can be assigned a phase by the well-
known formula

$\displaystyle\varphi(R,\theta)=\theta-\alpha\ln\left(\frac{R}{\sqrt{\mu}}\right)\,,$
(36)

for which one can verify that $\dot{\varphi}=\eta-\alpha\mu=:\omega$. Since
the SL system is two-dimensional, one additional isostable $\psi$ suffices to
describe the dynamics in the vicinity of the limit cycle completely. Hence,
for the SL example, the isostable coordinate system is a one-to-one
transformation for $\mathbb{R}\setminus\\{0\\}$. For the case of radial
isochrons ($\alpha=0$), $\psi$ has already been derived in Ref. Wilson and
Ermentrout (2018b) as

$\displaystyle\psi(R,\theta)=c(1-\frac{\mu}{R^{2}})\,.$ (37)

Moreover, this coordinate transformation yields
$\dot{\psi}=-2\mu\psi:=\kappa\psi$ also in the case of non-radial isochrons
($\alpha\neq 0$). From this expression, one can check the system’s dynamical
properties: the roots of $\psi$ at $R=\sqrt{\mu}$ yield the coordinates of the
limit cycle, and the divergence at $R=0$ indicates the border of the limit
cycle’s basin of attraction. Also, as $\psi$ converges to a finite value as
$R\rightarrow\infty$, we deduce that it takes a finite time to reach the
vicinity of the limit cycle for initial conditions with arbitrarily large $R$.

The SL system can be re-written using the dynamical parameters $\omega$ and
$\kappa$ instead of $\eta$ and $\mu$ as

$\displaystyle\dot{x}$ $\displaystyle=-\omega
y-(x^{2}+y^{2}+\frac{\kappa}{2})(x-\alpha y)+G_{x}(x,y)p(t)\;,$ (38)
$\displaystyle\dot{y}$ $\displaystyle=\leavevmode\nobreak\
\leavevmode\nobreak\ \omega x-(x^{2}+y^{2}+\frac{\kappa}{2})(y+\alpha
x)+G_{y}(x,y)p(t)\;,$

where the functions $G_{x,y}(x,y)$ specify how the perturbation acts on the
oscillator.

The knowledge of the analytical form of the coordinate transformation of
$\psi$ and $\varphi$ allows us to compute the PRC $Z(\varphi)$ and IRC
$I(\varphi)$ also in a closed analytical form. Those are, by definition, the
response curves of phase and isostable variables evaluated at the limit cycle,
where $\psi=0$. In general, the response curves are computed via

$\displaystyle\begin{pmatrix}I(\varphi)\\\ Z(\varphi)\\\
\end{pmatrix}=\begin{pmatrix}G_{\psi}\\\ G_{\varphi}\\\
\end{pmatrix}\Big{|}_{\psi=0}=J_{i\leftarrow p}\big{|}_{\psi=0}\cdot
J_{p\leftarrow c}\big{|}_{\psi=0}\cdot\begin{pmatrix}G_{x}\\\ G_{y}\\\
\end{pmatrix}\big{|}_{\psi=0}\;,$ (39)

where $J_{i\leftarrow p}$ is the Jacobian of the coordinate transformation
from polar to isostable coordinates, and $J_{p\leftarrow c}$ is the well-known
Jacobian of the coordinate transformation from Cartesian to polar coordinates:

$\displaystyle J_{i\leftarrow
p}=\begin{pmatrix}\partial_{R}\psi&\partial_{\theta}\psi\\\
\partial_{R}\varphi&\partial_{\theta}\varphi\\\
\end{pmatrix}=\begin{pmatrix}\frac{2\mu c}{R^{3}}&0\\\ -\frac{\alpha}{R}&1\\\
\end{pmatrix}\quad\;,\quad J_{p\leftarrow
c}=\begin{pmatrix}\partial_{x}R&\partial_{y}R\\\
\partial_{x}\theta&\partial_{y}\theta\\\
\end{pmatrix}=\begin{pmatrix}\cos(\theta)&\sin(\theta)\\\
-\frac{\sin(\theta)}{R}&\frac{\cos(\theta)}{R}\\\ \end{pmatrix}\,.$ (40)

Throughout this paper, we assume a state-independent stimulation, namely,
$G_{x}=\cos(\beta)$ and $G_{y}=\sin(\beta)$. Hence we arrive at the following
phase and isostable response curves for the SL system:

$\displaystyle I(\varphi)$
$\displaystyle=\frac{2c}{\sqrt{-\frac{\kappa}{2}}}\cos(\varphi-\beta)\;,$ (41)
$\displaystyle Z(\varphi)$
$\displaystyle=-\frac{1}{\sqrt{-\frac{\kappa}{2}}}\left(\sin(\varphi-\beta)+\alpha\cos(\varphi-\beta)\right)\;.$
(42)

For the rotationally invariant SL system, it is convenient to interpret the
polar radius $R$ as the amplitude of the limit cycle oscillation. Thus, from
the response curve of $R$, denoted as $G_{R}$, we infer the ARC as

$\displaystyle
A(\varphi)=\frac{R_{a}}{R_{b}}\approx\frac{R_{b}+fG_{R}(\varphi)}{R_{b}}=1+\frac{fG_{R}(\varphi)}{R_{b}}\;,$
(43)

approximating a short pulse with action $f$ by a Dirac delta function, hitting
the system at phase $\varphi$. For the SL system, the initial amplitude
$R_{b}=\sqrt{\mu}$ is independent of $\varphi$. For a stimulation that is
state-independent in Cartesian coordinates, the response in $R$ is given by
$G_{R}(\varphi)=\cos(\varphi-\beta)$ and thus the ARC simplifies to the
expression

$\displaystyle A(\varphi)\approx 1+\frac{f\cos(\varphi-\beta)}{\sqrt{\mu}}\;.$
(44)

###  A.2 Test model with a non-sinusoidal solution and known phase and
isostable response

In order to construct more complex two-dimensional test models with known
properties, we use the isostable coordinate system. Instead of trying to
derive the phase and isostable variable from known systems, we give an
explicit analytical expression of $\varphi$ and $\psi$ in the first place.
Since the roots of $\psi(\mathbf{x})$ are the location of the stable invariant
set for $\kappa<0$, we can control the position and shape of the limit cycle.
In particular, the coordinate transformation

$\displaystyle\psi(R,\theta)=c\Big{(}1-\frac{q(\theta)}{R^{2}}\Big{)}$ (45)

describes a limit cycle at $R_{0}=\sqrt{q(\theta)}$ for any strictly positive,
$2\pi$-periodic function $q$ of $\theta$. Together with the phase defined as

$\displaystyle\varphi(R,\theta)=\theta-\alpha\ln\left(\frac{R}{\sqrt{q(\theta)}}\right)\;,$
(46)

such that $\varphi=\theta$ on the limit cycle, we obtain a generalized Stuart-
Landau system. For $q(\theta)=-\frac{\kappa}{2}=\mu$ we recover the original
SL system. Since the Jacobian $J_{i\leftarrow p}$ can be computed in the same
straightforward manner as for the SL system as

$\displaystyle J_{i\leftarrow
p}\big{|}_{\psi=0}=\begin{pmatrix}\partial_{R}\psi&\partial_{\theta}\psi\\\
\partial_{R}\varphi&\partial_{\theta}\varphi\\\
\end{pmatrix}\bigg{|}_{\psi=0}=\begin{pmatrix}\frac{2cq(\theta)}{R^{3}}&-\frac{cq^{\prime}(\theta)}{R^{2}}\\\
-\frac{\alpha}{R}&1+\frac{\alpha q^{\prime}(\theta)}{2q(\theta)}\\\
\end{pmatrix}\bigg{|}_{\psi=0}=\begin{pmatrix}\frac{2c}{\sqrt{q(\theta)}}&-\frac{cq^{\prime}(\theta)}{q(\theta)}\\\
-\frac{\alpha}{\sqrt{q(\theta)}}&1+\frac{\alpha
q^{\prime}(\theta)}{2q(\theta)}\\\ \end{pmatrix}\quad\;,$ (47)

we can explicitly write the dynamical equations, both in polar

$\displaystyle\dot{R}$
$\displaystyle=\frac{R}{2q(\theta)}\left(\kappa\left(1+\frac{\alpha
q^{\prime}(\theta)}{2q(\theta)}\right)\left(R^{2}-q(\theta)\right)+\omega_{0}q^{\prime}(\theta)\right)\;,$
(48) $\displaystyle\dot{\theta}$
$\displaystyle=\omega_{0}+\frac{\alpha\kappa}{2q(\theta)}(R^{2}-q(\theta))\;,$
(49)

and in Cartesian coordinates (where $q$ and $q^{\prime}$ have to be read as
functions of $\theta(x,y)$)

$\displaystyle\dot{x}$
$\displaystyle=\omega\left(\frac{xq^{\prime}}{2q}-y\right)+\frac{\kappa}{2}\left(\frac{x^{2}+y^{2}}{q}-1\right)\left(x+\alpha\left(\frac{xq^{\prime}}{2q}-y\right)\right)\;,$
(50) $\displaystyle\dot{y}$
$\displaystyle=\omega\left(\frac{yq^{\prime}}{2q}+x\right)+\frac{\kappa}{2}\left(\frac{x^{2}+y^{2}}{q}-1\right)\left(y+\alpha\left(\frac{yq^{\prime}}{2q}+x\right)\right)\;.$
(51)

Having the explicit Jacobians for the coordinate transformation, we obtain the
PRC and IRC generally given by the response curves in polar angle $G_{\theta}$
and polar radius $G_{R}$ as

$\displaystyle I(\varphi)$ $\displaystyle=2c\mathcal{I}(\varphi)\;,$ (52)
$\displaystyle Z(\varphi)$
$\displaystyle=G_{\theta}\big{|}_{\psi=0}(\varphi)-\alpha\mathcal{I}(\varphi)\;,$
(53)

with the function $\mathcal{I}$ defined as

$\displaystyle\mathcal{I}(\varphi)$
$\displaystyle=q^{-\frac{1}{2}}(\varphi)G_{R}\big{|}_{\psi=0}(\varphi)-\frac{q^{\prime}(\varphi)}{2q(\varphi)}G_{\theta}\big{|}_{\psi=0}(\varphi)\;.$
(54)

The polar angle $\theta$ serves as a protophase in this model since, on the
limit cycle, $\theta$ and $\varphi$ coincide despite not being equal in
general. Thus, on the limit cycle where $\psi=0$, we have $\theta=\varphi$ and
we obtain the shape of the limit cycle $R_{0}=\sqrt{q(\varphi)}$ as a function
of phase $\varphi$.

Note the relation from Eqs. (53), stating that the IRC is proportional to the
difference of PRC and the response curve of the polar angle
$G_{\theta}\big{|}_{\psi=0}$ for non-zero non-isochronicity parameter
$\alpha$:

$\displaystyle Z(\varphi)-G_{\theta}\big{|}_{\psi=0}(\varphi)$
$\displaystyle=-\frac{\alpha}{2c}I(\varphi)\;.$ (55)

We emphasize this relation between PRC and IRC since it exists without further
assumptions made about $G_{\theta}$ and $G_{R}$, i.e., how external
stimulation enters the system, but is exclusively due to the autonomous
dynamics of the system. However, despite the generalizing extension this model
has added to the simple SL model, it still holds strong assumptions that lead
to this conclusion, e.g., that phase $\varphi$ and polar angle $\theta$ are
equal on the limit cycle.

Assuming that the external stimulation enters state-independently in Cartesian
coordinates, we compute the polar angle response function as

$\displaystyle
G_{\theta}\big{|}_{\psi=0}(\varphi)=-\frac{\sin(\varphi-\beta)}{\sqrt{q(\varphi)}}$
(56)

and $\mathcal{I}$ as

$\displaystyle\mathcal{I}(\varphi)$
$\displaystyle=\frac{2\cos(\varphi-\beta)q(\varphi)+\sin(\varphi-\beta)q^{\prime}(\varphi)}{2q^{\frac{3}{2}}(\varphi)}\;.$
(57)

For the particular example in Sec. II, we use $q(\theta)=r+2\cos^{2}(\theta)$,
where $r$ is a positive parameter. We denote this particular model as the
modified Stuart-Landau (mSL) model since it does not contain the original SL
model as a special case. Its polar angle response curve is given by

$\displaystyle
G_{\theta}^{\text{mSL}}\big{|}_{\psi=0}(\varphi)=-\frac{\sin(\varphi-\beta)}{\sqrt{r+2\cos^{2}(\varphi)}}$
(58)

and the function $\mathcal{I}$ as

$\displaystyle\mathcal{I}(\varphi)=\frac{(r+1)\cos(\varphi-\beta)+\cos(3\varphi-\beta)}{(r+2\cos^{2}(\varphi))^{\frac{3}{2}}}\;.$
(59)

Having these two function, the PRC in Eq. (12) and IRC in Eq. (11) are derived
using the relations given in Eq. (53) for $c=1$.

###  A.3 Relationship of phase and amplitude response curves

Duchet et al. Duchet  _et al._ (2020) hypothesized that the amplitude response
is proportional to $\partial_{\varphi}{Z(\varphi)}$. Indeed, Figure 14a
demonstrates that this hypothesis holds for the isochronous SL oscillator. In
the following, we sketch how geometrical assumptions about an oscillating
system can lead to this hypothesis on its own: Assume (as we also did
throughout this paper) that in a two-dimensional system, the external
stimulation is state-independent in Cartesian coordinates $x$ and $y$. Then,
without loss of generality, the response functions can be written as
$G_{x}=\rho\cos(\beta)$ and $G_{y}=\rho\sin(\beta)$, where the angle $\beta$
determines the direction of the stimulation in the $x$-$y$-plane and the non-
negative value $\rho$ determines the strength. Without loss of generality, we
can set $\rho=1$ by absorbing the stimulation strength as a factor into the
stimulation shape function $p(t)$. Using a coordinate transformation from
Cartesian to polar coordinates, see Eq. (40), the response functions in the
polar angle $\theta$ and the polar radius $R$ follow as

$\displaystyle G_{\theta}=-\frac{1}{R}\sin(\theta-\beta)\qquad\text{and}\qquad
G_{R}=\cos(\theta-\beta)\;.$ (60)

These two functions establish a number of relations, e.g.,

$\displaystyle G_{R}=-R\cdot\partial_{\theta}G_{\theta}\qquad\text{or}\qquad
R^{2}G^{2}_{\theta}+G_{R}^{2}=1\,.$ (61)

Thus, we obtain the relation that $G_{R}$ is proportional to the derivative of
$G_{\theta}$. However, we emphasize that up to this point, the considerations
are purely geometrical and are valid regardless of the system’s dynamics.

If the limit cycle rotates around the origin, one sometimes tries to
approximate the system’s phase $\varphi$ by the polar angle $\theta$ and the
amplitude by the polar radius $R$. Indeed, for the standard SL system with its
circular limit cycle, the polar radius $R$ can serve as a measure of
amplitude, leaving the isostable response curve proportional to the response
in $R$ and proportional to the ARC shifted by one: $I\propto G_{R}\propto
A-1$. In addition, for the isochronous SL system ($\alpha=0$, Fig. 14a) the
phase $\varphi$ equals the polar angle $\theta$ and the geometric relation
(61) directly translates into a relation between PRC and IRC:
$I\propto\partial_{\varphi}Z$.

For the non-isochronous SL system illustrated in Fig. 14b, this relation is no
longer strict since phase $\varphi$ and polar angle $\theta$ differ outside
the limit cycle. Consequently, $Z\neq G_{\theta}$ and $\partial_{\varphi}Z$ is
only correlated instead of strictly proportional to $I$. The modified SL
system (Fig. 14c) has no circular limit cycle, and hence no straightforward
representation of amplitude by $R$, thus providing a counter-example to the
hypothesis.

Figure 14: Relationships of phase and isostable response curves for the SL
(Eqs. 5 and 6) and mSL (Eqs. 12 and 11) systems. In all panels, the blue
dashed curve depicts the PRC and the red solid curve its derivative as a
function of the IRC, parametrized by phase $\varphi$. Panel (a): isochronous
SL system with $\mu=0.2$ and $\alpha=0$. Panel (b): non-sochronous SL system
with $\mu=0.2$ and $\alpha=0.6$. Panel (c): mSL system with $r=0.75$ and
$\alpha=0$. For the SL oscillator, the amplitude can be expressed in terms of
the polar radius $R$. Thus, in panel (a) and (b) the IRC is proportional to
$A-1$ where $A$ is the ARC from Eq. (7). For the mSL depicted in panel (c), a
straightforward definition of amplitude does not exist.

##  Appendix B Relationship between empirical and infinitesimal PRC

Here, we discuss how to recompute the empirical PRC $Z_{\cal P}$ into
theoretical infinitesimal PRC $Z$. The phase shift evoked by a pulse of finite
width $\delta$ applied at the instant $t_{s}$ is
$\Delta\varphi=\varphi(t_{e})-\varphi(t_{s})-\omega\delta$, where
$t_{e}=t_{s}+\delta$. Hence, the Winfree Eq. (1) yields
$\Delta\varphi=\int_{t_{s}}^{t_{e}}Z(\varphi(t)){\cal P}(t)\,\mathrm{d}t$. If
forcing is weak: $Z(\varphi)p(t)\ll\omega$, then in the first approximation we
take $\varphi\approx\omega t$ and obtain:

$Z_{\cal
P}\Big{(}\varphi(t_{s})\Big{)}=\int_{t_{s}}^{t_{e}}Z\Big{(}\varphi(t_{s})+\omega
t\Big{)}{\cal P}(t)\,\mathrm{d}t\;.$ (62)

This relation is a convolution of the true PRC $Z$ and the pulse’s waveform
${\cal P}$: $Z_{\cal P}=Z\star{\cal P}$. The convolution operation can be
expressed with a product relation in Fourier space: ${\cal F}(Z_{\cal
P})={\cal F}(Z)\cdot{\cal F}({\cal P})$, which we utilize for expressing
$Z(\varphi)$:

$Z(\varphi)={\cal F}^{-1}\Big{(}{\cal F}(Z_{\cal P})/{\cal F}({\cal
P})\Big{)}\;,$ (63)

where ${\cal F}$ and ${\cal F}^{-1}$ denote the direct and inverse Fourier
Transform, respectively.

We mention, that for a charge-balanced pulse we do not recover the constant
term in $Z$. Indeed, substituting in Eq. (62) $Z$ by $Z+C$, where $C$ is an
arbitrary constant, we obtain the same function $Z_{\cal P}$ due to the charge
balance condition $\int_{t_{s}}^{t_{e}}{\cal P}\,\mathrm{d}t=0$. To solve this
problem, we set the mean value of $Z$ equal to the mean value of $Z_{\cal P}$.

An alternative solution for the inverse problem that does not assume
$\varphi\approx\omega t$ is as follows. We represent the yet unknown function
$Z(\varphi)$ as a finite Fourier series with coefficients $a_{k},b_{k}$. Then
using $\mathcal{P}(t)$ we express $Z_{\mathcal{P}}(\varphi)$ in terms of
$a_{k},b_{k}$ by numerically solving the Winfree equation. The nonlinear
problem of equating the empirical curve with its symbolic representation:
$Z_{\cal P}=Z_{\cal P}(\varphi;a_{k},b_{k})$, is then solved using the
Levenberg-Marquardt technique .

##  Appendix C Thresholding a signal

Here, we discuss how to determine events corresponding to the same asymptotic
phase from the observed signal $s(t)$. This task is easy if the system moves
on the limit cycle; then the signal $s(t)$ is a projection of a one-
dimensional trajectory, and any defining event that occurs once per period
will correspond to the same state and, therefore, the same phase. Generally,
defining events can be threshold crossings, extrema, and inflection points.

If the system is not on the limit cycle, determining the asymptotic phase is
challenging. Generally, one would have to know the isochronal structure and
then determine events as crossings of a particular isochron. (We remind that
the isochrons are the manifolds of constant phase.) However, even with a
priori knowledge of isochrons, one cannot determine the crossing of a
potentially high-dimensional manifold while observing a scalar signal $s(t)$
only. Therefore, we restrict our consideration to 2D oscillators or higher
dimensional oscillators that have one Floquet exponent much smaller in
absolute value than the rest (having a 2D slow manifold). In general, for
higher dimensional systems the approximation will be less accurate.
Additionally we mention, that one does not need to find the events exactly,
but just their good approximation.

Let us consider threshold crossings such as, e.g.,
$s(t)=s_{\text{thr}},\frac{d}{dt}s(t)>0$. Fig. 15 illustrates an example.
Since we are thresholding a scalar signal $s(t)$ we only have one parameter:
the threshold value $s_{\text{thr}}$. Different thresholds correspond to
different parallel straight lines (see two examples in Fig. 15 depicted in red
and blue). The best approximation in this situation is when the line is
tangential to a particular isochron at the limit cycle. Then, the
corresponding threshold provides a first-order approximation for the states
with the same phase. While we cannot change the inclination of the line, we
can choose its vertical position. For 2D systems, there always exist at least
two such values, one when the signal is growing (see red dashed line in Fig.
15) and one when it is falling (blue dash-dotted). From this illustration, one
can deduce that extremal points are generally poor approximations for equal-
phase states since the corresponding lines are (close to) tangential to the
limit cycle and thus unlikely to be tangential to any isochron.

Figure 15: An illustration of choosing a proper threshold. (a) The isostable
structure of a 2D limit-cycle oscillator. The limit cycle is depicted with
black and the isochrons are shown with yellow. The two thresholds that are
tangential to isochrons are shown in red and blue. (b) The corresponding
signal $x(t)$ is depicted with black. The two signals in gray correspond to a
motion along a non-zero isochron. Note that the signal and two non-zero
isochron signals all coincide at the ideal threshold crossings. Panels (a,b)
share the vertical axis.

##  Appendix D Test perturbation used in our experiments

Here, we discuss the test signals we apply to the test oscillators in our
experiments. Except for the cases where we mention it explicitly, we
consistently use a biphasic charge-balanced pulse with a period of positive
forcing lasting 0.2, no forcing lasting 0.4, and a period of negative forcing
lasting 1.0. The overall duration of the pulse is therefore 1.6, and since
both our test oscillators (4) and (8) are chosen to have a period 2$\pi$ each
pulse lasts approximately $1/4$ of a period, see Fig. 3 for a depiction of the
pulse relative to the oscillators phase response curves. We use different
forcing action $f=\frac{1}{2}\int_{0}^{\tau}|\mathcal{P}|dt$ for different
oscillators: for the Stuart-Landau system (4) we used $f=0.01$ while for the
generalization (8) we used $f=0.07$ since it has a larger limit cycle and thus
needs a stronger perturbation for roughly the same effect.

For the test with the standard inference technique, we introduced the pulses
periodically with a much larger period (roughly five times) than the period of
oscillators 2$\pi$, ensuring that the two periods are incommensurate; hence
the oscillator would eventually be stimulated at almost every phase. Such
perturbation is ideal for inferring the phase response since the pulses appear
rarely, giving the oscillator time to relax to the limit cycle where the phase
shift is easy to read. In contrast, for our introduced techniques
III.5.1,III.5.2 we devised a more realistic experimental example where the
pulses are coming randomly and often. We sample the inter-pulse times from a
Poisson distribution such that, on average, one oscillator period receives 1.6
stimulating pulses, see Fig 9b,h for a short signal vs. forcing example.

##  Appendix E Causal phase and amplitude estimation exploiting virtual
oscillators

Here we show how to estimate phase and amplitude through a virtual oscillator.
We explain how to obtain the amplitude first. We choose the oscillator’s
frequency $\eta$ to be much larger than the characteristic frequency $\nu$ of
$s(t)$ and the damping parameter $\alpha_{a}$ sufficiently large so that the
amplitude response is practically independent of $\nu$. Then, we recompute the
amplitude of the virtual oscillator (that we can monitor) into the amplitude
of the input signal, using the well-known expression for the resonance curve.
Namely, we obtain the input’s amplitude as

$a(t)=\sqrt{x^{2}+(\dot{x}/\nu)^{2}}\sqrt{(\eta^{2}-\nu^{2})^{2}+(\alpha_{a}\nu)^{2}}\;.$

For phase estimation, we choose the damping parameter $\alpha_{\varphi}$ to be
small so that the phase shift
$\beta=\arctan[-\alpha_{\varphi}\nu/(\eta^{2}-\nu^{2})]$ between the virtual
oscillator and its input is nearly zero and weakly dependent on $\nu$ for
$\nu\ll\eta$. Then, we compute input’s phase as
$\varphi(t)=\arctan(-\dot{x}/\nu x)-\beta$. For further detail and an
efficient numerical scheme, we refer to Ref. Rosenblum  _et al._ (2021).

##  Appendix F Isostable inference step-by step summary

  1. 1.

Determine events $\tau_{i}$ by thresholding the phase, e.g.,
$\varphi(\tau_{i})=\varphi_{\text{thr}}$, $\frac{d}{dt}s(\tau_{i})>0$. The
choice of threshold $\varphi_{\text{thr}}$ is not crucial but generally should
be picked such that the signal in events $s(\tau_{i})$, has a wide range of
values.

  2. 2.

Write system (23) by linearly approximating the isostable integral
$\int_{\tau_{i}}^{\tau_{i+1}}(\psi(t)+s_{0})dt\approx[s(\tau_{i})+s(\tau_{i+1})]/2$:

$s(\tau_{i+1})-s(\tau_{i})=-\kappa
s_{0}(\tau_{i+1}-\tau_{i})+\kappa\frac{s(\tau_{i})+s(\tau_{i+1})}{2}+\int\limits_{\tau_{i}}^{\tau_{i+1}}I(\varphi(t))p(t)\,\mathrm{d}t\;.$

  3. 3.

Express the response function $I(\varphi)$ as a Fourier series:
$I(\varphi)=\sum_{n=0}^{N_{F}}[u_{n}^{\cos}\cos(n\varphi)+u_{n}^{\sin}\sin(n\varphi)]$,
and replace the order of integration and summation in the last term to express
it as a series of computable integrals:

$\int\limits_{\tau_{i}}^{\tau_{i+1}}I(\varphi(t))p(t)\,\mathrm{d}t=\sum\limits_{n=0}^{N_{F}}[u_{n}^{\cos}P_{n,i}^{\cos}+u_{n}^{\sin}P_{n,i}^{\sin}]\;,$

where we denoted the Fourier integrals with
$P_{n,i}^{\cos}=\int_{\tau_{i}}^{\tau_{i+1}}\cos(n\varphi(t))p(t)\,\mathrm{d}t$
and
$P_{n,i}^{\sin}=\int_{\tau_{i}}^{\tau_{i+1}}\sin(n\varphi(t))p(t)\,\mathrm{d}t$.
Numerically compute the integrals for each inter-event interval.

  4. 4.

Combine the approximations to express system (23) as a linear system with
known coefficients:

$s(\tau_{i+1})-s(\tau_{i})=-\kappa
s_{0}(\tau_{i+1}-\tau_{i})+\kappa\frac{s(\tau_{i})+s(\tau_{i+1})}{2}+\sum_{n=0}^{N_{F}}[u_{n}^{\cos}P_{n,i}^{\cos}+u_{n}^{\sin}P_{n,i}^{\sin}]\;.$

There are as many equations as there are inter-event intervals. The known
coefficients are: $s(\tau_{i+1})-s(\tau_{i})$, $\tau_{i+1}-\tau_{i}$,
$(s(\tau_{i})+s(\tau_{i+1}))/2$, and the Fourier integrals $P_{n,i}$. The
unknown quantities are: $\kappa s_{0},\ \kappa$ and the Fourier modes
$u_{n}^{\cos},u_{n}^{\sin}$, representing the response $I(\varphi)$. Thus
proceed to minimize this system via least squares (or a similar method) to
obtain the first approximated solution:
$s_{0}^{(1)},\kappa^{(1)},I^{(1)}(\varphi)$.

  5. 5.

Iterate the following: approximate the continuous amplitude $\psi(t)$ by
integrating the dynamics according to Eq. (25) using the current approximation
$s_{0}^{(m)},\kappa^{(m)},I^{(m)}(\varphi)$. Numerically compute the integral
$S_{i}=\int_{\tau_{i}}^{\tau_{i+1}}(\psi(t)+s_{0})\,\mathrm{d}t$ with the
approximated amplitude. Solve system (23) with the better approximated
integral:

$s(\tau_{i+1})-s(\tau_{i})=-\kappa s_{0}(\tau_{i+1}-\tau_{i})+\kappa
S_{i}+\sum_{n=0}^{N_{F}}[u_{n}^{\cos}P_{n,i}^{\cos}+u_{n}^{\sin}P_{n,i}^{\sin}]$

to obtain the next approximate solution
$s_{0}^{(m+1)},\kappa^{(m+1)},I^{(m+1)}(\varphi)$. The deviation of the
reconstructed amplitude from the true amplitude estimated at events $\tau_{i}$
can be used as an error measure.

## References

  * Winfree (1980) A. T. Winfree,  _The Geometry of Biological Time_ (Springer, Berlin, 1980).
  * Mackey and Glass (1988) M. C. Mackey and L. Glass,  _From Clock to Chaos: The Rhythms of Life_ (Princeton Univ. Press, Princeton, NJ, 1988).
  * Rinzel and Ermentrout (1998) J. Rinzel and B. Ermentrout, in  _Methods of Neuronal Modeling_ , edited by C. Koch and I. Segev (MIT Press, Cambridge, 1998) pp. 251–292.
  * Canavier (2006) C. C. Canavier, Scholarpedia 1(12), 1332 (2006).
  * Achuthan and Canavier (2007) S. Achuthan and C. C. Canavier, J. Neurosci. 29, 5218 (2007).
  * Ermentrout (2002) B. Ermentrout,  _Simulating, Analyzing, and Animating Dynamical Systems: A Guide To Xppaut for Researchers and Students_ (Society for Industrial and Applied Mathematics, USA, 2002).
  * Tass (1999) P. A. Tass,  _Phase Resetting in Medicine and Biology. Stochastic Modelling and Data Analysis._ (Springer-Verlag, Berlin, 1999).
  * Tass (2000) P. A. Tass, Prog. Theor. Phys. Suppl. 139, 301–313 (2000).
  * Tass (2001a) P. A. Tass, Europhys Lett. 53, 15 (2001a).
  * Tass (2001b) P. Tass, Europhys. Lett. 55 (2001b).
  * Tass (2002) P. Tass, Phys. Rev. E 66, 036226 (2002).
  * Rosenblum and Pikovsky (2004a) M. G. Rosenblum and A. S. Pikovsky, Phys. Rev. Lett. 92, 114102 (2004a).
  * Rosenblum and Pikovsky (2004b) M. G. Rosenblum and A. S. Pikovsky, Phys. Rev. E. 70, 041904 (2004b).
  * Popovych  _et al._ (2005) O. Popovych, C. Hauptmann, and P. A. Tass, Phys. Rev. Lett. 94, 164102 (2005).
  * Tukhlina  _et al._ (2007) N. Tukhlina, M. Rosenblum, A. Pikovsky, and J. Kurths, Phys. Rev. E. 75, 011019 (2007).
  * Hauptmann and Tass (2009) C. Hauptmann and P. A. Tass, J Neural Eng. 6, 016004 (2009).
  * Wilson  _et al._ (2011) C. J. Wilson, B. Beverlin, and T. Netoff, Frontiers in Systems Neuroscience 5 (2011), 10.3389/fnsys.2011.00050.
  * Popovych and Tass (2012) O. V. Popovych and P. A. Tass, Front Hum Neurosci. 6, 58 (2012).
  * Lin  _et al._ (2013) W. Lin, Y. Pu, Y. Guo, and J. Kurths, EPL (Europhysics Letters) 102, 20003 (2013).
  * Zhou  _et al._ (2017) S. Zhou, P. Ji, Q. Zhou, J. Feng, J. Kurths, and W. Lin, New Journal of Physics 19, 083004 (2017).
  * Wilson and Moehlis (2016a) D. Wilson and J. Moehlis, PLOS Computational Biology 11, 1 (2016a).
  * Popovych  _et al._ (2017) O. Popovych, B. Lysyansky, M. Rosenblum, A. Pikovsky, and P. Tass, PLOS One 12, e0173363 (2017).
  * Krylov  _et al._ (2020) D. Krylov, D. Dylov, and M. Rosenblum, Chaos: An Interdisciplinary Journal of Nonlinear Science 30, 033126 (2020).
  * Montaseri  _et al._ (2013) G. Montaseri, M. Javad Yazdanpanah, A. Pikovsky, and M. Rosenblum, Chaos 23, 033122 (2013).
  * Holt  _et al._ (2016) A. Holt, D. Wilson, M. Shinn, J. Moehlis, and T. Netoff, PLoS Comput Biol. 12, e1005011 (2016).
  * Rosenblum (2020) M. Rosenblum, Chaos: An Interdisciplinary Journal of Nonlinear Science 30, 093131 (2020).
  * Duchet  _et al._ (2020) B. Duchet, G. Weerasinghe, H. Cagnan, P. Brown, C. Bick, and R. Bogacz, J. Math. Neurosci. 10, 4 (2020).
  * Benabid  _et al._ (1991) A. Benabid, P. Pollak, C. Gervason, D. Hoffmann, D. Gao, M. Hommel, J. Perret, and J. De Rougemont, Lancet 337, 403 (1991).
  * Benabid  _et al._ (2009) A. Benabid, S. Chabardes, J. Mitrofanis, and P. Pollak, Lancet Neurol. 8, 67 (2009).
  * Kühn and Volkmann (2017) A. Kühn and J. Volkmann, Mov. Disorders. 32, 11 (2017).
  * Holt and Netoff (2014) A. Holt and T. Netoff, J Comput Neurosci 37, 505 (2014).
  * Note (1) We highlight that in the context of electrical stimulation, e.g., in neuroscience applications, $p(t)$ is a current or a voltage and the introduced action has the physical meaning of electrical charge.
  * Note (2) Instantaneous phase and amplitude of a real-valued signal $s(t)$ are the argument and absolute value of the complex-valued analytical signal $s(t)+\mathrm{i}s_{H}(t)$, where $s_{H}$ is the Hilbert Transform of $s$, see, e.g., pikovsky2001; King-09; Feldman (2011).
  * Wilson and Moehlis (2016b) D. Wilson and J. Moehlis, Physical Review E 94, 052213 (2016b).
  * Wilson and Ermentrout (2018a) D. Wilson and B. Ermentrout, SIAM Journal on Applied Dynamical Systems 17, 2516 (2018a).
  * Wilson (2019) D. Wilson, Physical Review E 99, 022210 (2019).
  * Rosenblum and Pikovsky (2019) M. Rosenblum and A. Pikovsky, Phil. Trans. R. Soc. A 377, 20190093 (2019).
  * Note (3) Alternatively, one can detect the signal’s maxima; this approach corresponds to threshold-crossing for the derivative.
  * Note (4) The empirical curve is an $8^{\text{th}}$-order Fourier fit of experimental points.
  * Note (5) Note that the oscillation frequency $\omega$ shall be either determined by fit or set to be equal to the central frequency of the bandpass filter, see Holt and Netoff (2014); Holt  _et al._ (2016).
  * Note (6) In their experiments, Duchet et al. Duchet  _et al._ (2020) used a complex stimulus consisting of 25 bursts of high-frequency pulses, of the total length of 5 s.
  * Feldman (2011) M. Feldman,  _Hilbert Transform Applications in Mechanical Vibration_ (Wiley, NJ, USA, 2011).
  * Note (7) For PRC estimation, a possible approach would be computing the error measure $E_{Z}$ according to Eq. (15) as a function of $\delta_{\text{off}}$ and searching for the minimum. However, this approach does not apply to ARC estimation since Hilbert technique does not yield the Floquet exponent. While we can use the same $\delta_{\text{off}}$ value to infer both curves, there is no guarantee that the optimal value is the same.
  * Kralemann  _et al._ (2007) B. Kralemann, L. Cimponeriu, M. Rosenblum, A. Pikovsky, and R. Mrowka, Phys. Rev. E 76, 055201 (2007).
  * Kralemann  _et al._ (2008) B. Kralemann, L. Cimponeriu, M. Rosenblum, A. Pikovsky, and R. Mrowka, Phys. Rev. E 77, 066205 (2008).
  * Rosenblum  _et al._ (2021) M. Rosenblum, A. Pikovsky, A. Kühn, and J. Busch, Scientific Reports 11, 18037 (2021).
  * Cestnik and Rosenblum (2018) R. Cestnik and M. Rosenblum, Scientific Reports 8 (2018), 10.1038/s41598-018-32069-y.
  * Note (8) Each phase approximation $\varphi^{(m)}$ is used to compute the integrals that represent coefficients of the linear system: $\intop\ilimits@\mathop{cos}\nolimits(n\varphi^{(m)})p(t)\tmspace+{.1667em}\mathrm{d}t$, $\intop\ilimits@\mathop{sin}\nolimits(n\varphi^{(m)})p(t)\tmspace+{.1667em}\mathrm{d}t$. Then when solving the linear system and obtaining a better PRC approximation $Z^{(m+1)}(\varphi)$, it is used to recompute the next approximation of the phase $\varphi^{(m+1)}$ by integrating the Winfree Eq. (1).
  * Note (9) Note that the events $\tau_{i}$ used for isostable inference generally differ from those used for the PRC reconstruction, though we keep the same notation. Practically, to determine $\tau_{i}$ we use linear interpolation of the monotonically growing function $\varphi(t)$ given in a discrete set of points.
  * Note (10) Since the forcing term is just a function of time, the isostable dynamics can be integrated via variation of constant, in which case the time integral can be computed explicitly: $\psi^{(m)}(t)=(s(\tau_{i})-s_{0}^{(m)})e^{\kappa^{(m)}(t-\tau_{i})}+\intop\ilimits@_{\tau_{i}}^{t}I^{(m)}(\varphi(t^{\prime}))p(t^{\prime})e^{\kappa^{(m)}(t-t^{\prime})}\tmspace+{.1667em}\mathrm{d}t^{\prime}$.
  * Note (11) The optimal threshold corresponds to thresholding with a Poincaré section that is tangential to a local isochron, see Fig. 15. We find it by performing a direct search over a reasonable range of threshold values and choosing the one that corresponds to the lowest error value (15), as explained in the method paper Cestnik and Rosenblum (2018).
  * Schwabedal  _et al._ (2012) J. T. C. Schwabedal, A. Pikovsky, B. Kralemann, and M. Rosenblum, Phys. Rev. E 85, 026216 (2012).
  * Note (12) Ignoring the phase dependence of $c(\varphi)$ is similar to ignoring the information of local isochrons in the phase response, as most techniques do, e.g., the Hilbert-based one.
  * Wilson and Ermentrout (2018b) D. Wilson and B. Ermentrout, Journal of Mathematical Biology 76, 37 (2018b).