# Numerical phase reduction beyond the first order approximation

###### Abstract

We develop a numerical approach to reconstruct the phase dynamics of driven or coupled self-sustained oscillators. Employing a simple algorithm for computation of the phase of a perturbed system, we construct numerically the equation for the evolution of the phase. Our simulations demonstrate that the description of the dynamics solely by phase variables can be valid for rather strong coupling strengths and large deviations from the limit cycle. Coupling functions depend crucially on the coupling and are generally non-decomposable in phase response and forcing terms. We also discuss limitations of the approach.

###### pacs:

05.45.Xt Synchronization; coupled oscillators, phase dynamicsIt is widely accepted that the dynamics of weakly interacting limit cycle oscillators can be fully described by their phases, while the amplitudes can be considered as enslaved and thus irrelevant. This idea is commonly used both in theoretical and experimental studies. While in the latter case the phase dynamics is inferred from data, in the former one, the corresponding equation shall be derived by a perturbation technique. However, the existing theory provides only a recipe for computation of the first order approximation. In spite of numerous recent efforts, phase dynamics description beyond the limit case of weak coupling is still lacking. Here we provide a numerical approach for the phase dynamics reduction beyond the first order approximation.

## I Introduction

Phase dynamics approximation plays an important role in the analysis of coupled self-sustained oscillators Winfree-80 (); Kuramoto-84 (); Hoppensteadt-Izhikevich-97 (); Pikovsky-Rosenblum-Kurths-01 (); Ermentrout-Terman-10 (); Monga_Wilson-Matchen-Moehlis-18 (). Reduction to the phase description essentially decreases the dimension of the problem and in some cases provides an analytical solution. Well-known examples of the solvable phase dynamics are the Adler equation Adler-46 () for a periodically driven oscillator and the Kuramoto model of (infinitely) many globally coupled systems Kuramoto-84 (). However, a practical application of the theory and a derivation of the phase dynamics equation from equations formulated in state variables for an arbitrary driven or coupled oscillator remains an unsolved problem: this derivation requires an analytical expression for isochrons Guckenheimer-75 () or for the phase response curves, and therefore can be accomplished analytically only in exceptional cases, e.g., for the Stuart-Landau system. Moreover, this equation can be obtained only in the first approximation, as the first-order term of the perturbative expansion with respect to small coupling parameter Kuramoto-84 (). In spite of recent attempts (see a discussion in Refs. Monga_Wilson-Matchen-Moehlis-18, ; Wilson-Ermentrout-18, ), there exist no practical algorithms for inferring phase description in a general case when a perturbed trajectory is not very close to the unperturbed limit cycle.

In this communication we propose an efficient numerical approach which provides the desired phase dynamics equation. The evolution of the phase is described by a coupling function that is either defined on a fine grid or represented by a finite Fourier series. Next, we use this approach to study the phase dynamics for strongly driven systems, where we cannot expect validity of the first order approximation. For two different examples we reconstruct the second- and the third-order terms of the coupling functions, and verify that this description works with a good precision.

We start by introducing our notations and recalling main ideas of the phase reduction theory. Let the autonomous system be described by , where is the state vector in -dimensional phase space, . This equation defines the flow . Our goal is to find such a function of the phase space , mapping the state of the system onto a unit circle, that allows to attribute to every trajectory the corresponding phase evolution . Depending on the context, we will either treat as an wrapped phase (points on a unit circle), or as unwrapped phase with .

If the system possesses a stable limit cycle with period , then the phase definition can be constructed for all states in the basin of attraction of . First, the phase is easily introduced on . One assigns zero value of to an arbitrary point on the limit cycle. Then, for , the phase on the limit cycle is defined as , with .

Extension of this definition to the whole basin of attraction of the cycle is based on the isochrons Guckenheimer-75 (). For a point in the basin of attraction of , the phase is defined as , where (); obviously this limiting point belongs to the limit cycle, . In other words, isochrons, as the sets of equal phases, are stable manifolds of the corresponding points on the limit cycle under the action of the stroboscopic map .

Knowledge of isochrons, i.e. of the relation in a vicinity of the limit cycle, allows one to obtain phase equations for a perturbed system

(1) |

where describes external forcing or coupling with the strength . Then, for small , exploiting the perturbation approach, using the unperturbed definition of the phase, and neglecting deviations from the limit cycle (see Kuramoto-84 (); Pikovsky-Rosenblum-Kurths-01 () for details), one obtains

(2) |

where is called the coupling function and the subscript emphasizes that the latter is obtained in the first order approximation in . In typical cases, e.g., for periodic external forcing or coupling to another limit cycle oscillator, can be parameterized by its phase, (it means that we can write ), and then Eq. (2) can be written as

(3) |

Generally, the perturbation approach implies that the dynamics can be represented by a power series in :

(4) |

One can expect that this description is valid not only for small coupling, but as long as the dynamics of the driven system occurs on a smooth attractive torus. Although the methods for computation of the high-order terms are not known, the function , and the corresponding contributions , , can be obtained numerically, as is demonstrated below.

## Ii Numerical computation of coupling functions

### ii.1 Inference of the coupling function

The key issue is determination of the phase for any state of the perturbed oscillator.
Suppose that for some instant of time the perturbed oscillator
(1) has the state
and our goal is to compute .
Recall that the theory operates with the phases, defined for the autonomous system.
Thus, the problem reduces to determination of the phase of the latter.
As a preparatory step we choose a zero-phase point
on the limit cycle and compute period of the oscillation
^{1}^{1}1We use the known Henon-trick Henon-82 () for precise
determination of the times when the trajectory returns to ..
Next, we consider a copy of the autonomous part of Eq. (1),
i.e. we integrate equations , taking
for the initial condition.
The integration time is chosen to be ,
where an integer shall be taken sufficiently large so that we can assume
that the trajectory is already attracted to the limit cycle. (In fact, is the
only parameter of our procedure, it crucially depends on the second largest
multiplier of the limit cycle.)
Finally, we just have to compute the phase of the
point . For this goal we determine the time required for the
trajectory started from to achieve .
Then
.

Having introduced this crucial step of phase determination, we now present
the algorithm for coupling function inference.
We integrate the equations
of the perturbed system with some discrete time step and, after a sufficiently
large transient, obtain for each point
its phase as described above.
Next, we evaluate numerically .
Finally, we exploit the technique used in data analysis for reconstruction
of coupling functions from observations
^{2}^{2}2Practically, we use the Savitzky-Golay smoothing filter
for numerical differentiation of
and kernel density function estimation Chen_KernelTutorium-17 () for the fit,
see Kralemann_et_al-13 (); Rosenblum-Pikovsky-18 () for details. .
Namely, having the time series of ,
, and , with known frequency , we fit the r.h.s. of
Eq. (4) and thus obtain on the torus .
Practically, we obtain on a grid ; alternatively it can be
obtained as a finite double Fourier series.
Notice that the fit is possible if the
systems remains in a quasiperiodic state so that the trajectory densely fills
the torus.

### ii.2 Coupling function in the Winfree form

If only one component of the force in (1) is present, and this force does not depend on the state of the system, i.e., , the coupling function in the first approximation, according to (2), can be represented as a product

(5) |

where is the phase sensitivity function, also called the phase response curve (PRC). For a known force , determination of the coupling function thus reduces to finding PRC. This can be accomplished by expanding PRC in a Fourier series and obtaining the coefficients of this series via linear fit.

### ii.3 Decomposition in powers of the coupling strength

Here we shortly discuss how a general coupling function can be decomposed in a power series of , cf. Eq. (4). In our simulations we restrict the order of the series to three; hence, we are looking for a third-degree polynomial representation of and determine . For this goal it suffices to compute for at least three different values of coupling . Practically, in our examples shown below, we use more than 10 values of . Then, with account of the condition , the polynomial -dependence can be found by least mean squares fit. Notice that this procedure can be performed point-wisely for coupling functions given on a grid (below we use this option), or for each Fourier harmonics for the Fourier representation of . As a result, the partial coupling functions defined in (4) are obtained.

## Iii Results

In order to illustrate the advantages and limitations of our approach we analyze the phase dynamics of two simple paradigmatic models.

### iii.1 The Rayleigh oscillator

In the first example we consider the harmonically forced Rayleigh oscillator

(6) |

The nonlinearity parameter is ; for this value the limit cycle is strongly stable (its transversal Lyapunov exponent is , what for the period of the cycle yields a multiplier .) Notice that in this particular case the perturbation is scalar and therefore the first-order coupling function can be written in the Winfree form (5) with .

We fix the frequency of the forcing at and increase the amplitude from a small value (for which we can expect that the first-order phase approximation is valid) up to , which is not small at all: this forcing is comparable to the amplitude of . Notice that for larger values of the system gets locked to the force and determination of the coupling function with our method becomes impossible.

First, we demonstrate that the phase description is well-defined for the whole explored range of . For this purpose we compute, for all available , , , the rest term of our full model,

(7) |

and quantify the quality of the fit by

(8) |

where and bar denotes the time averaging over the available time series. The dependence is shown in Fig. 1a, together with a similar measure for the Winfree form (5). The latter increases with coupling strength, what clearly demonstrates that the Winfree form is valid for small couplings only. In contradistinction, the full coupling function describes the dynamics with accuracy in the whole range of couplings.

In Fig. 1(b,c) we also show the inferred coupling functions for weak, (here nearly coincides with ), and strong, , forcing amplitudes. Although qualitatively similar, the coupling functions are clearly different, indicating importance of higher-order components and .

Next, we obtain functions as discussed in Section II.3 and verify the validity of the series representation by Eq. (4). To this end we compute the error of the first-, second-, and third-order approximations

(9) |

where and . Here the averaging is performed over the torus on which the coupling function is defined:

These errors are shown in Fig. 2, while the functions are given in Fig. 3.

Fig. 2 shows that, while the first order approximation is quite accurate for , for larger coupling strengths the higher order terms become relevant. The third-order approximation provides a rather good accuracy (better than 1%) in the whole explored range of coupling strengths. Remarkably, the coupling functions are very different. While the first-order coupling function reproduces the known theoretical result (5), the second-order term contains strong second harmonics in both phases. Preliminary computations also show that is independent of the driving frequency , as expected, while the higher-order functions depend on ; this issue shall be analyzed in details separately.

### iii.2 The Rössler oscillator

For the second example we take a three-dimensional system, namely the harmonically driven Rössler oscillator:

(10) |

For the chosen parameter values, the unforced system has a limit cycle with complex multiplicators . Thus, the cycle is weakly stable; moreover, complex-valued multipliers mean that the dynamics cannot be reduced to a two-dimensional one. However, our approach successfully constructs the phase dynamics model. To demonstrate this, we infer the coupling functions for and forcing strengths varied from to . The quality of the model for different is illustrated in Fig. 4, where we show the normalized standard deviation , defined according to (8). One can see that the model remains valid for the coupling values as large as . We emphasize, that for such a driving the deviation of the trajectory from the unperturbed limit cycle is not small at all, see Fig. 5, where the attracting torus at this regime is depicted. For the coupling functions for weak and strong coupling see Fig. 6. Qualitatively, the results are similar to that for the Rayleigh oscillator: The coupling function is strongly nonlinear, as one can conclude from the comparison of two panels of Fig. 6.

We now discuss, why the technique fails for larger coupling, . In the case of the Rayleigh oscillator, the value of was limited by the effect of locking to the drive, what makes the inference of the coupling function impossible. In the case of the Rössler system, the reason for failure is different. Here, the torus of the driven system in the phase space becomes strongly shifted with respect to the original cycle, so that the trajectory starts to cross ”wrong” isochrons, see Fig. 7. In other words, the isochrons do not anymore represent a proper phase coordinate on the attracting torus. This results in a non-monotonic growth of the isochron-based phase, manifested by the violation of the “good-phase-condition” , cf. Fig. 4b.

## Iv Discussion

In summary, we have developed a simple and efficient technique for numerical determination of the phase of a perturbed self-sustained oscillator from the equations formulated in state space variables. In its turn, computation of the phase provides an easy way to reconstruct the equation (4) for the phase dynamics, where the coupling function is given on an equidistant grid within or/and can be given by a double Fourier series in phases . For transparency of the presentation we have considered driven systems, but extension to the case of two coupled oscillators is straightforward. The phases can be in the same way obtained for a network of more than two oscillators as well, but reconstruction of high-dimensional coupling functions requires extremely long data sets and becomes unfeasible, see a discussion in Ref. Kralemann-Pikovsky-Rosenblum-14, .

We have demonstrated that reduction to the phase description is valid even for quite strong driving/coupling and have confirmed numerically that the coupling function can be represented as a power series of the driving strength. The suggested algorithm for correct computation of the phase offers also an easy way to obtain the PRC of the system. However, our simulations show that description in terms of an (amplitude-dependent) PRC fails for strong coupling: one cannot generalize the Winfree-type coupling (5) to the full coupling function as .

We have illustrated two effects that limit applicability of the approach. The first one is synchronization of the system to the external perturbation, as in the case of the Rayleigh oscillator. In the example of the Roessler system, the technique does not work if a strong perturbation makes the mapping not unique within one cycle. Finally, we expect that the approach can also fail if a very strong coupling destroys the smooth torus in the phase space, see (Afraimovich-Shilnikov-83, ) for a rigorous treatment and Pikovsky-Rosenblum-Kurths-01 () for a qualitative discussion of this mechanism.

Our results open a way for further studies on phase reduction, e.g. they provide a numerical tool for an analysis of the dependence of high-order terms of the coupling functions on parameters of the driving, with the utmost goal to extend the analytical techniques. Finally, our approach can be used as a benchmark for testing technique for coupling functions reconstruction from data Kralemann_et_al-07 (); Kralemann_et_al-08 (); Kralemann_et_al-13 (); PhysRevLett.109.024101 () where only one observable of the system is available.

###### Acknowledgements.

The work was supported by ITN COSMOS (funded by the European Union Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No. 642563). Numerical part of this work was supported by the Russian Science Foundation (Rayleigh model: Grant No. 17-12-01534; Roessler model: Grant No. 14-12-00811).## References

- (1) A. T. Winfree, The Geometry of Biological Time (Springer, Berlin, 1980).
- (2) Y. Kuramoto, Chemical Oscillations, Waves and Turbulence (Springer, Berlin, 1984).
- (3) F. C. Hoppensteadt and E. M. Izhikevich, Weakly Connected Neural Networks (Springer, New York, 1980).
- (4) A. Pikovsky, M. Rosenblum, and J. Kurths, Synchronization. A Universal Concept in Nonlinear Sciences. (Cambridge University Press, Cambridge, 2001).
- (5) G. B. Ermentrout and D. H. Terman, Mathematical Foundations of Neuroscience (Springer, New York, 2010).
- (6) B. Monga, D. Wilson, T. Matchen, and J. Moehlis, “Phase reduction and phase-based optimal control for biological systems: a tutorial,” Biological Cybernetics(2018).
- (7) R. Adler, “A study of locking phenomena in oscillators,” Proc. IRE 34, 351–357 (1946), reprinted in Proc. IEEE, 61(10):1380â1385, 1973.
- (8) J. Guckenheimer, “Isochrons and phaseless sets,” Journal of Mathematical Biology 1, 259–273 (1975).
- (9) D. Wilson and B. Ermentrout, “Greater accuracy and broadened applicability of phase reduction using isostable coordinates,” Biological Cybernetics 76, 37–66 (2018).
- (10) We use the known Henon-trick Henon-82 () for precise determination of the times when the trajectory returns to .
- (11) Practically, we use the Savitzky-Golay smoothing filter for numerical differentiation of and kernel density function estimation Chen_KernelTutorium-17 () for the fit, see Kralemann_et_al-13 (); Rosenblum-Pikovsky-18 () for details.
- (12) B. Kralemann, A. Pikovsky, and M. Rosenblum, “Reconstructing effective phase connectivity of oscillator networks from observations,” New Journal of Physics 16, 085013 (2014).
- (13) V. S. Afraimovich and L. P. Shilnikov, “Invariant two-dimensional tori, their destroying and stochasticity,” in Methods of Qualitative Theory of Differential Equations (Gorki, 1983) p. 3â28, in Russian; English translation Amer. Math. Soc. Transl. Ser. 2 149:201â212, 1991.
- (14) B. Kralemann, L. Cimponeriu, M. Rosenblum, A. Pikovsky, and R. Mrowka, “Uncovering interaction of coupled oscillators from data,” Phys. Rev. E 76, 055201 (2007).
- (15) B. Kralemann, L. Cimponeriu, M. Rosenblum, A. Pikovsky, and R. Mrowka, “Phase dynamics of coupled oscillators reconstructed from data,” Phys. Rev. E 77, 066205 (2008).
- (16) B. Kralemann, M. Frühwirth, A. Pikovsky, M. Rosenblum, T. Kenner, J. Schaefer, and M. Moser, “In vivo cardiac phase response curve elucidates human respiratory heart rate variability,” Nature Communications 4, 2418 (2013).
- (17) T. Stankovski, A. Duggento, P. V. E. McClintock, and A. Stefanovska, ‘‘Inference of time-evolving coupled dynamical systems in the presence of noise,” Phys. Rev. Lett. 109, 024101 (2012).
- (18) M. Henon, “On the numerical computation of Poincaré maps,” Physica D: Nonlinear Phenomena 5, 412–414 (1982).
- (19) Yen-Chi Chen, “A tutorial on kernel density estimation and recent advances,” Biostatistics & Epidemiology 1, 161–187 (2017).
- (20) M. Rosenblum and A. Pikovsky, “Efficient determination of synchronization domains from observations of asynchronous dynamics,” Chaos 28, 106301 (2018).