High-Order Phase Reduction for Coupled Oscillators

We explore the phase reduction in networks of coupled oscillators in the higher orders of the coupling parameter. For coupled Stuart-Landau oscillators, where the phase can be introduced explicitly, we develop an analytic perturbation procedure to allow for the obtaining of the higher-order approximation explicitly. We demonstrate this by deriving the second-order phase equations for a network of three Stuart-Landau oscillators. For systems where explicit expressions of the phase are not available, we present a numerical procedure that constructs the phase dynamics equations for a small network of coupled units. We apply this approach to a network of three van der Pol oscillators and reveal components in the coupling with different scaling in the interaction strength.

One of the most famous theoretical tools for the analysis of coupled oscillators is the phase reduction [15][16][17][18].This approach provides a recipe for a description of complex high-dimensional oscillators by a single cyclic variable, phase ϕ, so that the dynamics of a network of N generally high-dimensional elements reduces to a set of N coupled one-dimensional differential equations for the phases.This reduction is based on the slaving principle (the corresponding notion in mathematical literature is the normal hyperbolicity): While the phases correspond to the neutral directions with zero Lyapunov exponents, the "amplitudes" (i.e.all other variables except for the phases) are stable and thus follow the phase dynamics.As a result, the full system's dynamics reduces to that on the N -dimensional torus spanned by the phases.
This reduction allows studying general behaviours of coupled oscillators through phase dynamics models, a prominent example here is the analytically tractable Kuramoto model of all-to-all interconnected units, derived for a population of weakly coupled oscillators close to the Hopf bifurcation point [15,19,20].Generally, phase dynamics models are expected to be valid for arbitrary oscillators and weak to moderate coupling.However, the existing theory employs a perturbative approach and provides phase equations only in the first-order approximation in the coupling strength.Derivation of high-order corrections would extend the validity of the approach beyond the weak-coupling limit and, thus, essentially increase our ability to analyse complex networks.Yet, for the moment it remains a theoretical challenge, in spite of numerous attempts [21][22][23][24][25][26].
In addition to obvious advantages for theoretical studies, phase reduction provides a framework for model reconstruction from measurements and, hence, plays an essential role in experimental research.Reconstruction of the phase dynamics equations from data is much more simple than recovery of equations for the state variables because (i) phase equations are low-dimensional and (ii) have a simple universal form because their right-hand sides are 2π-periodic functions of the phases.Though precise estimation of the phases from scalar signals remains a topic of ongoing studies [27][28][29], phase models have been successfully recovered for several laboratory experiments as well as for physiological systems [30][31][32][33][34][35][36][37][38].In particular, this approach yields a way to tackle the connectivity problem, i.e. to recover directional causal links between oscillatory sources solely from multivariate data.
However, for more than two oscillators and moderate coupling the connectivity, as it appears in phase equations, generally differs from the connectivity as defined by physical connections, e.g. the units that are not directly coupled appear as interacting with each other or the interaction of pairwise physically connected elements appears as nonpairwise in the phase description.These effects are due to the terms that do not appear in the first-order phase approximation.Thus, reconstruction and interpretation of phase models requires better understanding of the high-order phase reduction [8,39,40].
In this paper, we take a step to explore high-order phase dynamics.Our main goal is to match theoretical derivation of higher-order phase coupling terms with the data analysis.After a brief introduction to phase reduction approach in Section 2.1, in the rest of Section 2 we outline the perturbation approach and apply it to a particular example of three coupled Stuart-Landau oscillators.Our technique is close to the procedure of Ref. [21], but is better suited for generic networks and different natural frequencies of the oscillators.For this system, we derive and discuss phase dynamics equations in the higher orders in the coupling parameter (we present the second-order terms in full details, while for the higher orders we discuss only the phase dependence of the terms).Next, in Section 3 we present a numerical framework for the reconstruction of the phase equations from the simulations.First, we apply this framework to the Stuart-Landau network and in this way verify our analytical results.Second, we analyse a network of three van der Pol oscillators for which only a numerical study is currently possible, and identify the scaling of different coupling terms with the coupling strength parameter.In Section 4 we discuss our results.

Coupled oscillators and first-order phase reduction
We first briefly recall the basic ideas of the phase reduction in the first order in the coupling parameter ε (see Refs. [16,18] for more details).One starts by considering an oscillatory dynamical system dy dt = f (y) possessing a stable limit cycle Y(t) = Y(t + T ) in its state space.As is well-known, on the limit cycle and in its basin of attraction one can define the phase ϕ = Φ(y) which grows uniformly in time with basic frequency ω = 2π/T .Notice, that on the limit cycle the system's state is uniquely defined by the phase: y = Y(ϕ).Outside of the limit cycle, this is not true: one also has to know the deviation from the cycle.
In the context of oscillatory networks, one considers many coupled oscillators that we label by index k: Here ε is the small parameter governing the strength of the coupling.The equation for the phases is obtained by exploiting their definition ϕ k = Φ k (y k ): This equation is exact, but it contains the full state space trajectories {y j }.However, for a small perturbation these trajectories are close to the limit cycle, up to deviations of order ∼ ε.Thus, substituting the zero-order approximation y j ≈ Y j into the r.h.s. of (2), we obtain equations, where only the states on the limit cycles appear, and these states are unambiguously determined by the phases: As is clear from the presented consideration, to extend the reduction beyond the first-order approximation, one has to calculate the deviations from the limit cycle of orders ∼ ε, ε 2 , . . .and to express these deviations in terms of the phases.Below, we do not provide a general approach valid for arbitrary systems, but restrict ourselves to the simplest case where the phase can be introduced analytically.

A network of Stuart-Landau oscillators
It is instructive to introduce the Stuart-Landau oscillators in dimensional variables, to understand the meaning of dimensionless parameters in the problem.We write the equation for the complex amplitude of oscillations Z for a particular oscillator as where G is the coupling term depending on the states of other elements of the network.Notice that for brevity and without loss of generality we omit the index for the considered oscillator and label other units by the index k = 1, 2, . ... In Eq. ( 3), all the parameters are dimensional, and it is convenient to perform a transformation to dimensionless variables and parameters.Here one should also take into account, that generally the parameters might be different for different oscillators.Using only local parameters, one can introduce the dimensionless amplitude Z = µ/βA, so that all uncoupled oscillators will have amplitude one.Another variable that we can scale is the time, and here one has to choose some common scaling for all oscillators.It appears convenient to assume that the growth rate of linear oscillations µ is the same for all oscillators (the results can be straightforwardly generalized also to the case of different µ), and use it to introduce dimensionless time as t = µτ .Then we obtain where Ȧ now means derivative with respect to t.Here ω = ν/µ−γ/β is the dimensionless frequency of the limit cycle oscillation and the limit cycle has unit amplitude |A (0) | = 1.
Parameter α = γ/β measures non-isochronicity, as we will see below, it determines the phase definition function Φ(A).The dimensionless coupling parameter ε depends on the scaling of the function G.If G contains first powers of the amplitudes Z, then ε = /µ, where we assume that µ/β ≈ µ k /β k , i.e. all the coupled oscillators Z have similar amplitudes.As we see, all the dimensionless parameters entering the problem, namely ω, α, ε, are the original parameters of the system normalized by the linear growth rate of oscillations µ.
For further analysis it is instructive to write the Stuart-Landau equations in polar coordinates, with the amplitude R and the angle θ, so that A = Re iθ and It is straightforward to check that, in the absence of coupling, the quantity ϕ = θ−α ln R grows uniformly in time and, hence, is the true phase.Therefore, we rewrite Eqs.(5,6) in terms of R, ϕ: Here we have explicitly used that the coupling terms depend on the amplitudes and the phases of other oscillators.

Outline of the perturbation method
Now, we outline the exploited perturbation procedure, while in the next subsection we will elaborate on it for a particular example.The idea is to look for a dynamics, in which the amplitudes are enslaved by the phases.Namely, we assume that the amplitude R is a function of the phases: and that the dynamics of the phases is also represented as a power series in ε: We substitute these expressions, together with the similar expressions for R k , ϕ k , into Eqs.(7,8).Equating the terms in each power of ε, we obtain the unknown functions r (1) , r (2) , Ψ (1) , Ψ (2) , . ... We illustrate here the first few steps.In the equation for the phase, the first-order approximation, as discussed above, corresponds to taking into account only the leading term in Eq. ( 9), this yields The equation for the amplitude in the first order is obtained by substituting Eq. ( 9) in Eq. ( 7): We express the time derivative of r (1) via partial derivatives with respect to the phases, and insert the zero-order expressions for these derivatives: dr (1)  dt = ∂r (1)  ∂ϕ φ + ∂r (1)  ∂ϕ 1 φ1 + ∂r (1)  ∂ϕ 2 φ2 + . . .≈ ≈ ∂r (1)  ∂ϕ ω + ∂r (1)  ∂ϕ 1 ω 1 + ∂r (1)  ∂ϕ 2 ω 2 + . . . .
Thus, the problem of determining first-order correction to the amplitude reduces to the partial differential equation 2r (1) + ∂r (1)  ∂ϕ ω + ∂r (1)  ∂ϕ Here, we used that the function G is 2π-periodic with respect to the phases and, hence, the r.h.s.can be written as a multiple Fourier series.Similarly, we expand the function r (1) as: which finally yields an expression for the Fourier coefficients of r (1) These are the basic steps in the perturbation expansion.The expressions for r (1) , r k , being substituted in Eq. ( 8) provide phase equations in order ∼ ε 2 .Equations for r (2)  are partial differential equations of type (12) with r.h.s containing also Ψ (1) , r (1) , etc.

Example: Three coupled Stuart-Landau oscillators
In this Section, we exemplify the perturbative procedure for the derivation of phase dynamics equations in a higher order of the perturbation parameter ε with a particular configuration of three coupled Stuart-Landau oscillators.

Configuration of the network.
We consider an array of three elements with the coupling structure 1 ↔ 2 ↔ 3 and write the equations in the form (4): Notice that the oscillators have different frequencies.Introducing the amplitudes and the angle variables according to A k = R k exp[iθ k ] and the phases ϕ k = θ k − α ln R k , we obtain a system of equations in the form (7-8): 2.4.2.Small parameter expansion.Now we expand the amplitude deviations in powers of ε (below k = 1, 2, 3): According to this, the angles θ k can be represented as Also the ratios of the amplitudes, entering (15), can be expressed as Substituting these expansions in Eq. ( 15), we obtain the following expressions for the dynamics of the phases, up to the order ε 3 : 2 ) 2 + 0.5(r Next, we have to evaluate corrections r k , . ... This is accomplished by substituting expressions (16) in the equations for the amplitudes in (15).Here the time derivatives are calculated according to the chain rule, as the corrections r k , . . .are assumed to be functions of the phases ϕ k .For the sake of brevity, we present only the formulas for r 1 : (2) The partial differential equations for r k , . . .are straightforwardly solved in the Fourier representation, because the variations of the amplitudes are 2π-periodic functions of the phases, as outlined in discussion of Eq. ( 13) above.As a result, we obtain in the 1st order in ε: Substitution of these expression in Eq. ( 19) completes the second-order phase reduction and yields closed equations: The coefficients of the second-order coupling terms, denoted in Eqs.(22)(23)(24) by a k;l , are listed in Tables A1, A2, A3).Here the 3-component vector l = (l 1 , l 2 , l 3 ) is used to signify the term with the combination of the phases l 1 ϕ 1 + l 2 ϕ 2 + l 3 ϕ 3 .Furthermore, in Table B1 we list the terms (without coupling coefficients) appearing in orders ε 3 , ε 4 .
We now shortly discuss the physical meaning of the terms that appear in higher orders in ε.
(i) In the second-order approximation, there are no correction terms to the first-order couplings.These corrections appear in the third order (and, presumably, in all odd orders).
(ii) There are terms, which can be roughly described as "squares" of the basic coupling terms; for the dynamics of the 1st oscillator ϕ 1 these are constant terms and terms containing the second harmonics of the phase difference, e.g.∼ sin(2ϕ 2 − 2ϕ 1 ).These high-order terms do not arise from the interaction within the whole network; they appear already in a system of two coupled oscillators.
(iv) Terms containing phase differences for not directly coupled oscillators, (e.g., the term ∼ sin(ϕ 3 − ϕ 1 ) on the r.h.s. of the equation for φ1 ) mean that connections in terms of phase dynamics do not coincide with the "structural" connections in the original formulation.
(v) While the first-order coupling terms are frequency-independent, the second-order terms depend explicitly on frequency differences.In a more general setup (cf. a system of coupled van der Pol equations treated below) we expect that coupling will depend on the frequencies themselves.

Numerical phase reduction
Stuart-Landau oscillator represents an exceptional case when the phase and the instantaneous frequency can be directly obtained from equations.We exploited this feature to derive the second-order phase dynamics equation in the previous Section.For a general oscillator, one has to evaluate the phases numerically.This immediately provides the first-order approximation of the phase dynamics via numerically calculated phase sensitivity functions.
In this Section, we describe a numerical procedure for determining coupling functions in higher orders, by virtue of the phase analysis of numerically obtained trajectories of the full system.We will first verify it by comparing the results for the Stuart-Landau model (14) with theory in Section 2, and then apply it to a network of three interacting van der Pol oscillators.

Numerical computation of phases and their derivatives
The first step in numerical analysis is the determination of phases ϕ i and their derivatives φi for all elements of the analyzed network, i = 1, . . ., N .For this purpose, we extend the technique suggested in [43], where it was described how to obtain phases from an arbitrary trajectory.
Consider a particular oscillator within a network described by Eq. ( 1).Omitting the index of this unit for simplicity of presentation, we write where the last term describes the coupling to all other units.As a preparatory step we compute the autonomous period T .It is, for ε = 0 we take an arbitrary point on the limit cycle Y, and assign to it ϕ = 0; the return time to this point is exactly T .Now, for any other point on the limit cycle we can compute the time τ (Y) required to reach the zero point where ϕ = 0 or, equivalently, ϕ = 2π.Since the true phase grows linearly in time and gains 2π with one revolution around the cycle, its value can be obtained as T .The next step is to compute ϕ(t) and φ(t) for an arbitrary trajectory, on which at some time t the oscillator has the state u = y(t), and the time derivative of this state is v = ẏ(t).To this end, we introduce an autonomous copy of the investigated unit: and let this auxiliary system evolve from initial conditions w(0) = u, for the time interval nT , where the integer n shall be large enough to ensure that the trajectory attracts to the limit cycle.(Practically, we stop the evolution when w ((n − 1)T ) − w(nT ) is smaller than a given tolerance.)Since the time interval of the evolution is a multiple of the period, the initial point w(0) = u and the end point w(nT ) = w have same value of the phase.The point w is on the limit cycle, and, hence, its phase ϕ can be easily computed as described above, and ϕ(u T is exactly the desired phase of the oscillator at the state u. In order to obtain the phase derivative we also have to follow the autonomous evolution (26) towards the limit cycle of the initial condition u + vdt.(Practically, it can be performed simultaneously with the evolution of the point u.) Since dt is (infinitely) small, this can be done by tracing the linear evolution of vdt to vdt within time interval nT .The law of this linear evolution is given by the Jacobian of the original equations (26).Thus, two states u and u + vdt of the coupled system (25) map to two points w and w + vdt on the limit cycle of the autonomous system (26).These points are characterized by the phases ϕ and ϕ + dϕ, respectively.On the other hand, let us consider the evolution of the autonomous system, from the point w to w + vdt.This evolution occurs along the limit cycle of the system (26) within time interval dt.(Notice that generally dt = dt).The evolution is governed by the flow on the cycle, i.e. w + vdt = w + f ( w)dt, what yields dt = dt (v • f ( w)) / f ( w) 2 .Accordingly, phase growth is determined by the natural frequency ω = 2π/T , i.e. dϕ = ωdt, which finally yields Thus, with the presented algorithm we can compute phases and their derivatives as time series with an arbitrary time step and of sufficient length.Certainly, this can be done for all elements of the network.

Numerical reconstruction of the phase dynamics equations
Now, we discuss how the phase dynamics equations of a network can be constructed numerically from given phases and their derivatives.To explain the procedure, it is convenient to denote the r.h.s. of Eq. ( 10) as Q k (ϕ 1 , ϕ 2 , . . ., ϕ N ), where the subscript k = 1, 2, . . ., N labels the oscillator.Q k are commonly referred to as the coupling functions.Since each Q k is 2π-periodic with respect to all its arguments, we can write it as a multiple Fourier series where l = (l 1 , l 2 , . . ., l N ) and ϕ = (ϕ 1 , ϕ 2 , . . ., ϕ N ) are N -dimensional vectors of integer indices and phases, respectively.ϕ • l = N j=1 l j ϕ j denotes the scalar product between the vectors of phases and the mode indices.Notice the relation between the Fourier coefficients a k;l , b k;l , and coefficients in Eqs.(22)(23)(24): k;l . . . .Thus, we have time series of the phases ϕ 1,n , ϕ 2,n , . . ., ϕ N,n and their derivatives φ1,n , φ2,n , . . ., φN,n , where n = 1, 2, . . ., L is the point index, Eq. ( 27) becomes a system of L linear equations for the unknown coupling coefficients a, b.It is natural to approximate Q k by a finite Fourier series, preserving only M L terms in the sum in Eq. ( 27).Practically, we vary the mode indices in the range |l j | ≤ m, which leaves M = 2(m + 1)(2m(2m + 1) + 1) − 1 unknown coefficients.We apply the least square fit method to find the Fourier modes in the coupling from the over-determined linear system.Practically this is accomplished via the Singular Value Decomposition as described in [44].
We denote the Fourier coefficients of the truncated series obtained from this procedure by capitals A k;l , B k;l .Thus, as a result of the numerical evaluation we obtain an approximation of the phase dynamics as The sum in Eq. ( 27) goes over the combination of mode indices such that either l or −l is counted.The success of this numerical approach depends on how how interdependent the series φ k are.Here we use two different protocols.
Asynchronous case.Suppose that the network does not synchronize.It means that the trajectories of the system (1) are dense on the N -dimensional torus.Then one just has to calculate one long trajectory of the system (1) starting from some initial conditions and compute phases and instantaneous frequencies for each point of the solution, as described above.If the data series is sufficiently long, the computed trajectory covers the surface of the torus spanned by ϕ, ϕ 1 , ϕ 2 , . . ., ϕ N .Hence, we obtain enough information to recover the function Q of N variables and to determine the Fourier coefficients.
Synchronous case.If the network synchronizes, the trajectory on the torus collapses to a closed line and Eqs.(27) for the Fourier coefficients become dependent and cannot be solved.However, there is a way to obtain enough data to solve the problem even in this case.For this goal one considers not one long trajectory, but an ensemble of short transients (cf.[45,46]).Namely, one starts numerical integration with some asynchronous initial conditions and follows the trajectory unless it is on the invariant torus.The corresponding transient time should be larger than the amplitude relaxation time, but smaller than the synchronization time of the phases.Then the procedure is repeated with new initial conditions.Thus, instead of one long record one collects many dynamical states on the torus, unless the sufficient number of points is obtained.Certainly, this protocol can be use for the asynchronous case as well.
We exemplify these two protocols in the next Section, but before proceeding with the examples, we have to discuss an important issue.As we mentioned, the least squares optimization requires that data points fill the surface of an N -dimensional torus.Thus, we face the curse of dimensionality: the data requirement grows fast with the network size N and becomes hardly feasible already for N > 3. Some information about the network, e.g., strength of directed links, can, however, be revealed for N > 3 as well.A possible approach is to perform the triplet-based analysis [40,[47][48][49].
Results for the asynchronous configuration of three coupled Stuart-Landau oscillators, see Eq. (15).Differences between the numerically calculated coupling coefficients and the theoretical values given by Eqs.(22)(23)(24) are shown as a function of the coupling strength, for all oscillators.Panel (a) presents this difference for the terms appearing in the 1st order in ε; the black dashed line here corresponds to ∼ ε 3 .Panel (b) presents the terms appearing in the 2nd order in ε; the magenta dashed line here corresponds to ∼ ε 4 .Red squares (green circles) represent cosine (sine) coefficients.Additional black markers in (b) show the differences between the zero-order terms A k;0 , k = 1, 2, 3 and their theoretical values in the second-order approximation.
our theoretical conclusion that there are no second-order correction to the first-order terms.Thus, numerical results exhibit a good correspondence with the theory.Figures 3,4 present an overview of all Fourier coefficients for which we do not have theoretical values.Namely, we show here all coefficients except for those entering Eqs.(22)(23)(24) and analyzed in Figs.1,2.The results are in full agreement with the powerseries representation of the coupling terms.Indeed, we see four groups of coefficients that scale as ε 3 , ε 4 , ε 5 and ε 6 , respectively.
Finally, we mention that the overall precision of the numerical procedure can be estimated by computing the rest term ξ k of the Fourier series representation in Eq. ( 28).The results presented in Fig. 5 show that the rest term grows with the coupling strength as ε 9 .This is an indication that all the terms of orders from 0 to 8 are within the set of chosen Fourier modes.

A network of van der Pol oscillators
Up to this point, our analysis was restricted to the case of Stuart-Landau oscillators where expressions for phases and their derivatives are known.In this Section we present  a purely numerical evaluation of high-order coupling terms for a network of three nonidentical van der Pol oscillators: ẍ1 We fixed parameters ω 1 = 1, ω 2 = 1.324715957, ω 3 = ω 2 2 (here ω 2 is the spiral mean, the root of cubic equation ω 3  2 − ω 2 − 1 = 0).We fixed µ = 1 and varied the coupling constant ε in the range 0.001 ≤ ε ≤ 0.3, in this range the dynamics of the network is asynchronous.From a trajectory of the system (30) we obtained time series ϕ 1,2,3 , φ1,2,3 of length L = 10 6 points each.We used this set to estimate the coefficients of the truncated Fourier series in Eq. ( 28) for m = 4. Thus, 729 coupling constants A k;l , B k;l were calculated for each ε.Below we restrict our attention only to the strengths of the coupling and ignore the relative phase, so we look on the properties of 364 coupling constants H k,l = (A 2 k;l + B 2 k;l ) 1/2 .Together with the free term A k;0 this constitutes a set of 365 numerically obtained coefficients for each oscillator and for each coupling constant ε.
For presentation of the results we perform a preliminary sorting.According to the theory, we expect to obtain terms with leading dependencies ∼ ε q , with q = 0, 1, 2, . ... Therefore, for each set of indices l and each oscillator, we tried to approximate the coefficients H k,l (ε) by the function ∼ ε q , and if the fitted value q was close to an integer and the reliability of the fit was high, we attributed the corresponding power to the coupling term.Additionally, for ε 0 , ε, ε 2 we used the five small values of ε = 0.001, 0.002, 0.005, 0.01, 0.02, while for powers ε 3 and ε 4 we used larger coupling constants ε = 0.04, 0.06, 0.08, 0.1, 0.15.We did not look for powers 5 and larger.
Figure 6 illustrates these findings.It is instructive to look which coupling modes appear in each power of ε.For modes up to ∼ ε 2 we summarize this in Table 1.The modes A k;0 , describing corrections to the natural frequency are not listed.Notice that in each cell of the table the modes are ordered according to their amplitudes.Below we summarize the properties of the coupling modes of the van der Pol oscillators: (i) Looking at the modes that appear in the first order in ε, we notice that the terms with phase differences and the terms with phase sums have nearly the same amplitude.This means that the coupling terms ∼ ε have nearly the Winfree form: they are products of two functions of the oscillator phase and of the driving oscillator phase, in full agreement with the first-order theory.Recall that it is different for the Stuart-Landau oscillators, where only terms with phase differences appear due to the condition (29).
(ii) Because the variable x of the van der Pol equation possesses odd harmonics of the phase (and the same is true for the phase response curve), terms with the third harmonics appear already in the 1st order in ε.
(iii) Inspection of terms ∼ ε in Table 1 reveals some terms that are not expected in the first-order analysis, we show them in italic font.Data in Fig. 6 show that the amplitudes of these terms are extremely small, comparable to the errors in terms reconstruction.We conclude that these terms are probably spurious and just occasionally possess a scaling ∼ ε, what is not surprising due to the fact that the total number of terms to be found is large.
(iv) With the size of the time series explored, we could not reliably detect coupling terms appearing in order ε 5 .However, as the figures show, the terms ∼ 4 can be detected with good confidence.

Conclusion
In this work, we have presented an analytic perturbation approach, allowing the derivation of the equations for the phase dynamics for general networks of Stuart-Landau oscillators.We exemplified this framework with calculations for a particular system of three units.We demonstrated explicitly that already in the second order in coupling strength, there exist coupling terms that are not present in the structural coupling configuration.This result confirms a general statement that phase connectivity generally differs from the structural connectivity of a network.In particular, in higher-order approximation we find triplet coupling terms, which are characteristic for a hypernetwork.
Since analytic derivations of the phase dynamics equations for generic oscillators remain a theoretical challenge, we developed a numerical method to compute the phase dynamics of an oscillatory network and to extract the coupling terms from time series of the obtained phases.These data could be one long trajectory if the dynamics is quasiperiodic, or multiple small pieces (or points) if the dynamics is synchronous.The result of the numerical procedure is a set of Fourier modes of the coupling functions.Analysing these sets for different coupling strengths, we found terms with various powerlaw dependencies on this strength.We first tested this approach on the coupled Stuart-Landau oscillators, where we demonstrated an excellent agreement with the theory.This numerical approach has also been applied to a network of van der Pol oscillators, coupling functions of which are much more involved.

Figure 4 .
Figure4.The same as in Fig.3, but for the synchronous configuration.

Figure 5 .
Figure 5. Accuracy of the phase dynamics reconstruction for all oscillators.Red squares and blue circles represent the rest term ξ k of the Fourier representation Eq. (28) for asynchronous and synchronous cases, respectively.The dashed line shows power law ∼ ε 9 . a

Figure C5 .
Figure C5.Accuracy of the phase dynamics reconstruction.Red squares and blue circles illustrate asynchronous and synchronous cases, respectively.The dashed line shows power law ∼ ε 5 .