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

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 -- a task relevant to 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.


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 concentrate on 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) 1-4 , a crucial characteristic of a self-sustained oscillator. Knowing the PRC one describes the oscillator's phase dynamics via the Winfree equatioṅ where ϕ is the phase, ω is the natural frequency, Z(ϕ) 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 5 . For known oscillator equations, one can obtain PRC by solving the adjoint problem 6 . 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 2,4 . Indeed, if the pulse at phase ϕ is Dirac's delta function, then the induced phase shift equals precisely Z(ϕ). 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 living tissue should 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 [7][8][9][10][11] of the clinical technique known as deep brain stimulation (DBS) 12 . DBS aims to modulate the 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 9,11,13 . 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 a signal's instantaneous amplitude and phase 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 in terms of phases and isostable variables [14][15][16] 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 in 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 inefficient.
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 the tools for testing and comparing inference techniques. In Section III we compare methods and their performance, including the novel reconstruction of the phase-isostable dynamics which is described in Section III C. Section IV summarizes and discusses the results. The appendix section introduces the phase-isostables representation and presents the details of the test systems and inference techniques.
A. State of the art and problems to be solved Here, we specify the problems and introduce some notations.

Estimating phase response in experiments
The traditional approach relies on measuring the phase shift ∆ϕ evoked by a pulse applied at phase ϕ. However, in practice, the stimuli are not Dirac's delta pulses, so the inferred response function differs from Z(ϕ). Therefore, we consider general pulses P(t) of length δ, P(t) = 0 for t / ∈ [0, δ]. In particular, we are interested in the charge-balanced stimuli that additionally fulfill the condition δ 0 P(t) dt = 0 . ( For neuroscience applications, such stimuli are required to avoid charge accumulation in living tissue. We denote the response to non-Dirac stimuli as empirical PRC Z P (ϕ). Thus, the first problem is to relate Z(ϕ) and Z P (ϕ) for a given P(t). In the following, we infer the oscillator's response by perturbing it with the pulse train We specify the choice of stimulation times t k in the test examples below. The phase shift caused by a stimulus also depends on the stimulus amplitude and hence, has to be normalized. For unipolar stimuli, it is natural to normalize by the integral f = δ 0 P dt, commonly referred to as the stimulus' "action" 17 . For bipolar charge-balanced pulses the normalization is ambiguous; we choose f = 1 2 δ 0 |P| dt. Thus, Z P = ∆ϕ/f . The success of the standard technique is due to the presence of well-defined spikes, typical for neuronal oscillators. Consider now a system that exhibits nearly sinusoidal oscillation. Although no well-defined markers exist, thresholdcrossing events can also determine the points of a constant phase value 18 . 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. 9,13 . 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. 11 approached this problem using the Hilbert Transform: they estimated the instantaneous phase before and after the stimulus and, in this way, computed the PRC. 19 Unfortunately, neither of the techniques published in 9,11,13 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.

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 the 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 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 11 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. 11 in this way introduced the amplitude response curve, ARC.
Another approach to characterizing the amplitude variations exploits the notion of isostables 16,20,21 . In this framework, the Winfree Eq. (1) is complemented by an equation for the isostable variable ψ that can be interpreted as a deviation from the limit cycle:ψ = κψ + I(ϕ)p(t) . ( The real-valued parameter κ (Floquet exponent) determines the stability of the cycle and is negative for stable cycles. The function I(ϕ) quantifies the phase dependence of the stimulation's effect. We denote I(ϕ) as the isostable response curve, IRC. For an introduction to the phase-isostable representation, see Appendix A. We remind that Eqs. (1,3) represent the 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 some improvements. Next, we present the method for inference of the phase-isostable Eqs. (1,3) for the first time.

II. TOOLS FOR TESTING AND COMPARING INFERENCE TECHNIQUES
In this section, we first suggest a simple approach for relating the response curve to stimuli of arbitrary shape Z P (ϕ) to the infinitesimal PRC Z(ϕ). Then we design and present the test systems. And finally, we introduce the error measures we will be using.

A. Relation between empirical ZP and infinitesimal Z PRCs
From the infinitesimal PRC Z(ϕ), one can compute the empirical PRC Z P (ϕ) as a response to an arbitrary stimulus P. One simply has to evaluate the effective phase shift for the chosen stimulation P at every phase by integrating the Winfree phase equation (1). For example, let us focus on evaluating the phase response for stimuli P at a particular phase ϕ * . One first needs to compute the phase in time for the duration of the pulse by solving the Winfree equation (1) using initial condition ϕ(t = 0) = ϕ * . Then the empirical PRC at the chosen phase corresponds to the phase shift that is induced by P over the duration of the stimulation: where f is the action of the pulse. The inverse problem of how to determine the infinitesimal response Z(ϕ) when given the empirical one Z P (ϕ) and stimulation shape P, is in general, much harder. However, for pulses with small action f , the empirical response Z P (ϕ) is well approximated as the convolution of the infinitesimal curve Z(ϕ) and the pulse shape P(t). In Fourier space, convolution corresponds to multiplication, meaning the inverse operation is represented by division. One can, therefore, easily deconvolve an empirical response to the infinitesimal one by dividing the Fourier representations of Z P and 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 P → Z is generic and can be combined with any technique that infers an empirical response. B. Test models

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: where the parameter µ governs the limit cycle's stability as it is directly related to its Floquet exponent by κ = −2µ and α is the non-isochronicity parameter. The frequency of the limit-cycle oscillation is given by ω := η − αµ. σ is the noise strength, and ξ x,y are two realizations of the Gaussian white noise with zero mean and unit variance. p(t) is the external perturbation specified in the test examples below, and parameter β 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 √ µ. For definiteness, throughout this paper, we assume that the perturbation is state-independent in the chosen Cartesian coordinates (5). Then, the SL system (5) has PRC (see, e.g., 22 ) and IRC For the SL system, we also define the empirical amplitude response A(ϕ) = R e /R s , where R s,e are the values of the amplitude variable R immediately before and after the applied pulse. For a narrow unipolar pulse, this expression reads: where f is the pulse's action. For the 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 ω, the Floquet exponent κ, and non-isochronicity α as parameters. With perturbations, these equations read: where the functions C and D are: and r is a positive parameter determining the shape of the limit cycle in state space. System (9) has IRC and PRC r + 2 cos 2 (ϕ) − α (r + 1) cos(ϕ − β) + cos(3ϕ − β) (r + 2 cos 2 (ϕ)) 3 2 .
The system's dynamics are illustrated in Fig. 1.   (14). (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. 11 .

High-dimensional test model
Our third test system is a simple model of macroscopic neuronal oscillation. We take a system of N globally coupled Bonhoeffer-van der Pol oscillators:ẋ where k is the oscillator index, k = 1, . . . , N , and the term εX describes the mean-field coupling, where X = N −1 k x k . In the following, we take N = 1000. The oscillators' frequencies are determined by the parameter J k that is Gaussian-distributed with mean 0.6 and standard deviation 0.1. The coefficient ε explicitly describes the interaction between the ensemble elements. We choose ε = 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 k y k in (a) and X(t) in (b). With this test example, we imitate the natural variability of real-world signals.

C. 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 where · = (2π) −1 2π 0 (·) dϕ denotes the average value of a function, and Z rec is the inferred (recovered) PRC. The error value is L Z,A,I = 0 if the two curves coincide, it is of order 1 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 C 2 when introducing the method we also specify an analog for the isostable variable error. Note that the inferred phase response curve Z(ϕ) in conjunction with the natural frequency ω can be used in Eq. (1) to reproduce ϕ(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 τ 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 τ i we reproduce the phase evolution up to the time τ i+1 . In the ideal noise-free case, if the stimulation is weak and the inferred curve is exact, the reproduced phase ϕ(τ i+1 ) = Φ i equals 2π (or, equivalently, zero since we consider the wrapped phase). In practice, this reproduced phase Φ i deviates from 2π, 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 Φ i from 2π: where · 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 = τ i+1 − τ i (instantaneous periods): 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 our inferred model and the data itself. This makes it a helpful tool for experimental data analysis.

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(ϕ) and IRC I(ϕ) or ARC A(ϕ) 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 review the standard PRC inference technique, followed by testing the performance of other methods in use. Here we also describe and test our approach to the problem.

A. Standard approach (phase response only)
The standard technique is very efficient for neuronal oscillators exhibiting slow and fast motion. A spike corresponds to an epoch of fast motion, where the isochron density in the phase space is low. Therefore, since spike detection via 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. (9). Orange (full) and gray (dashed) bold transparent curves depict the theoretical curves Z(ϕ) and ZP (ϕ). Blue (dashed) and red (full) 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 = ϕ 2π T (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.
threshold-crossing is weakly dependent on the threshold, the precision of the phase estimation is high. Consider now a system without slow-fast motion that exhibits nearly sinusoidal oscillation. To trace the variation of the period, we determine states with equal phases employing threshold-crossing 23 , cf. 18 . 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 τ n of the n th threshold-crossing event and compute where τ 0 is the threshold-crossing event preceding the stimulus, and the phase of the stimulus application is ϕ = 2π(t s − τ 0 )/T where T is the natural period. The choice of n depends on the stability of the limit cycle, quantified with the Floquet exponent κ; as a rule of thumb, we suggest nκT 1. However, in practice, κ is unknown, and n shall be chosen by trial and error. To obtain the function Z P (ϕ), we repeat the perturbation P for different ϕ. To this end, we either choose random intervals between stimuli or stimulate periodically with the period incommensurate with T . 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. (5,9). For the SL system (5) the parameter values are: ω = 1, κ = −0.1, α = −0.3 and for the mSL system (9) they are ω = 1, κ = −0.1, α = 0, r = 0.75. There is no noise, σ = 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 = 1 2 |P| dt 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(ϕ) according to Eqs. (6) and (13), (ii) theoretical response to charge-balanced pulses Z P -the ground truth, (iii) inferred empirical PRC Z P obtained via Eq. (18) ("direct inference") 24 , and (iv) the result of the deconvolution of the empirical PRC Z P as described in Appendix B ("deconvolved"). For this simple case, we see a nearly perfect reconstruction of the effective PRC Z P (compare blue and gray curves in Fig. 3) as well as a successful deconvolution to the infinitesimal curve (compare orange and red). Since we used bipolar pulses, the empirical PRC resembles the negative derivative of the infinitesimal curve (compare gray with orange).
We emphasize that there is no analog to the standard technique for inferring the amplitude response.

B. 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 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 ϕ s,e and a s,e , respectively, where indices s and e stand for "start" and "end". From these quantities, we obtain the empirical PRC Z P (ϕ s ) = (ϕ e − ϕ s − ωδ)/f , where δ is the pulse's width, ω is the frequency of the unperturbed oscillation, and f is the normalization factor. Next, following Duchet et al. 11 , we introduce the empirical amplitude response curve (ARC) as A(ϕ) = a e /a s , to be distinguished from IRC. Below, we explore the performance of three different techniques.

Fitting the signal by a sine (phase response only)
Holt et al. 9,13 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 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 (5) in Fig. 4, computing the phase shift in the time domain. Suppose the pulse of length δ occurs at t s . Let the fitted functions be s(t) ≈ a 1 cos[ω(t − (t s − T fit )) + χ 1 ] for t s − T fit < t < t s and s(t) ≈ a 2 cos[ω(t − (t s − T fit )) + χ 2 ] for t s + δ < t < t s + δ + T fit . The phase of the second cosine at t = t s + δ is ϕ = ω(T fit + δ) + χ 2 . If there were no stimulus, the phase of the first cosine at this point would be ϕ = ω(δ + T fit ) + χ 1 . Hence, the phase shift is ∆ϕ = χ 2 − χ 1 . 25 Naturally, this technique applies to the signals close to 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 unipolar rectangular pulses with amplitude 0.1 and width δ = 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.

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). Duchet et al. 11 exploited HT to compute the phase and amplitude response of tremor oscillation with DBS of the thalamus. They took the instantaneous phase and amplitude immediately before and after the stimulus for this Here, the black dashed-dotted curve shows the evolution of the oscillator's phase ϕ, while the red curve illustrates the HT-based phase ϕH (for both phases, we subtract the average growth with the frequency ω). The vertical lines in all panels indicate the beginning and end of the rectangular pulse. In the absence of stimulation, the HT amplitude aH (t) reproduces R(t), but the discrepancy is essential before and after the pulse. For the phases, the effect is even more pronounced.
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, δ = 0.6). The parameters of the SL system are the same as before. Figure 5 illustrates the results, see also Fig. 3.5 in 27 . 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. 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 −δ off , t s ) and (t s +δ, t s +δ +δ off ), where t s is the instant of the stimulus's onset and δ off is the offset time. Next, we obtain the phase ϕ s by extrapolating the linear fit of ϕ(t) over the interval (t s −δ off −δ fit , t s −δ off ) to the instant t s ; the parameter δ fit is the length of the fitting interval. Note that fitting requires the unwrapped phase. Similarly, extrapolating the fit over the interval (t s + δ + δ off , t s + δ + δ off + δ fit ) to t s + δ yields ϕ e . Moreover, linear fit provides the frequency ω. Hence, we compute Z P (ϕ s ). The amplitude response (ARC) we compute as We emphasize that the algorithm has two parameters, δ off and δ fit .
We demonstrate the performance of the described algorithm perturbing the SL system (5) by rectangular pulses, taking the inter-pulse interval incommensurate with the natural period T -in this way, we ensure that stimulation occurs at different phases. We set the observable s(t) = x(t) and first choose β = 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(ϕ), A(ϕ) according to Eqs. (7,8) To compute the latter, we choose δ fit = 50δ and δ off = 2δ 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 δ fit , the choice of δ off is crucial. We demonstrate this in Fig. 6c,d by showing the dependence of the error of inference on δ off . We compute this error according to Eq. (15). We see, that a proper choice of the offset δ off essentially affects the inference; the reasonable results in Fig. 6a,b are due to the optimal value δ off /δ = 2. Unfortunately, we do not see any practical way to choose δ off when the true curve is unknown 28 . Moreover, the results strongly depend on how the stimulation enters the Panels (c,d) show the inference errors as a function of the offset parameter δ off for β = 0 (black circles), β = π/4 (red crosses), and β = π/2 (blue stars). An optimal choice of δ 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.
system's equations. To demonstrate this, we plot in Fig. 6c,d the corresponding curves for β = π/4 and β = π/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 11 , these inference results may be loosely related to the theory. Additionally, we mention that HT serves merely as a signal embedding technique, yielding the angle variable or protophase; see 29,30 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.

Estimating phase and amplitude using a virtual auxiliary oscillator
Finally, we adapted the technique 31 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ẍ + α a,ϕẋ + η 2 x = s(t). We choose the oscillators' frequencies η to be much larger than the characteristic frequency ν of s(t) so that the systems are far from resonance and the response weakly depends on ν. Next, we take the damping parameters α ϕ and α a for the phase and amplitude measurement, to ensure a simple relation between the phase ϕ(t) = arctan(−ẋ/νx) and amplitude x 2 + (ẋ/ν) 2 of the forced oscillation and those of the investigated signal. The implementation is simple and boils down to the 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 and 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 η 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/phase twice, first for the Next, we compute the amplitude ofŝ k via the same algorithm and flip it in time to obtain the flipped amplitude denoted byâ(t).â(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,φ(t) → −φ(t). We illustrate this algorithm in Fig. 7; parameters of the virtual oscillators are η = 5, α a = 6, α ϕ = 0.2.
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 in the angle variable, it generally suffers significant errors because it does not include information on local isochrons.
C. 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 32 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.

Inference of the phase response
The key step in the inference is determining time instants τ i corresponding to the same asymptotic phase ϕ; 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(τ i ) = s thr , choosing one crossing per period, e.g., with an additional condition d dt s(τ i ) > 0. The choice of the threshold s 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 (16).
The core of the inference is fitting the Winfree phase equation (1) integrated over individual periods, determined as the time interval between two threshold-crossing events, [τ i , τ i+1 ]: The left-hand side equals 2π due to the definition of events τ i having the same phase. On the right-hand side, we approximate the PRC by a finite Fourier series of order N F : Z(ϕ) = N F n=0 z cos n cos(nϕ) + z sin n sin(nϕ) . By interchanging the order of integration and summation, we obtain from Eqs. (19) a linear system for the unknown Fourier coefficients z n and frequency ω. Since we have as many equations (19) 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 cos(nϕ)p(t) dt, sin(nϕ)p(t) dt, however, cannot be evaluated yet because the phase ϕ(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 33 .
Approximating the phase initially as linearly growing between the events, ϕ (0) (t) = 2π t−τi τi+1−τi , we obtain the first-approximation solution of Eq. (19), namely ω (1) and Z (1) (ϕ). (The superscripts denote the iteration number.) Then, we obtain the next approximation of the phase by integrating Winfree Eq. (1) between events. In general, the m th -order approximation for the phase is obtained as: Estimated in this way, the phase at the end of a period generally differs from 2π: Φ we additionally re-scale the phase with the factor Φ (m) i to ensure ϕ (m) (τ i+1 ) ≡ 2π. As a result, the approximations gradually improve through iterations. Estimated phases at the end of periods Φ i indicate how well our inference fits the observations, see error measure (16). For each iteration (m) we can compute the error and monitor the convergence. For further details of the technique, we refer to Ref. 32 .

Inference of the isostable variable 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 have already inferred the PRC, we also have the instantaneous phase ϕ(t).
Like the PRC technique, the isostable inference relies on time events of equal phase τ i , and additionally, on estimating the isostable variable at those events, ψ i ≡ ψ(τ i ). The time events τ i are straightforwardly obtained from the instantaneous phase ϕ(t) 34 , 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 straight lines can approximate isochrons. Thus in the first-order approximation, the isostable variable linearly depends on signal s(t): where factors c and s 0 are generally phase-dependent. However, since all crossing events τ 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 ≡ 1. Note that approximation via Eq. (21) 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 [τ i , τ i+1 ]: This system of equations is similar to Eq. (19) in that it can be approximated with a linear system by expanding the unknown function in a Fourier series, I(ϕ) = N F n=0 u cos n cos(nϕ) + u sin n sin(nϕ) , 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 (21) we evaluate the left-hand side as ψ i+1 − ψ i = s(τ i+1 ) − s(τ i ) because the still unknown constant s 0 cancels in the difference. The integral ψ(t) dt is challenging since it includes a time-dependent isostable variable, which we do not have.
There are two issues with evaluating the isostable variable integral τi+1 τi ψ(t) dt: (i) we do not know the constant s 0 a priori, and (ii) we want to use approximation (21) only in discrete events τ 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: In this representation κs 0 becomes just another unknown variable that will be determined while solving the linear system (22), which together with determining κ allows expressing s 0 as the ratio s 0 = κs0 κ . The remaining integral involves the continuous quantity ψ(t)+s 0 , which coincides with our observable s(t) at events τ i (generally this quantity differs from the signal for t = τ i : s(t) = ψ(t) + s 0 ). The property s(τ i ) = ψ(τ i ) + s 0 also helps us resolve the second issue of only using discrete values s(τ i ) to estimate the integral of the continuous isostable variable. Just as in the phase response method, we approach this iteratively. First, we linearly interpolate ψ(t) + s 0 between known events so that we can evaluate the last term: 0 , κ (1) , I (1) (ϕ). Next, we use this solution to improve the estimate of the continuous isostable variable by integrating the underlying dynamical equation (3). In general, from the approximate solution of order m, we obtain the subsequent isostable variable estimation via integration 35 : Estimated in this way, the isostable variable ψ (m) starts the interval [τ i , τ i+1 ] as a direct estimation from the signal: 0 , but further along the interval it can deviate due to the approximation in integrands κ (m) , I (m) . Thus, at the end of the interval this quantity does not exactly correspond to s(τ i+1 ) − s (m) 0 , which is the value that we take to start the next interval [τ i+1 , τ i+2 ]. The estimated time series of the isostable variable obtained in this way is thus discontinuous at the events τ 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 isostable variable estimated at the end of a period as Ψ We define the error of the fit as the standard deviation of the difference: Just like the PRC error was compared to the irregularity of inter-event intervals, see Eq. (17), this error should be compared to the irregularity of the isostable variable at the events: The entire procedure of inferring the isostable variable response is condensed into step-by-step instructions in Appendix F.

Performance of the IPID-1 phase-isostable reconstruction
We test the performance of the introduced IPID-1 method on the two example oscillators (5) and (9). We infer the PRC Z(ϕ) and IRC I(ϕ) with the approaches explained in Sections III C 1 and III C 2 respectively. The results can be seen in Fig. 9. Parameters are the same as used in Sec. III A and Fig. 3. We stimulate with a bipolar charge-balanced pulse with a period of positive stimulation lasting 0.2, a 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 a relatively long time series (1500 periods), and there was no noise or other unknown inputs (see Appendix D for details). Additionally, we determined events of equal phase with the optimal threshold 36 by performing a direct search over possible threshold values, minimizing the error (16). 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 shown in a later Figure. 13.

D. 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 (5) and (9)  σ 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 unipolar 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 σ, for three different values of β, cf. Eq. (5). As expected, for this method the PRC inference is more sensitive to noise than the ARC reconstruction. Next, we proceed with testing the most promising IPID-1 technique, based on the phase-isostable dynamics reconstruction. We keep the forcing pulse action constant, and then for each σ generate several trajectory realizations, and for each one, we compute the L error measure using Eq. (15). 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 with the technique introduced in Section III C 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 marginally 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 a much shorter time series, see Appendix section D for details.
Finally, we test our phase-isostable reconstruction technique on the complicated case of the modified SL oscillator (9) that provides a non-sinusoidal observable, cf. Fig 1b. The results in Fig. 11c,d confirm the advantage of this technique, namely, it is not restricted to sinusoidal signals. Moreover, the inference errors L Z , L I for the modified SL are not much different from those for the standard SL that provides a sinusoidal signal. The results in panel (c) clearly demonstrate that our technique is more stable with respect to noise than the standard one.

E. 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. (14). 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 37 . We, therefore, estimate the goodness of inference solely based on how well the inferred model reproduces observations; namely, we consider the ratio of errors (16) and (17): 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 inferred the empirical PRCs using the mean field in the x variable as our observable, s = X, and exploiting bipolar pulses as before. Next, we deconvolved all curves to obtain the theoretical phase response curves. We also performed the inference on a different observable s = X + 2Y of the same system and obtained similar results.
Only three techniques can infer the amplitude response: the Hilbert-based, auxiliary oscillator-based, and IPID-1.  Figure 12. Application of all phase response methods on a high dimensional example of a chaotic oscillatory ensemble (14). (a) Inferred PRCs and (b) the ratio of errors (16) and (17), 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).
The first two methods yield the effective ARC. In contrast, the IPID-1technique estimates the infinitesimal isostable curve I(ϕ). Hence, for comparison, we recompute I(ϕ) into the ARC in the following way: first, using the infinitesimal curve I(ϕ) and pulse shape P(t) we evaluate the isostable shift ∆ψ. Then we relate it to the ARC. For a particular phase ϕ * , we integrate the first-order dynamics (1,3) with the considered pulse shape P(t) to obtain phase as a function of time, ϕ(t), by solving the Winfree equation with the initial condition ϕ(t = 0) = ϕ * , and then compute the effective isostable shift as: ∆ψ(ϕ * ) = δ 0 (κψ(t) + I(ϕ(t))P(t)) dt .
Next, we have to relate the isostable variable ψ to amplitude a. As already mentioned, locally, the isostable linearly depends on the distance from the limit cycle, see Eq. (21): ψ ≈ c(ϕ)(R − R 0 (ϕ)). Here R represents the distance from the origin, while c(ϕ) and R 0 (ϕ) 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 (ϕ) + ψ/c(ϕ).
Recalling the definition of ARC as a ratio of amplitudes before and after a pulse, we write: We cannot estimate c(ϕ) nor R 0 (ϕ) 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 max 38 . 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(ϕ) 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 by quantifying how well the dynamical model fits the observation. Namely, we compute the ratio of errors (25) and (26). 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 the 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 ψ 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 an 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.

IV. DISCUSSION AND CONCLUSIONS
This paper addressed the problem of inferring an oscillator's phase and amplitude response from observations using test stimulation. Below, we discuss these two tasks of phase and amplitude inference separately.
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 concentrated on the case where we apply a specially designed stimulation. However, this technique can be exploited for almost any input, as long as we can observe it and it does not entrain the oscillator 39 .
A. 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 P (ϕ) → Z(ϕ). This technique can be used alongside 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 error estimation. Without knowing the ground truth and using only the observations, we quantify the inferred PRC's accuracy in reproducing the observations. 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 0 π/2 π 3π/2 2π  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 a 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.

B. 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, which provides a universal description independent of 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 20 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 remain a subject of further studies.

C. 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 8,10 . The 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 11 . 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 11 , generally, this relationship does not hold, see Appendix A 3 for details. Another idea is stimulating around zeros of Z(ϕ). 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 the 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. N -dimensional system can be represented in the isostable coordinate systeṁ where ω is the frequency, and κ i are Floquet exponents of the limit cycle. For a system with a stable limit cycle, all Floquet exponents have a negative real part, (κ i ) < 0. However, in many cases, it is sufficient to consider only one isostable variable ψ 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 ψ i decay considerably faster than ψ max . Restricting our consideration to only one isostable coordinate ψ max complementing ϕ, 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 κ := max i (κ i ) and below we omit the subscript max in the notation ψ 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 16,20,21 .
As one can check from the autonomous dynamical Eqs. (A2), all variables ψ i are invariant to scaling with some non-zero factor c -the dynamical equations remain unchanged after a change of variablesψ := cψ. As soon as external perturbations p(t) are introduced to the system viȧ the scaling of isostable variable leads to the adjustment in the stimulation function Thus, the IRC (defined as the isostable amplitude response evaluated at the limit cycle) is also scaled as a result of the scaling of ψ: The PRC is invariant under rescaling of the isostable ψ. In the following derivations, we keep the scaling factor c to show where it appears in the equations.

Isostable reduction for the SL system
The Stuart-Landau system, Eq. (5), is the analytically solvable normal form of a Hopf bifurcation: For µ > 0, a radial limit cycle at R 0 = √ µ 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 θ, can be assigned a phase by the well-known formula for which one can verify thatφ = η − αµ =: ω. Since the SL system is two-dimensional, one additional isostable ψ 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 R \ {0}. For the case of radial isochrons (α = 0), ψ has already been derived in Ref. 15 as Moreover, this coordinate transformation yieldsψ = −2µψ := κψ also in the case of non-radial isochrons (α = 0). From this expression, one can check the system's dynamical properties: the roots of ψ at R = √ µ 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 ψ converges to a finite value as R → ∞, 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 ω and κ instead of η and µ aṡ 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 ψ and ϕ allows us to compute the PRC Z(ϕ) and IRC I(ϕ) also in a closed analytical form. Those are, by definition, the response curves of phase and isostable variables evaluated at the limit cycle, where ψ = 0. In general, the response curves are computed via where J isostable←polar is the Jacobian of the coordinate transformation from polar to isostable coordinates, and J polar←Cartesian is the well-known Jacobian of the coordinate transformation from Cartesian to polar coordinates: Throughout this paper, we assume a state-independent stimulation, namely, G x = cos(β) and G y = sin(β). Hence we arrive at the following phase and isostable response curves for the SL system: 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 approximating a short pulse with action f by a Dirac delta function, hitting the system at phase ϕ. For the SL system, the initial amplitude R s = √ µ is independent of ϕ. For a stimulation that is state-independent in Cartesian coordinates, the response in R is given by G R (ϕ) = cos(ϕ − β) and thus the ARC simplifies to the expression

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 ϕ and ψ in the first place. Since the roots of ψ(x) are the location of the stable invariant set for κ < 0, we can control the position and shape of the limit cycle. In particular, the coordinate transformation describes a limit cycle at R 0 = q(θ) for any strictly positive, 2π-periodic function q of θ. Together with the phase defined as such that ϕ = θ on the limit cycle, we obtain a generalized Stuart-Landau system. For q(θ) = − κ 2 = µ we recover the original SL system. Since the Jacobian J isostable←polar can be computed in the same straightforward manner as for the SL system as we can explicitly write the dynamical equations, both in polaṙ and in Cartesian coordinates (where q and q have to be read as functions of θ(x, y)) Having the explicit Jacobians for the coordinate transformation, we obtain the PRC and IRC generally given by the response curves in polar angle G θ and polar radius G R as with the function I defined as The polar angle θ serves as a protophase in this model since, on the limit cycle, θ and ϕ coincide despite not being equal in general. Thus, on the limit cycle where ψ = 0, we have θ = ϕ and we obtain the shape of the limit cycle R 0 = q(ϕ) as a function of phase ϕ. Note the relation from Eqs. (A24), stating that the IRC is proportional to the difference of PRC and the response curve of the polar angle G θ ψ=0 for non-zero non-isochronicity parameter α: We emphasize this relation between PRC and IRC since it exists without further assumptions made about G θ 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 ϕ and polar angle θ 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 and I as For the particular example in Sec. II, we use q(θ) = r + 2 cos 2 (θ), 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 r + 2 cos 2 (ϕ) (A29) and the function I as I(ϕ) = (r + 1) cos(ϕ − β) + cos(3ϕ − β) (r + 2 cos 2 (ϕ)) 3 2 . (A30) Having these two function, the PRC in Eq. (13) and IRC in Eq. (12) are derived using the relations given in Eq. (A24) for c = 1.

Relationship of phase and amplitude response curves
Duchet et al. 11 hypothesized that the amplitude response is proportional to ∂ ϕ Z(ϕ). 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 = ρ cos(β) and G y = ρ sin(β), where the angle β determines the direction of the stimulation in the x-y-plane and the non-negative value ρ determines the strength. Without loss of generality, we can set ρ = 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. (A11), the response functions in the polar angle θ and the polar radius R follow as These two functions establish several relations, e.g., Thus, we obtain the relation that G R is proportional to the derivative of G θ . 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 ϕ by the polar angle θ 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 ∝ G R ∝ A − 1. In addition, for the isochronous SL system (α = 0, Fig. 14a) the phase ϕ equals the polar angle θ and the geometric relation (A32) directly translates into a relation between PRC and IRC: I ∝ ∂ ϕ Z.
For the non-isochronous SL system illustrated in Fig. 14b, this relation is no longer strict since phase ϕ and polar angle θ differ outside the limit cycle. Consequently, Z = G θ and ∂ ϕ 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. where F and 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. (B1) Z by Z + C, where C is an arbitrary constant, we obtain the same function Z P due to the charge balance condition te ts P dt = 0. To solve this problem, we set the mean value of Z equal to that of Z P . An alternative solution for the inverse problem that does not assume ϕ ≈ ω(t − t s ) is as follows. We represent the yet unknown function Z(ϕ) as a finite Fourier series with coefficients z cos n , z sin n . Then using P(t) we express Z P (ϕ) in terms of z cos n , z sin n by numerically solving the Winfree equation. The nonlinear problem of equating the empirical curve with its symbolic representation: Z P = Z P (ϕ; z cos n , z sin n ), is then solved using the Levenberg-Marquardt technique. The isostable structure of a 2D limit-cycle oscillator. The limit cycle is depicted in black, and the isochrons are shown in yellow. The two thresholds tangential to isochrons are shown in red and blue. (b) The corresponding signal x(t) is depicted in 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 F: Isostable inference step-by-step summary 1. Determine events τ i by thresholding the phase, e.g., ϕ(τ i ) = ϕ thr , d dt s(τ i ) > 0. The choice of threshold ϕ thr is not crucial but generally should be picked such that the signal in events s(τ i ), has a wide range of values. where we denoted the Fourier integrals with P cos n,i = τi+1 τi cos(nϕ(t))p(t) dt and P sin n,i = τi+1 τi sin(nϕ(t))p(t) dt. Numerically compute the integrals for each inter-event interval.

Write system (22) by linearly approximating the isostable integral
4. Combine the approximations to express system (22) as a linear system with known coefficients: [u cos n P cos n,i + u sin n P sin n,i ] .
There are as many equations as there are inter-event intervals. The known coefficients are: s(τ i+1 ) − s(τ i ), τ i+1 − τ i , (s(τ i ) + s(τ i+1 ))/2, and the Fourier integrals P n,i . The unknown quantities are: κs 0 , κ and the Fourier modes u cos n , u sin n , representing the response I(ϕ). Thus proceed to minimize this system via least squares (or a similar method) to obtain the first approximated solution: s (ψ(t) + s 0 ) dt with the approximated amplitude. Solve system (22) with the better approximated integral: [u cos n P cos n,i + u sin n P sin n,i ] to obtain the next approximate solution s (m+1) 0 , κ (m+1) , I (m+1) (ϕ). The deviation of the reconstructed amplitude from the true amplitude estimated at events τ i can be used as an error measure.