Inverse problems for dynamic patterns in coupled oscillator networks: when larger networks are simpler

TL;DR

This paper introduces a method to infer model parameters in large networks of oscillators by analyzing their synchronized dynamics.

Contribution

The paper derives statistical equilibrium relations for large oscillator networks and uses them for parameter reconstruction.

Findings

Statistical equilibrium relations can be derived for partially synchronized patterns in large oscillator networks.

The method is effective for large networks, noisy data, and partial observations.

The approach is demonstrated on chimera states and extended to other network topologies.

Abstract

Networks of coupled phase oscillators are one of the most studied dynamical systems with numerous applications in physics, chemistry, biology, and engineering. Their behaviour is often characterized by the emergence of various partially synchronized dynamic patterns, which in the case of large networks can be analysed using a variant of the mean-field approach. This method allows to predict what type of network dynamics can be observed for different system parameters. But it is less known that for different partially synchronized patterns it also allows to obtain statistical equilibrium relations that express the dependence of some time-averaged observable quantities of individual oscillators on the internal parameters of these oscillators and the interaction functions between them. In this paper, we show how such relations can be derived, what their typical accuracy is for finite-size networks, and how they can be used to reconstruct the parameters of the corresponding model. The proposed method is particularly effective for large networks, for unevenly sampled or noisy observables, and for partial observations. Its possibilities are demonstrated by application to chimera states in networks of phase oscillator with nonlocal coupling. The extension of the method to other systems with all-to-all and random network topologies is also described.

Inferring model parameters from partially synchronized dynamics is a major challenge in networks of coupled oscillators. Here, the author derives statistical equilibrium relations valid for large networks and uses them to reconstruct system parameters from time-averaged observables.

Figures

Application of the parameter reconstruction algorithm to partial data.For a chimera state in model (5) with top-hat coupling and *N* = 8192 oscillators, we used 1996 randomly selected oscillators to reconstruct the parameters *ω*, *β* and *G*(*x*). The quantities shown in (a–c) are the same as in Fig. [4](#Fig4). Other model parameters are given in Fig. [2](#Fig2).

Schematic representation of the proposed parameter reconstruction method.Given a complex spatio-temporal pattern in a system of coupled phase oscillators, the model parameters can be reconstructed by calculating a small number of averages and using statistical equilibrium relations relevant to this model.

Time-averaged characteristics of chimera states.A typical chimera state in system (5) for a top-hat coupling function with *σ* = 0.7, *ω* = 1, *α* = *π*/2 − 0.1, and *N* = 1024. a Space-time plot of *θ*_*j*_(*t*). b, c Effective frequencies Ω_*j*_ and local order parameters *ζ*_*j*_ obtained by averaging over 2000 time units. Every 16th point *x*_*j*_ is shown.

Accuracy of statistical equilibrium relations for finite *N.*Mean discrepancies of SERs (7)–(9) for the chimera state from Fig. [2](#Fig2). The discrepancies *δ*_1_, *δ*_2_ and *δ*_3_ defined in the text are shown separately in (a–c) respectively. The five curves in each panel show the results for different system sizes *N* mentioned in panel (a).

Reconstruction of model parameters using statistical equilibrium relations (SER).a SER (7) for the chimera state in Fig. [2](#Fig2). The circles show the averages calculated from the numerical trajectory (only every 16th point is shown), the line shows a linear fit. b The graph of the function *J*_incoh_(*β*). The dashed line shows the position of the minimum $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta

Efficiency of the proposed method in reconstructing the phase lag *β.*The method was applied to the chimera state observed in model (5) with top-hat coupling. The dots show the values of $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{\min }$$\end{document}$ found as the global minimum of the function *J*_incoh_(*β*) for different *β*. Parameters: *N* = 2048, *ω* = 1 and *σ* = 0.7.

Efficiency of the proposed method in reconstructing the coupling function *G*(*x*).The six leading Fourier coefficients *c*_*m*_ in formula (14) were recovered from the chimera state observed in model (5) with top-hat coupling. The curves show the theoretical values given by formula (16) and the symbols show the reconstructed values. Parameters: *N* = 2048, *ω* = 1 and *β* = 0.1.

Other examples of application of the parameter reconstruction algorithm to chimera states in model (5) with *N* = 2048.a–c Exponential coupling function with *κ* = 0.5, *ω* = 1, and *α* = *π*/2 − 0.1. d–f Cosine coupling function with *A* = 0.9, *ω* = 1, and *α* = *π*/2 − 0.15. In both cases, the coupling function *G*(*x*) was approximated by the ansatz (14) containing *M* = 10 spatial Fourier modes. In (c, f), the red/dark curves show reconstructed coupling functions, while the grey/light curves show the original coupling functions. Other notations are the same as in Fig. [4](#Fig4).Fig. 8Rec

Application of the parameter reconstruction algorithm to partial data.For a chimera state in model (5) with top-hat coupling and *N* = 8192 oscillators, we used 1996 randomly selected oscillators to reconstruct the parameters *ω*, *β* and *G*(*x*). The quantities shown in (a–c) are the same as in Fig. [4](#Fig4). Other model parameters are given in Fig. [2](#Fig2).

Full text

Introduction

One of the main goals of natural sciences is to predict the behaviour of a given system, assuming that changes in its state are determined by certain dynamical rules expressed by differential equations. In some cases, these equations can be derived from first principles and the results of specially designed experiments, but more often they have to be obtained from uncontrolled observational data. This duality is reflected in the coexistence of two general approaches to the identification of dynamical systems: model-based and data-driven. At first glance, the latter approach seems to be more versatile, as it relies on a minimal amount of information about the system, such as the assumption of sparsity of the governing equations^1^. However, its implementation typically requires a large amount of data (e.g., many trajectories passing through different parts of the phase space) and can become computationally cumbersome as the system size increases. To overcome these difficulties, a number of more sophisticated methods have been proposed, including an equation-free method for inferring coarse-grained multiscale dynamics^2^, automated adaptive model inference^3^, data-driven discovery of intrinsic lower-dimensional dynamics^4,5^, machine learning techniques based on reduced order models^6^, and others^7^. A common feature of all these methods is that they attempt to approximate the behaviour of a complex large-scale system using a phenomenological lower-dimensional model, although they utilize this simplification ansatz almost heuristically.

A similar dimensionality reduction scheme also exists in the model-based approach. But there it is better justified and can be used more effectively and purposefully. Roughly speaking, it is well-known that large systems of many interacting agents have the property of coordinating their behaviour in such a way that it is described by the laws of statistical physics. This means that no matter how complex the microscopic dynamics of the system is, it is characterized by a certain statistical balance between the dynamics of the constituent agents and their intrinsic properties. Usually, this relationship is described at the macroscopic level using global coarse-grained variables and some form of mean-field analysis, while the detailed balance at the microscopic level remains in the shadows. In this paper, we show that mathematical formulas expressing this detailed balance can actually be very useful, in particular, for reconstructing the parameters of the corresponding high-dimensional dynamical systems. The general scheme of the proposed approach is described in the context of its application to complex dynamic patterns in networks of coupled phase oscillators. Using it, we formulate a parameter reconstruction algorithm that is non-invasive, fast, easy to compute, suitable for partial observation and robust to measurement noise.

Mathematical models describing the collective behaviour of large populations of coupled phase oscillators can be found in various fields of physics, chemistry, and biology^8,9^. They play a key role in the study of synchronization phenomena^10–12^ and have a direct connection to more complex real-world models through the standard phase reduction procedure^10,13–15^. Even without a rigorous justification from first principles, such models are often used in theoretical biology and neuroscience to explain observed properties of dynamical quorum sensing^16,17^, circadian rhythm generators^18,19^, metachronal waves in cilia carpets^20^, brain disorders^21,22^, and other physiological processes related to synchrony^23^. In general, these models are defined as follows. The population consists of N oscillators, the state of each of which is described by a scalar quantity, its phase θj. Each oscillator has a label \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${p}_{j}\in {{\mathbb{R}}}^{{N}_{{{{\rm{p}}}}}}$$\end{document} containing information about its intrinsic properties (e.g., natural frequency, position in space, etc.), which remain unchanged over time. Accordingly, the dynamics of this oscillator is determined by a differential equation

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{d{\theta }_{j}}{dt}=F({\theta }_{j},{p}_{j},{W}_{j}),$$\end{document}

where

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${W}_{j}=\frac{1}{N}{\sum }_{k=1}^{N}Q({p}_{j},{p}_{k},{\theta }_{k})$$\end{document}

is a mean-field acting on the jth oscillator due to the influence of all other oscillators. Note that, despite their simple structure, Eqs. (1), (2) describe a broad class of coupled oscillator networks, including fully connected and spatially extended networks, as well as annealed approximations of random networks.

In the thermodynamic limit, when the number of oscillators N tends to infinity and the distribution of labels pj converges to some probability density g(p), it is often observed that after a sufficiently long transient, the state of system (1) approaches some statistical equilibrium. In the mean-field approximation, this equilibrium is characterized by a single particle probability density function ρ(θ, p, t). Using this function, we can replace sum (2) with the integral

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${W}_{j}\,\mapsto \,{{{\mathcal{W}}}}[\rho ]({p}_{j})={\int }_{{\!\!\!\!{\mathbb{R}}}^{{N}_{{{{\rm{p}}}}}}}{\int }_{\!\!\!\!-\pi }^{\pi }Q({p}_{j},p,\theta )\rho (\theta,p,t)d\theta \,dp$$\end{document}

and write a nonlinear integro-differential continuity equation

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{\partial \rho }{\partial t}+\frac{\partial }{\partial \theta }\left(F(\theta,p,{{{\mathcal{W}}}}[\rho ])\rho \right)=0$$\end{document}

which describes the evolution of ρ(θ, p, t).

Although Eq. (3) looks more complicated than the original oscillator system (1), its solution representing the statistical equilibrium of (1) has usually a much simpler form than the corresponding oscillators’ trajectory. In many cases, this solution ρse(θ, p, t) can be written in analytical (but not necessarily explicit) form, using some kind of self-consistency analysis^10^. Then due to the ergodicity property of statistical equilibrium the solution ρse(θ, p, t) can be used to derive statistical equilibrium relations, i.e. formulas relating the time-averaged observables in system (1) and the parameters of this system. For example, one of the most common quantities characterizing the dynamics of the jth oscillator is its effective frequency, which is defined as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\Omega }_{j}=\left\langle \frac{d{\theta }_{j}}{dt}\right\rangle,$$\end{document}

where 〈 ⋅ 〉 denotes time average. Using Eq. (1), the same value can also be written as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\Omega }_{j}=\left\langle {\int }_{\!\!\!\!-\pi }^{\pi }\frac{{\rho }_{{{{\rm{se}}}}}(\theta,{p}_{j},t)}{g({p}_{j})}F(\theta,{p}_{j},{{{\mathcal{W}}}}[\rho ]({p}_{j}))d\theta \right\rangle,$$\end{document}

where the use of the conditional probability density ρse(θ, pj, t)/g(pj), is due to the fact that we are considering an oscillator with label pj. Formula (4) gives an algebraic relationship between the time-averaged observable Ωj and the system parameters {pj}. In other words, it expresses the microscopic balance between the dynamics of individual oscillators and their intrinsic properties, and can therefore be considered as a statistical equilibrium relation.

Similar relations, but for other time-averaged quantities, will be described below. In addition, we will show how they can be used to reconstruct the parameters of model (1), (2). For clarity, we will focus on a special but important case — the Kuramoto-Battogtokh system of nonlocally coupled phase oscillators^24^. It is famous as a prototype system for chimera states^25–27^, which are dynamic patterns with self-organized domains of synchronized (coherent) and desynchronized (incoherent) behaviour.

Results

The structure of this section is graphically presented in Fig. 1. First, we describe the Kuramoto-Battogtokh system and show a typical example of a chimera state. Then, we define additional time-averaged quantities, the local order parameters, and write down statistical equilibrium relations for them. (The mathematical details of their derivation can be found in the section Methods.) Finally, we describe our parameter reconstruction algorithm and demonstrate its effectiveness on various examples.Fig. 1. Schematic representation of the proposed parameter reconstruction method.Given a complex spatio-temporal pattern in a system of coupled phase oscillators, the model parameters can be reconstructed by calculating a small number of averages and using statistical equilibrium relations relevant to this model.

Model

We consider a ring of N nonlocally coupled identical phase oscillators

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{d{\theta }_{j}}{dt}=\omega -\frac{2\pi }{N}{\sum }_{k=1}^{N}G({x}_{j}-{x}_{k})\sin ({\theta }_{j}-{\theta }_{k}+\alpha ).$$\end{document}

Here, ω is the natural frequency of all oscillators and α is the Kuramoto-Sakaguchi phase lag parameter. The position of the jth oscillator is given by xj ∈ [ − π, π] and the nonlocal interaction between oscillators is determined by a scalar symmetric 2π-periodic coupling function G(x). More precisely, the positions xj are assumed to be uniformly distributed on the interval [ − π, π], although in most numerical examples below we will use a special deterministic choice xj = − π + 2π**j/N.

Pattern

It is well-known^24,25^ that for a wide range of parameters in (5) this system exhibits peculiar spatio-temporal patterns, where some oscillators rotate almost synchronously, while others exhibit mutually asynchronous behaviour, see Fig. 2a. In the literature, they are usually called coherence-incoherence patterns or chimera states. Nonlocal couplings for which chimera states have been found include exponential function^24^

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G(x)=\frac{\kappa }{2(1-{e}^{-\pi \kappa })}{e}^{-\kappa \arccos (\cos x)},\,\kappa > 0,$$\end{document}

cosine function^25^

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G(x)=\frac{1}{2\pi }(1+A\cos x),\,A > 0,$$\end{document}

top-hat function^28^

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G(x)=\frac{1}{4\pi \sigma }\left(1+\frac{\pi \sigma -\arccos (\cos x)}{| \pi \sigma -\arccos (\cos x)| }\right),\,0 < \sigma < 1,$$\end{document}

and many others^27^. (Note that due to the periodic boundary conditions in model (5), above we used the expression \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\arccos (\cos x)$$\end{document} , which is equal to ∣x∣ if ∣x∣≤π, and defines a 2π-periodic extension of ∣x∣ if ∣x∣ > π.)Fig. 2. Time-averaged characteristics of chimera states.A typical chimera state in system (5) for a top-hat coupling function with σ = 0.7, ω = 1, α = π/2 − 0.1, and N = 1024. a Space-time plot of θj(t). b, c Effective frequencies Ω_j_ and local order parameters ζj obtained by averaging over 2000 time units. Every 16th point xj is shown.

The initial interest in chimera states was purely theoretical. But later their existence was confirmed experimentally in systems of chemical Belousov-Zhabotinsky oscillators^29,30^ and in systems of electrochemical oscillators^31^. In addition, their similarity to dynamic patterns in various biological systems has been established. These include synchronization patterns of elastic cilia^32,33^, collective states of coupled inner-ear hair cells^34^, and epileptic seizures^35^. Although the functional role of chimera states remains unclear, one could consider using them to obtain information about the chemical or biological system in which they occur. For example, the function G(x) in Eq. (5) contains important characteristics of the nonlocal coupling, such as its range, its monotonic (or not) dependence on distance, and its decay rate. The phase lag α is a measure of the nonreciprocity of the interaction between oscillators^36^, while the natural frequency ω is related to the properties of the oscillators in isolation. So, what can we do to find all these interesting parameters in a situation where they cannot be measured directly? More specifically, we can ask the following questions. (i) Can these parameters be determined from the observation of a single chimera state in system (5)? (ii) And if so, how can this be done effectively? Below we will give an affirmative answer to the first question and propose a relatively simple algorithm for solving the second question.

At first glance, the following approach seems to be the most natural to address the parameter reconstruction problem for system (5). Insert the observed trajectory {θj(t)} and its derivative \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left\{\frac{d{\theta }_{j}}{dt}(t)\right\}$$\end{document} into Eq. (5) and solve the resulting system with respect to the unknown parameters^37^. However, this method has a number of disadvantages. First, its implementation requires knowledge of the trajectory with high time resolution for accurate calculation of derivatives. Second, the calculations use the entire trajectory {θj(t)} as a rectangular matrix, which becomes extremely huge for large system sizes N. Third, the trajectory of system (5) must be complete, that is, the behaviour of all oscillators must be known.

Averages

Below we describe an alternative parameter reconstruction method that does not have the above drawbacks. It is based on statistical equilibrium relations for system (5) and only requires computing O(N) time-averaged quantities from the trajectory {θj(t)}. Knowledge of derivatives is not required at all. More precisely, for each oscillator θj(t), we only need to calculate its effective frequency Ωj and its local order parameter

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\zeta }_{j}=\left\langle {e}^{i{\theta }_{j}(t)}\frac{\overline{Z}(t)}{| Z(t)| }\right\rangle \in {\mathbb{C}},$$\end{document}

where

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Z(t)=\frac{1}{N}{\sum }_{k=1}^{N}{e}^{i{\theta }_{k}(t)}$$\end{document}

is the global order parameter of all oscillators and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{Z}(t)$$\end{document} is its complex conjugate value, see Fig. 2b, c. Note that after calculating Ωj and ζj, the oscillator trajectory {θj(t)} is no longer needed and does not need to be stored.

Statistical equilibrium relations (SER)

In the thermodynamic limit, chimera states have an analytic representation following from the corresponding continuity equation (3). Using it, we can derive statistical equilibrium relations (see Methods)

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{\omega -{\Omega }_{j}}{\omega -\Omega }=\frac{2| {\zeta }_{j}{| }^{2}}{1+| {\zeta }_{j}{| }^{2}},$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\rm{Re}}}}\left(\frac{{\xi }_{j}}{{\zeta }_{j}}\right)=\frac{2}{1+| {\zeta }_{j}{| }^{2}},$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\xi }_{j}=\frac{2{\zeta }_{j}}{1+| {\zeta }_{j}{| }^{2}}\,\,{{{\rm{for}}}}\,\,| {\zeta }_{j}| < 1,$$\end{document}

where

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\xi }_{j}=\frac{{e}^{i\beta }}{\omega -\Omega }{\sum }_{k=1}^{N}G({x}_{j}-{x}_{k}){\zeta }_{k}\frac{{x}_{k+1}-{x}_{k-1}}{2}$$\end{document}

and

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta=\frac{\pi }{2}-\alpha .$$\end{document}

(Note that in (10) the notations x2 − x0 = 2π + x2 − xN and xN+1 − xN−1 = 2π + x1 − xN−1 are used to represent periodic boundary conditions).

In other words, whatever stationary coherence-incoherence pattern we find in model (5) with N → ∞, relations (7)–(9) will always be satisfied for it, regardless of the natural frequency ω, the phase lag α and the coupling function G(x). Note that by stationarity we mean here that the pattern does not change its position on the network and does not change its shape at the macroscopic level. Therefore, dynamic patterns such as travelling chimera states^38^ and breathing chimera states^39^ are currently excluded from our consideration. Moreover, the requirement of pattern stationarity automatically imposes a reflection-symmetry requirement on G(x), since asymmetric coupling functions usually lead to pattern motion^38^.

In the following, we want to verify whether relations (7)–(9) can be used to recover the main system parameters based on the observation of a stationary chimera state in model (5). More precisely, we assume that the observation is limited to recording a relatively small dataset, including oscillator positions xj, effective frequencies Ωj, and local order parameters ζj.

Practical accuracy of SERs

Thermodynamic limit theory predicts that SERs (7)–(9) are only exact for an infinitely large system size N and for effective frequencies Ωj and local order parameters ζj calculated by infinitely long time averaging. But they also remain approximately accurate under much weaker constraints. For example, let us consider finite-time averages

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\Omega }_{j}(T)=\frac{1}{T}{\int }_{\!\!\!\!0}^{T}\frac{d{\theta }_{j}(t)}{dt}dt,$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\zeta }_{j}(T)=\frac{1}{T}{\int }_{\!\!\!\!0}^{T}{e}^{i{\theta }_{j}(t)}\frac{\overline{Z}(t)}{| Z(t)| }dt,$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega (T)=\frac{1}{T}{\int }_{0}^{T}{{{\rm{Im}}}}\left(\frac{1}{Z(t)}\frac{dZ(t)}{dt}\right)dt,$$\end{document}

and ξj(T) given by formula (10) with ζj(T) and Ω(T). Then, for each of the SERs (7)–(9), we can define its mean discrepancy

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\delta }_{1}(T) \;=\; \frac{1}{N}{\sum }_{j=1}^{N}\left|\frac{\omega -{\Omega }_{j}(T)}{\omega -\Omega (T)}-\frac{2| {\zeta }_{j}(T){| }^{2}}{1+| {\zeta }_{j}(T){| }^{2}}\right|,\\ {\delta }_{2}(T) \;=\; \frac{1}{N}{\sum }_{j=1}^{N}\left|{{{\rm{Re}}}}\left(\frac{{\xi }_{j}(T)}{{\zeta }_{j}(T)}\right)-\frac{2}{1+| {\zeta }_{j}(T){| }^{2}}\right|,\\ {\delta }_{3}(T) \;=\; \frac{1}{{N}_{ * }}{\sum }_{j:| {\zeta }_{j}(T)| < 1-1/\sqrt{N}}\left|{\xi }_{j}(T)-\frac{2{\zeta }_{j}(T)}{1+| {\zeta }_{j}(T){| }^{2}}\right|,$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${N}_{*}$$\end{document} is the number of indices j satisfying the inequality \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$| {\zeta }_{j}(T)| < 1-1/\sqrt{N}$$\end{document} . Calculating these mean discrepancies for the chimera state from Fig. 2, as well as for chimera states with the same parameters but different system sizes N, we see that SERs (7)–(9) are very accurate already for N > 1000 and averaging times T > 1000, Fig. 3. Thus, these relations can also be used in realistic situations where N and T are moderately large. This approach is roughly comparable to the application of the laws of thermodynamics, which are proven by statistical physics for infinitely large systems, but are used for systems consisting of a finite number of particles, provided that this number is large enough.Fig. 3. Accuracy of statistical equilibrium relations for finite *N.*Mean discrepancies of SERs (7)–(9) for the chimera state from Fig. 2. The discrepancies δ1, δ2 and δ3 defined in the text are shown separately in (a–c) respectively. The five curves in each panel show the results for different system sizes N mentioned in panel (a).

Parameter reconstruction algorithm

Suppose that the statistical equilibrium relations (7)–(9) are satisfied (with some discrepancy) for the observables Ωj and ζj. How can we use this fact to reconstruct the parameters ω, β and G(x) in model (5)? It is easy to see that the values of ω and ω − Ω can be found from the statistical equilibrium relations (7) written in the form

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\Omega }_{j}=\omega+(\Omega -\omega ){\eta }_{j}\,\,{{{\rm{with}}}}\,\,{\eta }_{j}=2| {\zeta }_{j}{| }^{2}/(1+| {\zeta }_{j}{| }^{2}).$$\end{document}

Indeed, standard linear regression yields, see Fig. 4a,

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega=\frac{{S}_{\Omega }{S}_{\eta \eta }-{S}_{\eta }{S}_{\Omega \eta }}{{S}_{\eta \eta }-{S}_{\eta }^{2}},\,\Omega -\omega=\frac{{S}_{\Omega \eta }-{S}_{\eta }{S}_{\Omega }}{{S}_{\eta \eta }-{S}_{\eta }^{2}},$$\end{document}

where

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${S}_{\eta }= \frac{1}{N}{\sum }_{j=1}^{N}{\eta }_{j},\,{S}_{\Omega }=\frac{1}{N}{\sum }_{j=1}^{N}{\Omega }_{j},\\ {S}_{\eta \eta } = \frac{1}{N}{\sum }_{j=1}^{N}{\eta }_{j}^{2},\,{S}_{\Omega \eta }=\frac{1}{N}{\sum }_{j=1}^{N}{\Omega }_{j}{\eta }_{j}.$$\end{document}

Fig. 4. Reconstruction of model parameters using statistical equilibrium relations (SER).a SER (7) for the chimera state in Fig. 2. The circles show the averages calculated from the numerical trajectory (only every 16th point is shown), the line shows a linear fit. b The graph of the function Jincoh(β). The dashed line shows the position of the minimum \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{\min }$$\end{document} . c The red/dark curve shows the reconstructed coupling function with M = 10 spatial Fourier modes. The grey/light curve shows the original coupling function G(x) and the dotted curve shows its exact Fourier approximation with 10 modes.

The remaining phase lag parameter β and the coupling function G(x) can be found as follows. First, we note that formula (10) implies

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\rm{Re}}}}\left(\frac{{\xi }_{j}}{{\zeta }_{j}}\right)={\sum }_{k=1}^{N}\frac{G({x}_{j}-{x}_{k})}{\omega -\Omega }{{{\rm{Re}}}}\left({e}^{i\beta }\frac{{\zeta }_{k}}{{\zeta }_{j}}\right)\frac{{x}_{k+1}-{x}_{k-1}}{2}$$\end{document}

for all j = 1, …, N. On the other hand, if G(x) is symmetric, i.e. G( − x) = G(x), then it can be approximated by a Fourier sum

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G(x)={\sum }_{m=0}^{M}{c}_{m}{q}_{m}(x)\,\,{{{\rm{with}}}}\,\,{q}_{m}(x)=\cos (mx).$$\end{document}

In Supplementary Note 1, we explain that for each chimera state there is an optimal number of spatial modes Mopt that corresponds to the best efficiency of the reconstruction algorithm. But it is not known a priori, so M in Eq. (14) is chosen empirically.

According to (8) we can expect that the vector ({cm}, β) is the minimizer of the functional

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J(\{{c}_{m}\},\beta )=\frac{1}{N}{\sum }_{j=1}^{N}{\left[\frac{2}{1+| {\zeta }_{j}{| }^{2}}-{\sum }_{m=0}^{M}{c}_{m}{Q}_{jm}(\beta )\right]}^{2},$$\end{document}

where

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${Q}_{jm}(\beta )={\sum }_{k=1}^{N}\frac{{q}_{m}({x}_{j}-{x}_{k})}{\omega -\Omega }{{{\rm{Re}}}}\left({e}^{i\beta }\frac{{\zeta }_{k}}{{\zeta }_{j}}\right)\frac{{x}_{k+1}-{x}_{k-1}}{2}.$$\end{document}

Note that in order not to lose the information provided by relations (9), we use the minimization problem for J({cm}, β) only to express the coefficients cm as functions of β. For this, we rewrite the corresponding local minimum condition

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\partial }_{{c}_{n}}J(\{{c}_{m}\},\beta )=- \frac{2}{N}{\sum }_{j=1}^{N}{Q}_{jn}(\beta )\left[\frac{2}{1+| {\zeta }_{j} | ^{2}}-{\sum }_{m=0}^{M}{c}_{m}{Q}_{jm}(\beta )\right]=0$$\end{document}

in the matrix form

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A(\beta )c=b(\beta ),$$\end{document}

where

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${A}_{nm}(\beta )={\sum }_{j=1}^{N}{Q}_{jn}(\beta ){Q}_{jm}(\beta )$$\end{document}

and

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${b}_{n}(\beta )={\sum }_{j=1}^{N}\frac{2{Q}_{jn}(\beta )}{1+| {\zeta }_{j}{| }^{2}}.$$\end{document}

Then, the solution of (15) reads

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\widetilde{c}(\beta )={A}^{-1}(\beta )b(\beta ).$$\end{document}

Now, using the statistical equilibrium relations (9) and formulas (10) and (14), we define a function

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${J}_{{{{\rm{incoh}}}}}(\beta )=\frac{1}{N}{\sum }_{j\,:\,| {\zeta }_{j}| < 1-\frac{1}{\sqrt{N}}}{\left|\frac{2{\zeta }_{j}}{1+| {\zeta }_{j}{| }^{2}}-{\sum }_{m=0}^{M}{e}^{i\beta }{\widetilde{c}}_{m}(\beta ){\widetilde{Q}}_{jm}\right|}^{2},$$\end{document}

where

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\widetilde{Q}}_{jm}={\sum }_{k=1}^{N}\frac{{q}_{m}({x}_{j}-{x}_{k})}{\omega -\Omega }{\zeta }_{k}\frac{{x}_{k+1}-{x}_{k-1}}{2},$$\end{document}

and look for its global minimum \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{\min }$$\end{document} , which in theory must coincide with the value of phase lag β in model (5). (Note that in the thermodynamic limit, statistical equilibrium relation (9) holds for all oscillators j with ∣ζj∣ < 1. But because of the finite-size fluctuations we replace this inequality with a more restrictive one \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$| {\zeta }_{j}| < 1-1/\sqrt{N}$$\end{document} in our definition of Jincoh(β).) To find the global minimum of Jincoh(β), we calculate this function at 20 points evenly spaced in the interval [0, π] and use the resulting approximate estimate of the minimizer as an initial guess to solve the equation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J'_{{{{\rm{incoh}}}}}(\beta )=0$$\end{document} using Newton’s method. The obtained global minimizer \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{\min }$$\end{document} is interpreted as an approximate value of the phase lag β in Eq. (5), see Fig. 4b. Respectively, formula (14) with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${c}_{m}={\widetilde{c}}_{m}({\beta }_{\min })$$\end{document} gives an approximate representation of the coupling function G(x), see Fig. 4c.

To complete the description of our algorithm, we will also estimate its computational cost. It is easy to see that the largest matrices used by the algorithm are the N × M matrices Qj**m(β) and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\widetilde{Q}}_{jm}(\beta )$$\end{document} . All their elements can be computed in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O({N}^{2}M)$$\end{document} operations, while the entire matrix A(β) in Eq. (15) can be computed in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O({NM}^{2})$$\end{document} operations. Moreover, solving the linear system (15) requires an additional \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O({M}^{3})$$\end{document} operations. Together, this means that the minimization of the functionals J({cm}, β) and Jincoh(β) requires \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O({N}^{2}M+{NM}^{2}+{M}^{3})$$\end{document} operations. Therefore, if M ≪ N, the resulting computational cost is much less than the standard number of operations \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O({N}^{3})$$\end{document} required to invert an N × N matrix.

Examples

To illustrate the possibilities of our parameter reconstruction algorithm, we apply it to the analysis of chimera states in model (5) with top-hat coupling and N = 2048. We choose the natural frequencies of all oscillators to be ω = 1. Then, for a fixed coupling range σ = 0.7, we vary the phase lag β in the range from 0.02 to 0.16, where stable chimera states can be observed. For each value of β, we first simulate system (5) for 10^5^ time units to allow it to reach statistical equilibrium, and then calculate the effective frequencies Ωj and local order parameters ζj using formulas (11) and (12) with T = 2000. Finally, we apply the above described parameter reconstruction algorithm with M = 10 spatial Fourier harmonics in formula (14). Figure 5 shows that in this way the value of β is reconstructed with an absolute accuracy of less than 0.0015. Similarly, the value of ω is reconstructed with an accuracy of less than 0.0008 (not shown).Fig. 5. Efficiency of the proposed method in reconstructing the phase lag *β.*The method was applied to the chimera state observed in model (5) with top-hat coupling. The dots show the values of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\beta }_{\min }$$\end{document} found as the global minimum of the function Jincoh(β) for different β. Parameters: N = 2048, ω = 1 and σ = 0.7.

In another round of simulations, we fix β = 0.1 and vary the coupling range σ from 0.6 to 0.78. For each value of σ, we repeat the same numerical protocol as above and check the accuracy with which our algorithm reconstructs the six leading Fourier coefficients cm in formula (14). Note that for the top-hat coupling function, these coefficients can be calculated analytically

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${c}_{m}=\left\{\begin{array}{ll}1/(2\pi ) \hfill & {{{\rm{for}}}}\,m=0, \hfill \\ \sin (\pi m\sigma )/({\pi }^{2}m\sigma ) & {{{\rm{for}}}}\,m=1,2,\ldots,\end{array}\right.$$\end{document}

so in Fig. 6 we compare the theoretical curves with several reconstructed parameter values (symbols), which turn out to be in excellent agreement with each other.Fig. 6. Efficiency of the proposed method in reconstructing the coupling function G(x).The six leading Fourier coefficients cm in formula (14) were recovered from the chimera state observed in model (5) with top-hat coupling. The curves show the theoretical values given by formula (16) and the symbols show the reconstructed values. Parameters: N = 2048, ω = 1 and β = 0.1.

Finally, in Figs. 7–9 we show that the proposed parameter reconstruction algorithm works equally well for other types of coupling functions in model (5). In particular, comparing the distribution of Fourier coefficients cm in Fig. 6, Fig. 8 and Fig. 9, we clearly see the possibility of distinguishing nonlocal couplings of the top-hat, exponential and cosine type. Moreover, from the value of the Fourier coefficient c1, we can uniquely determine the ranges and decay rates of the corresponding coupling functions, which confirms the reliability and efficiency of the proposed approach.Fig. 7. Other examples of application of the parameter reconstruction algorithm to chimera states in model (5) with N = 2048.a**–c Exponential coupling function with κ = 0.5, ω = 1, and α = π/2 − 0.1. d–f Cosine coupling function with A = 0.9, ω = 1, and α = π/2 − 0.15. In both cases, the coupling function G(x) was approximated by the ansatz (14) containing M = 10 spatial Fourier modes. In (c, f), the red/dark curves show reconstructed coupling functions, while the grey/light curves show the original coupling functions. Other notations are the same as in Fig. 4.Fig. 8. Reconstruction of the exponential coupling function from the chimera state observed in model (5).The Fourier coefficients, given by the formulas c0 = 1/(2π) and cm = (1 − (−1)^m^e^−*πκ*^)κ^2^/(π(1 − e^−π**κ^)(κ^2^ + m^2^)) for m≥1, are shown as curves, and the symbols indicate the reconstructed values. Other parameters: N = 2048, ω = 1 and β = 0.1.Fig. 9. Reconstruction of the cosine coupling function from the chimera state observed in model (5).The Fourier coefficients, given by the formulas c0 = 1/(2π), c1 = A/(2π) and cm = 0 for m≥2, are shown as curves, and the symbols indicate the reconstructed values. Other parameters: N = 2048, ω = 1 and β = 0.15.

It is important to note that although all the coupling functions G(x) considered above are normalized to unity, this is not a requirement of our method. We simply followed the established tradition in the literature on chimera states. Coupling functions with other normalizations can be analysed in the same way.

Parameter reconstruction algorithm with partial data

Even though the reconstruction algorithm described above requires performing calculations with only 2N variables Ωj and ζj, it can still become too resource-demanding if N is too large. This complication can be overcome by noting that statistical equilibrium relations (7)–(9) also remain valid if, instead of all values (xj, Ωj, ζj), j = 1, …, N, only a sufficiently large subset of them is used. Roughly speaking, from N points xj we can randomly select a smaller subset {xj: j ∈ S} with the number of elements #{S} < N. Then using the trapezoidal rule we can write an analogue of formula (10) that approximates the integral (21), albeit with worse accuracy than (10) (see Methods). This fact allows us to repeat all the steps of the above reconstruction algorithm, using only the indices j ∈ S in the linear regression, as well as in the definition of J({cm}, β) and Jincoh(β). Importantly, in this case, we need to calculate the effective frequencies and local order parameters only for j ∈ S. Moreover, the global order parameter Z(t) in (12) must be replaced with its “rarefied” analogue

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Z(t)=\frac{1}{\#\{S\}}{\sum }_{j\in S}{e}^{i{\theta }_{j}(t)}.$$\end{document}

Thus, the resulting reconstruction algorithm will use only observation of oscillators θj(t) with j ∈ S.

Figure 10 shows how such a modified algorithm works for a chimera state in model (5) with top-hat coupling and N = 8192 oscillators, if instead of all oscillators we randomly select 25% of them. Comparing Figs. 4 and 10, we see that our algorithm has good performance also with partial data.Fig. 10. Application of the parameter reconstruction algorithm to partial data.For a chimera state in model (5) with top-hat coupling and N = 8192 oscillators, we used 1996 randomly selected oscillators to reconstruct the parameters ω, β and G(x). The quantities shown in (a–c) are the same as in Fig. 4. Other model parameters are given in Fig. 2.

Time sampling and sensitivity to measurement noise

The only input data used in our algorithm are the time-averaged values of Ωj and ζj, which can be considered its advantage. Indeed, time averaging is a natural low-pass filter, so the algorithm is insensitive to the presence of noise in the phases θj, provided that the noise is unbiased (i.e., has a zero mean). On the other hand, for time averaging, phases do not need to be measured at evenly spaced time points (as is done in the map-based algorithms in refs. ^40,41^). More precisely, if we calculate the local order parameter ζj by formula

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\zeta }_{j}=\frac{1}{{N}_{T}}{\sum }_{k=1}^{{N}_{T}}{e}^{i{\theta }_{j}({t}_{k})}\frac{\overline{Z}({t}_{k})}{| Z({t}_{k})| },$$\end{document}

then we only need to worry that the number of points NT is large enough and that the points tk are uniformly distributed with respect to the oscillation period. As for calculating the effective frequencies Ωj, the effective formula is

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\Omega }_{j}=\frac{1}{{t}_{{N}_{T}}-{t}_{1}}{\sum }_{k=2}^{{N}_{T}}\arg {e}^{i({\theta }_{j}({t}_{k})-{\theta }_{j}({t}_{k-1}))}.$$\end{document}

Here, we need to satisfy the Nyquist criterion that the minimum of the intervals tk − tk−1 is less than the half-period of the corresponding oscillations. But the intervals tk − tk−1 do not have to be small, since we are not computing any time derivatives.

Chimera states with zero global order parameter

The definition of local order parameters exploits the fact that the global order parameter Z(t) does not vanish, and therefore its argument Z(t)/∣Z(t)∣ is well-defined. However, this condition is not satisfied for all types of chimera states. For example, it is violated for so called multi-headed chimera states^42^ with antiphase adjacent coherent regions. Fortunately, in this case, an alternative definition of local order parameters can be given. To do this, we first note that for patterns with antiphase regions, the Daido order parameter

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${Z}_{2}(t)=\frac{1}{N}{\sum }_{k=1}^{N}{e}^{2i{\theta }_{k}(t)}$$\end{document}

is usually non-vanishing. But the argument of this order parameter has twice the speed compared to the order parameter Z(t). Therefore, we define

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\widetilde{\zeta }}_{j}=\left\langle {e}^{i{\theta }_{j}(t)}\sqrt{\frac{{\overline{Z}}_{2}(t)}{| {Z}_{2}(t)| }}\right\rangle,$$\end{document}

where the complex square root is calculated in such a way that it varies as a continuous function of time. The value of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\widetilde{\zeta }}_{j}$$\end{document} may or may not coincide with the value of ζj defined in the Introduction, depending on the choice of the square root branch at t = 0. To avoid this ambiguity, we ultimately define

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\zeta }_{j}={\widetilde{\zeta }}_{j}\frac{\overline{\widetilde{Z}}}{| \widetilde{Z}| },\,\,{{{\rm{where}}}}\,\,\widetilde{Z}=\frac{1}{N}{\sum }_{k=1}^{N}{\widetilde{\zeta }}_{k}.$$\end{document}

The resulting formula gives the correct values of the local order parameters for multi-headed chimera states, using which the parameter reconstruction algorithm can be applied to these coherence-incoherence patterns as well.

Calculations with protophases

In applications, phase oscillator networks are used as simplified mathematical models to study the behaviour of weakly interacting limit cycle oscillators. The physical phase of a limit cycle oscillator is defined as a variable that increases uniformly from 0 to 2π during one cycle of its oscillations^10,11^. It is this phase that appears in the formulation of Kuramoto-type models obtained using the phase reduction method^10,13–15^. But the physical phase cannot be measured directly in an experiment. Therefore, it is usually extracted from the measured signal using Hilbert transform, projection onto a two-dimensional plane, marker appearance analysis, etc. In each of these approaches, the protophase is obtained, i.e., a quantity of the form ϕ = Φ(θ) where θ is the physical phase and Φ(θ) is an unknown 2π-periodic phase-protophase transformation^43^. This means that additional effort is usually required to reconstruct the function Φ(θ) along with other system parameters. Fortunately, for stationary coherence-incoherence patterns, such a step is not needed, provided that we only want to find the effective frequencies and local order parameters. This can be done as follows.

Suppose that we found a chimera state in Eq. (5), and instead of the physical phases θj(t) we have as available data the corresponding protophases ϕj(t) = Φ(θj(t)). Then, thanks to (32), the effective frequencies Ωj can be approximately calculated by

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\Omega }_{j}=\frac{1}{T}{\int }_{0}^{T}\frac{d{\phi }_{j}(t)}{dt}dt.$$\end{document}

Recalling that the coherent region of the chimera state corresponds to a plateau in the graph of Ωj versus xj, see Fig. 2, we can easily identify its frequency, for example, as a minimum \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\Omega }_{{{{\rm{c}}}}}={\min }_{j}{\Omega }_{j}$$\end{document} . Moreover, according to the continuum limit analysis presented in Methods, the frequency Ωc coincides with the frequency of the global order parameter Z(t), so that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Z(t)\approx {Z}_{0}{e}^{i{\Omega }_{{{{\rm{c}}}}}t}$$\end{document} .

Next, we consider the modified local order parameters

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\widehat{\zeta }}_{j}=\frac{1}{T}{\int }_{0}^{T}{e}^{i{\phi }_{j}(t)}{e}^{-i{\Omega }_{{{{\rm{c}}}}}t}dt.$$\end{document}

As explained in Methods, there are simple formulas (33), (34), independent of the explicit form of Φ(θ), which allow us to express the local order parameters ζj from their modified counterparts \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\widehat{\zeta }}_{j}$$\end{document} , provided that N ≫ 1 and the averaging time T is sufficiently long. Using these formulas as an approximation in the case of moderate size N and moderate averaging time T, we can write

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\zeta }_{j}=\frac{{\widehat{\zeta }}_{j}}{{\max }_{k}| {\widehat{\zeta }}_{k}| }\frac{\overline{\widehat{Z}}}{| \widehat{Z}| },\,\,{{{\rm{where}}}}\,\,\widehat{Z}=\frac{1}{N}{\sum }_{j=1}^{N}{\widehat{\zeta }}_{j}.$$\end{document}

Thus, all the input information needed for our parameter reconstruction algorithm can ultimately be obtained from the protophases ϕj(t) without knowing the explicit form of the phase-protophase transformation Φ(θ).

Discussion

The problem of model reconstruction for dynamical systems capable of exhibiting various patterns of synchrony and disorder has been a subject of research for a long time. Two general questions have usually been in focus. What is the coupling topology (i.e., network architecture) between individual agents in the system?^44^ And what form of coupling functions describes the interaction between these agents?^45^ Various methods have been proposed to answer these questions, including the finite-time mapping approach^40,41^, and random phase resetting method^46^, fixed points analysis^44^, and kernel density estimation^47^. In^48^ it was shown that for pulse-coupled oscillators the network topology can be reconstructed from spiking sequences. In addition, some methods were inspired by data assimilation approach, including dynamical Bayesian inference^49^, maximum likelihood estimation combined with multiple shooting^50,51^, and an ensemble Kalman filter with state space augmentation^52^. Each of these methods has its advantages and disadvantages, but they all become increasingly complex and resource-demanding as the system size increases, so they are usually applied to systems consisting of several dozen or hundreds of individual agents. On the contrary, in this paper we described a method that is much better suited for similar inverse problems, but in the case of large-size systems. We showed how it can be used to noninvasively reconstruct the parameters of the Kuramoto-Battogtokh model from a single observation of a chimera state transformed into a small dataset of time-averaged quantities. Although we have only examined this special example in detail, it seems that the method can also be generalized to a broader class of networks, including two-dimensional arrays with nonlocal coupling^53,54^, as well as networks with heterogeneous coupling coefficients^55^ and heterogeneous natural frequencies^56^. (Several examples of how this can be done are described in the section Methods below). Moreover, given the existing continuum limit theory for travelling and breathing chimera states^38,39^, the SER-based approach can also be extended to analyze these nonstationary coherence-incoherence patterns, which remained outside the scope of the present work. Furthermore, using our method, it is potentially possible to consider phase oscillator models with higher-harmonics^57^ and higher-order interactions^58,59^, although in this case other types of time-averaged observables and multiple observations would certainly be required.

From a more general perspective, analogues of statistical equilibrium relations can be written not only for phase oscillator networks, but also for many other systems, in particular for those that can be considered using the Ott-Antonsen theory or the self-consistency approach (see Methods). This suggests that the proposed model reconstruction scheme can be applied with appropriate modifications to neural networks (e.g., those consisting of theta neurons^60,61^ or quadratic integrate-and-fire neurons^62,63^) and Kuramoto-type models for power grids^64,65^.

Finally, we note that the knowledge of the existence of statistical equilibrium relations in a given system can be useful in itself. For example, it can be a natural clue to the lower-dimensional manifold or collective variables representing its long-term dynamics^4,5,66^. Such information, in turn, can facilitate or refine the development of a data-driven model reconstruction algorithms, reducing their memory usage and increasing their computational efficiency.

Methods

Derivation of SERs for nonlocally coupled oscillators

In this section, we will show how to derive statistical equilibrium relations (7)–(9) for the Kuramoto-Battogtokh system (5). Importantly, we do not make any special assumptions about the natural frequency ω, the phase lag α, or the coupling function G(x). But we assume that some stationary coherence-incoherence pattern arises in system (5). Roughly speaking, we carry out the following steps. First, we write the continuity equation corresponding to system (5) in the large-N limit and the general self-consistent ansatz of its stationary solutions. Using this ansatz, we obtain formulas (29), (30) and (31), which have the form of statistical equilibrium relations. Finally, to complete the definition of relations (29) and (30), we write down the Riemann sum approximation of formula (21), which gives the analytic expression (26) for the mean field ξj. Note that the above derivation scheme can be easily generalized to other types of coupled oscillator systems, for which the continuity equation and a self-consistent representation of its stationary solutions can be written.

First, we rewrite the Kuramoto-Battogtokh system (5) in the form (1), (2):

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{d{\theta }_{j}}{dt}=\omega -{{{\rm{Im}}}}\left({\overline{W}}_{j}{e}^{i{\theta }_{j}}{e}^{i\alpha }\right),$$\end{document}

where

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${W}_{j}=\frac{2\pi }{N}{\sum }_{k=1}^{N}G({x}_{j}-{x}_{k}){e}^{i{\theta }_{k}}$$\end{document}

and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\overline{W}}_{j}$$\end{document} denotes the complex-conjugate of Wj.

In the large-N limit, called also the continuum limit, the state of the system (5) can be represented by a probability density function ρ(θ, x, t) with x ∈ [ − π, π]. Then, the dynamics of ρ is determined by a continuity equation

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{\partial \rho }{\partial t} +\frac{\partial }{\partial \theta }\left(\left[\omega -{\int }_{-\pi }^{\pi }{\int }_{-\pi }^{\pi }G(x-x' )\sin (\theta -\theta'+\alpha )\right.\right.\\ \left.\left. \times \rho (\theta',x',t)d\theta' dx' \right]\rho \right)=0.$$\end{document}

It is known^27,67^ that the chimera patterns shown above behave like statistical equilibria, namely, each of them has a time-independent probability density in an appropriate corotating frame. In addition, it is known that these probability densities lie on a special Ott-Antonsen manifold^68,69^ consisting of functions of the form

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho (\theta,x,t)=\frac{1}{2\pi }\left(1+{\sum }_{n=1}^{\infty }\left[{\overline{z}}^{n}(x,t){e}^{in\theta }+{z}^{n}(x,t){e}^{-in\theta }\right]\right),$$\end{document}

where z(x, t) satisfies the integro-differential equation

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{dz}{dt}=i\omega z+\frac{1}{2}{e}^{-i\alpha }{{{\mathcal{G}}}}z-\frac{1}{2}{e}^{i\alpha }{z}^{2}{{{\mathcal{G}}}}\overline{z}$$\end{document}

with the integral operator

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$({{{\mathcal{G}}}}z)(x,t)={\int }_{-\pi }^{\pi }G(x-x' )z(x',t)dx',$$\end{document}

and moreover ∣z(x, t)∣≤1 for all x ∈ [ − π, π].

In ref. ^67^ it was shown that every stationary chimera state in the Kuramoto-Battogtokh system (5) corresponds to a rotating wave solution of Eq. (19) given by the formula

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z(x,t)=a(x){e}^{i\Omega t}.$$\end{document}

Inserting this ansatz into Eq. (19) and denoting

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$w(x)=\frac{1}{\omega -\Omega }{\int }_{-\pi }^{\pi }G(x-x')a(x')dx',$$\end{document}

we obtain

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${e}^{-i\beta }\overline{w}(x){a}^{2}(x)-2a(x)+{e}^{i\beta }w(x)=0.$$\end{document}

These equations allow us to justify the statistical equilibrium relations (7)–(9).

Note that the above ansatz for ρ(θ, x, t) implies

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\int }_{-\pi }^{\pi }{\int }_{-\pi }^{\pi }\rho (\theta,x,t)d\theta \,dx=2\pi,\\ {\int }_{-\pi }^{\pi }{e}^{ik\theta }\rho (\theta,x,t)d\theta={a}^{k}(x){e}^{ik\Omega t}\,\,{{{\rm{for}}}}\,\,k\in {\mathbb{N}}.$$\end{document}

Therefore, in the large-N limit, the definition of the global order parameter (6) can be written in the form

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Z(t)=\frac{1}{2\pi }{\int }_{-\pi }^{\pi }{\int }_{-\pi }^{\pi }{e}^{i\theta }\rho (\theta,x,t)d\theta \,dx={Z}_{0}{e}^{i\Omega t},$$\end{document}

where

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${Z}_{0}=\frac{1}{2\pi }{\int }_{-\pi }^{\pi }a(x)dx.$$\end{document}

Similarly, using the ergodicity property, we obtain

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\zeta }_{j}={\int }_{-\pi }^{\pi }{e}^{i\theta }\frac{{\overline{Z}}_{0}}{| {Z}_{0}| }{e}^{-i\Omega t}\rho (\theta,{x}_{j},t)d\theta=a({x}_{j})\frac{{\overline{Z}}_{0}}{| {Z}_{0}| },$$\end{document}

and therefore

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{1}{N}{\sum }_{j=1}^{N}{\zeta }_{j}=| {Z}_{0}| \,\,{{{\rm{as}}}}\,\,N\to \infty .$$\end{document}

Let us denote aj = a(xj) and wj = w(xj), then using the trapezoidal rule we write an approximate version of the definition (21)

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${w}_{j}=\frac{1}{\omega -\Omega }{\sum }_{k=1}^{N}G({x}_{j}-{x}_{k}){a}_{k}\frac{{x}_{k+1}-{x}_{k-1}}{2},$$\end{document}

where due to the periodicity of the variable x we assume x2 − x0 = 2π + x2 − xN and xN+1 − xN−1 = 2π + x1 − xN−1. Multiplying this by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${e}^{i\beta }{\overline{Z}}_{0}/| {Z}_{0}|$$\end{document} and defining

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\xi }_{j}={w}_{j}{e}^{i\beta }\frac{{\overline{Z}}_{0}}{| {Z}_{0}| },$$\end{document}

we obtain formula (10). On the other hand, from Eq. (22) it follows

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\overline{\xi }}_{j}{\zeta }_{j}^{2}-2{\zeta }_{j}+{\xi }_{j}=0.$$\end{document}

Proposition. Suppose that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\xi }_{j}\in {\mathbb{C}}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\zeta }_{j}\in {\mathbb{C}}$$\end{document} satisfy equation (28), then

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\xi }_{j}=\frac{2{\zeta }_{j}}{1+| {\zeta }_{j}{| }^{2}}\,\,{{{\rm{for}}}}\,\,| {\zeta }_{j}| \ne 1$$\end{document}

and

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\rm{Re}}}}({\xi }_{j}{\overline{\zeta }}_{j})=1\,\,{{{\rm{for}}}}\,\,| {\zeta }_{j}|=1.$$\end{document}

Moreover, for all values of ∣ζj∣ we have

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\rm{Re}}}}\left(\frac{{\xi }_{j}}{{\zeta }_{j}}\right)=\frac{2}{1+| {\zeta }_{j}{| }^{2}}$$\end{document}

and

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\rm{Re}}}}\left({\xi }_{j}{\overline{\zeta }}_{j}\right)=\frac{2| {\zeta }_{j}{| }^{2}}{1+| {\zeta }_{j}{| }^{2}}.$$\end{document}

Proof: The complex conjugate of Eq. (28) reads

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\xi }_{j}{\overline{\zeta }}_{j}^{2}-2{\overline{\zeta }}_{j}+{\overline{\xi }}_{j}=0,$$\end{document}

or equivalently \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\overline{\xi }}_{j}=2{\overline{\zeta }}_{j}-{\xi }_{j}{\overline{\zeta }}_{j}^{2}$$\end{document} . Inserting this into Eq. (28), we obtain

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2| {\zeta }_{j}{| }^{2}{\zeta }_{j}-{\xi }_{j}| {\zeta }_{j}{| }^{4}-2{\zeta }_{j}+{\xi }_{j}=0,$$\end{document}

or

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\xi }_{j}(1-| {\zeta }_{j}{| }^{4})=2{\zeta }_{j}(1-| {\zeta }_{j}{| }^{2}).$$\end{document}

If ∣ζj∣ ≠ 1, this yields (29).

On the other hand, if ∣ζj∣ = 1, then \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\overline{\zeta }}_{j}=1/{\zeta }_{j}$$\end{document} , and therefore dividing (28) by ζj, we obtain

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\overline{\xi }}_{j}{\zeta }_{j}-2+{\xi }_{j}{\overline{\zeta }}_{j}=0$$\end{document}

what is equivalent to (30).■

When N ≫ 1, formula (18) can be rewritten using (23), (26) and (27). This gives

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${W}_{j} = \frac{2\pi }{N}{\sum }_{k=1}^{N}G({x}_{j}-{x}_{k})a({x}_{k}){e}^{i\Omega t}=(\omega -\Omega ){w}_{j}{e}^{i\Omega t}\\ = (\omega -\Omega ){\xi }_{j}{e}^{-i\beta }\frac{{Z}_{0}}{| {Z}_{0}| }{e}^{i\Omega t}.$$\end{document}

Inserting this into (17), we obtain

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{d{\theta }_{j}}{dt}=\omega -(\omega -\Omega ){{{\rm{Im}}}}\left(i{\overline{\xi }}_{j}\frac{{\overline{Z}}_{0}}{| {Z}_{0}| }{e}^{-i\Omega t}{e}^{i{\theta }_{j}}\right).$$\end{document}

Therefore, thanks to the ergodicity property, identity (24) and the above Proposition, we have

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\Omega }_{j} =\left\langle \frac{d{\theta }_{j}}{dt}\right\rangle=\omega -(\omega -\Omega ){{{\rm{Im}}}}\left(i{\overline{\xi }}_{j}{\zeta }_{j}\right)\\ =\omega -(\omega -\Omega ){{{\rm{Re}}}}\left({\xi }_{j}{\overline{\zeta }}_{j}\right)=\omega -(\omega -\Omega )\frac{2| {\zeta }_{j}| ^{2}}{1+| {\zeta }_{j}| ^{2}},$$\end{document}

that is equivalent to the statistical equilibrium relation (7).

Formulas involving protophases

In addition to the identities we have already obtained, we can also derive similar formulas containing an arbitrary 2π-periodic function Φ(θ) that maps the interval [0, 2π) onto itself. More precisely, suppose that ϕj(t) = Φ(θj(t)), then

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left\langle \frac{d{\phi }_{j}}{dt}\right\rangle=\left\langle \frac{d{\theta }_{j}}{dt}\right\rangle={\Omega }_{j}.$$\end{document}

Roughly speaking, both of the above time-averages determine the same rotation number, which is invariant with respect to the form of the function Φ(θ).

Next, we consider a modified local order parameter

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\widehat{\zeta }}_{j}=\left\langle {e}^{i{\phi }_{j}(t)}{e}^{-i\Omega t}\right\rangle .$$\end{document}

We note that for any 2π-periodic function Φ(θ), it holds

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${e}^{i\Phi (\theta )}={\sum }_{k=-\infty }^{\infty }{\varphi }_{k}{e}^{ik\theta }\,\,{{{\rm{with\; some}}}}\,\,{\varphi }_{k}\in {\mathbb{C}},$$\end{document}

and therefore

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\widehat{\zeta }}_{j}=\left\langle {e}^{i\Phi ({\theta }_{j}(t))}{e}^{-i\Omega t}\right\rangle={\sum }_{k=-\infty }^{\infty }{\varphi }_{k}\left\langle {e}^{ik{\theta }_{j}(t)}{e}^{-i\Omega t}\right\rangle .$$\end{document}

Using the ergodicity property and formulas (23), we calculate

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left\langle {e}^{ik{\theta }_{j}(t)}{e}^{-i\Omega t}\right\rangle= \left\langle {\int }_{\!\!\!\!-\pi }^{\pi }{e}^{ik\theta }{e}^{-i\Omega t}\rho (\theta,{x}_{j},t)d\theta \right\rangle \,=\left\langle {a}^{k}({x}_{j}){e}^{i(k-1)\Omega t}\right\rangle \\= a({x}_{j}){\delta }_{k1}=\frac{{Z}_{0}}{| {Z}_{0}| }{\zeta }_{j}{\delta }_{k1}\,\,{{{\rm{for}}}}\,\,k\in {\mathbb{N}},$$\end{document}

and similarly

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left\langle {e}^{ik{\theta }_{j}(t)}{e}^{-i\Omega t}\right\rangle=0\,\,{{{\rm{for}}}}\,\,k=0,-1,-2,\ldots .$$\end{document}

This means

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\widehat{\zeta }}_{j}={\varphi }_{1}\frac{{Z}_{0}}{| {Z}_{0}| }{\zeta }_{j}.$$\end{document}

Defining

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\widehat{Z}=\frac{1}{N}{\sum }_{j=1}^{N}{\widehat{\zeta }}_{j}$$\end{document}

and using (25), we obtain \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\widehat{Z}={\varphi }_{1}{Z}_{0}$$\end{document} and

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\widehat{\zeta }}_{j}\frac{\overline{\widehat{Z}}}{| \widehat{Z}| }=| {\varphi }_{1}| {\zeta }_{j}.$$\end{document}

Given that for each coherent oscillator we have ∣ζj∣ = 1, and for each incoherent oscillator ∣ζj∣ < 1, we write \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$| {\varphi }_{1}|={\max }_{j}| {\widehat{\zeta }}_{j}|$$\end{document} . Thus, we obtain the formulas

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\zeta }_{j}=\frac{{\widehat{\zeta }}_{j}}{{\max }_{k}| {\widehat{\zeta }}_{k}| }\frac{\overline{\widehat{Z}}}{| \widehat{Z}| },$$\end{document}

which are valid for every stationary coherence-incoherence pattern in Eq. (5) with N ≫ 1, regardless of the choice of the function Φ(θ).

Towards other applications

To demonstrate the versatility of the proposed method, we show how it can be applied to other partially synchronized patterns and other coupled oscillator networks. At the same time, we will show that statistical equilibrium relations may look different in different models, and therefore the parameter reconstruction algorithm must be adapted to each specific case individually.

We recall that a plethora of thermodynamic limit results are known for heterogeneous phase oscillator networks, in particular for a network of the form

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{d{\theta }_{j}}{dt}={\omega }_{j}-\frac{1}{N}{\sum }_{k=1}^{N}{q}_{j}{q}_{k}\sin ({\theta }_{j}-{\theta }_{k}+\alpha ),$$\end{document}

where the natural frequencies ωj and coupling strengths qj > 0 are drawn randomly and independently from distributions Hω(ω) and Hq(q), respectively. In the large-N limit, this model allows a universal representation of all stationary partially synchronized states, which is given by the probability density^55,70,71^

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho (\theta,\omega,q,t)= \frac{{H}_{\omega }(\omega ){H}_{q}(q)}{2\pi } \\ \times \left(1+{\sum }_{n=1}^{\infty }\left[{\overline{a}}^{n}(\omega,q,t){e}^{in\theta }+{a}^{n}(\omega,q,t){e}^{-in\theta }\right]\right),$$\end{document}

where

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a(\omega,q,t)=h\left(\frac{\omega -\Omega }{pq}\right){e}^{i\Omega t}$$\end{document}

and

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h(s)=\left\{\begin{array}{ll}s-i\sqrt{1-{s}^{2}} & {{{\rm{for}}}}\,| s| \le 1,\\ (1-\sqrt{1-{s}^{-2}})s & {{{\rm{for}}}}\,| s| > 1,\end{array}\right.$$\end{document}

and the parameters p > 0 and  − ∞ < Ω < ∞ satisfy the self-consistency equation

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p=i{e}^{-i\alpha }\int^{\infty }_{-\infty }d\omega \int^{\infty }_{0}{H}_{\omega }(\omega ){H}_{q}(q)h\left(\frac{\omega -\Omega }{pq}\right)q\,dq.$$\end{document}

It is easy to see that the definitions of the effective frequencies Ωj and the local order parameters ζj given in the Introduction also make sense for model (35). Then, rewriting (35) in the form

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{d{\theta }_{j}}{dt}={\omega }_{j}-{{{\rm{Im}}}}(\overline{W}{e}^{i{\theta }_{j}}{e}^{i\alpha }),$$\end{document}

where

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$W(t)=\frac{1}{N}{\sum }_{k=1}^{N}{q}_{k}{e}^{i{\theta }_{k}(t)},$$\end{document}

and using the above formula for ρ(θ, ω, q, t), two sets of statistical equilibrium relations can be obtained^56^

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\zeta }_{k}=\frac{{\overline{{{{\mathcal{Z}}}}}}_{0}}{| {{{{\mathcal{Z}}}}}_{0}| }h\left(\frac{{\omega }_{k}-\Omega }{p{q}_{k}}\right),$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\Omega }_{k}=\Omega+p{q}_{k}Q\left(\frac{{\omega }_{k}-\Omega }{p{q}_{k}}\right),$$\end{document}

where

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{{\mathcal{Z}}}}}_{0}=\int\limits^{\infty }_{-\infty }d\omega \int\limits^{\infty }_{0}{H}_{\omega }(\omega ){H}_{q}(q)h\left(\frac{\omega -\Omega }{pq}\right)dq$$\end{document}

and

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Q(s)=\left\{\begin{array}{ll}0 \hfill & {{{\rm{for}}}}\,| s| \le 1,\\ s\sqrt{1-{s}^{-2}} & {{{\rm{for}}}}\,| s| > 1.\end{array}\right.$$\end{document}

According to their derivation, relations (38) and (39) are exact for N → ∞ and independent of the distributions Hω(ω) and Hq(q). So, from the perspective of the statistical physics, it is logical to assume that they remain valid with a small error also for finite but large N.

SERs for globally coupled oscillators

Our first example is the paradigmatic Kuramoto-Sakaguchi model^72^

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{d{\theta }_{j}}{dt}={\omega }_{j}-\frac{K}{N}{\sum }_{k=1}^{N}\sin ({\theta }_{j}-{\theta }_{k}+\alpha ),$$\end{document}

which describes the dynamics of all-to-all coupled phase oscillators with different frequencies ωj. In this case \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${H}_{q}(q)=\delta (q-\sqrt{K})$$\end{document} , and therefore the self-consistency equation (37) and formula (40) can be written in the form

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\widetilde{p}=iK{e}^{-i\alpha }{{{{\mathcal{Z}}}}}_{0},$$\end{document}

where

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\widetilde{p}=p\sqrt{K}\,\,{{{\rm{and}}}}\,\,{{{{\mathcal{Z}}}}}_{0}={\int }_{-\infty }^{\infty }{H}_{\omega }(\omega )h\left(\frac{\omega -\Omega }{\widetilde{p}}\right)d\omega .$$\end{document}

Accordingly, the statistical equilibrium relations (38) and (39) take the form

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\zeta }_{k}=i{e}^{-i\alpha }h\left(\frac{{\omega }_{k}-\Omega }{\widetilde{p}}\right),$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\Omega }_{k}=\Omega+\widetilde{p}Q\left(\frac{{\omega }_{k}-\Omega }{\widetilde{p}}\right).$$\end{document}

In^56^ it was shown that using relations (42)–(44) all parameters in model (41) can be reconstructed from the observed values of ζk and Ωk and from the time-averaged magnitude of the global order parameter Z(t). (Note that formula (43) differs from its counterpart in ref. ^56^ due to the presence of an additional prefactor  − i in the definition of ζk used in ref. ^56^). Roughly speaking, since the function h(s) satisfies the quadratic equation h^2^ − 2s**h + 1 = 0 and relations (43) hold, the parameter α can be found from an explicitly solvable minimization problem

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{\pi }{2}-\alpha={\arg\!\min }_{\beta \in [0,\pi ]}{\sum }_{k=1}^{N}{\left[{{{\rm{Im}}}}\left(\frac{{\zeta }_{k}^{2}{e}^{-2i\beta }+1}{2{\zeta }_{k}{e}^{-i\beta }}\right)\right]}^{2},$$\end{document}

and moreover

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\widetilde{s}}_{k}=\frac{{\omega }_{k}-\Omega }{\widetilde{p}}={{{\rm{Re}}}}\left(\frac{{\zeta }_{k}^{2}{e}^{2i\alpha }-1}{2i{\zeta }_{k}{e}^{i\alpha }}\right).$$\end{document}

Once the ratios \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\widetilde{s}}_{k}$$\end{document} are known, the parameters Ω and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\widetilde{p}$$\end{document} can be determined by linear regression from the set of equations (44). Then, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\omega }_{k}=\Omega+\widetilde{p}{\widetilde{s}}_{k}$$\end{document} . And finally, from (42) we get the coupling strength \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K=\widetilde{p}/| {{{{\mathcal{Z}}}}}_{0}|$$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$| {{{{\mathcal{Z}}}}}_{0}|$$\end{document} is replaced by the time-averaged magnitude of the global order parameter Z(t) defined by (6).

SERs for random networks

Partially synchronized patterns can also be found in random networks of identical coupled phase oscillators. In their simplest form, such networks are defined by equations^9^

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{d{\theta }_{j}}{dt}={\omega }_{0}-\frac{K}{{d}_{{{{\rm{mean}}}}}}{\sum }_{k=1}^{N}{a}_{jk}\sin ({\theta }_{j}-{\theta }_{k}+\alpha ),$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\{{a}_{jk}\}}_{j,k=1}^{N}$$\end{document} denotes an adjacency matrix such that aj**k = 1 if there is a link between the jth and kth oscillators, and aj**k = 0 otherwise. The number

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${d}_{j}={\sum }_{k=1}^{N}{a}_{jk}$$\end{document}

is called the degree of the oscillator j, and

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${d}_{{{{\rm{mean}}}}}=\frac{1}{N}{\sum }_{j=1}^{N}{d}_{j}$$\end{document}

denotes the mean degree. For densely connected random networks, their dynamics can be described using the annealed approximation^9,73^

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{d{\theta }_{j}}{dt}={\omega }_{0}-\frac{K}{N}{\sum }_{k=1}^{N}\frac{{d}_{j}}{{d}_{{{{\rm{mean}}}}}}\frac{{d}_{k}}{{d}_{{{{\rm{mean}}}}}}\sin ({\theta }_{j}-{\theta }_{k}+\alpha ),$$\end{document}

or equivalently

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{d{\theta }_{j}}{dt}={\omega }_{0}-\frac{1}{N}{\sum }_{k=1}^{N}{q}_{j}{q}_{k}\sin ({\theta }_{j}-{\theta }_{k}+\alpha ),$$\end{document}

where

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${q}_{j}=\frac{\sqrt{K}{d}_{j}}{{d}_{{{{\rm{mean}}}}}}.$$\end{document}

Now model (46) looks like a special case of the auxiliary model (35) with Hω(ω) = δ(ω − ω0). The corresponding statistical equilibrium relations have the form

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\zeta }_{k}=\frac{{\overline{{{{\mathcal{Z}}}}}}_{0}}{| {{{{\mathcal{Z}}}}}_{0}| }h({s}_{k}),$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\Omega }_{k}=\Omega+({\omega }_{0}-\Omega ){s}_{k}Q({s}_{k}),$$\end{document}

where

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${s}_{k}=\frac{{\omega }_{0}-\Omega }{p{q}_{k}}.$$\end{document}

Moreover, the self-consistency equation (37) also simplifies to

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p=i{e}^{-i\alpha }{{{{\mathcal{W}}}}}_{0},$$\end{document}

where

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{{\mathcal{W}}}}}_{0}={\int }_{0}^{\infty }{H}_{q}(q)h\left(\frac{{\omega }_{0}-\Omega }{pq}\right)q\,dq.$$\end{document}

Next, we will demonstrate that the parameters ω0, α and qj in Eq. (46) can be found from the observed values of ζj and Ωj. First, we note that the definitions (36) and (40) imply that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\rm{Im}}}}({{{{\mathcal{Z}}}}}_{0})\le 0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-\pi \le \arg {{{{\mathcal{Z}}}}}_{0}\le 0$$\end{document} . Recalling that h(s) satisfies the quadratic equation h^2^ − 2s**h + 1 = 0 and using relations (47), we write

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\arg {{{{\mathcal{Z}}}}}_{0}={\phi }_{ * } := {\arg\!\min }_{\phi \in [-\pi,0]}{\sum }_{k=1}^{N}{\left[{{{\rm{Im}}}}\left(\frac{{\zeta }_{k}^{2}{e}^{2i\phi }+1}{2{\zeta }_{k}{e}^{i\phi }}\right)\right]}^{2},$$\end{document}

and

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${s}_{k}={{{\rm{Re}}}}\left(\frac{{\zeta }_{k}^{2}{e}^{2i{\phi }_{ * }}+1}{2{\zeta }_{k}{e}^{i{\phi }_{ * }}}\right).$$\end{document}

Then, the parameter Ω and the difference ω0 − Ω can be determined by linear regression from equations (48). This allows us to find p**qk = (ω0 − Ω)/sk. Taking into account that formula (50) defines the mean values of the product q**h((ω0 − Ω)/p**q), we conclude

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p{{{{\mathcal{W}}}}}_{0}\approx {{{\mathcal{S}}}} := \frac{1}{N}{\sum }_{k=1}^{N}p{q}_{k}h({s}_{k})=\frac{1}{N}{\sum }_{k=1}^{N}\frac{{\omega }_{0}-\Omega }{{s}_{k}}h({s}_{k}).$$\end{document}

Substituting this into the self-consistency equation (49), we finally get

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p=\sqrt{{{{\mathcal{S}}}}},\,\alpha=\arg (i{{{\mathcal{S}}}}),\,\,{{{\rm{and}}}}\,\,{q}_{j}=\frac{{\omega }_{0}-\Omega }{p{s}_{k}}.$$\end{document}

Note that using the above approach, we can only find the values of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${q}_{k}=\sqrt{K}{d}_{k}/{d}_{{{{\rm{mean}}}}}$$\end{document} , but we cannot separate them into the coupling strength K and the degree dk. In any case, even this information is enough to decide, which oscillators are more connected and which are less connected, and also what type of distribution Hq(q) characterizes the network.

Remark about the partial data case

The statistical equilibrium relations (38) and (39) express the relationship between local observables and local parameters. Therefore, they remain valid also for partial observations of ζk and Ωk. However, to implement the above parameter reconstruction algorithms for models (41) and (46), it is necessary to know the value of a global order parameter such as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{{\mathcal{Z}}}}}_{0}$$\end{document} or \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{{\mathcal{W}}}}}_{0}$$\end{document} , see (40) and (50). Fortunately, both of these quantities are defined as averages across the entire network. Therefore, if measurements are only available for a part of the oscillators, their parameters can still be found (albeit with less precision), provided that these oscillators are uniformly distributed across the network. In this case, the global averages should be calculated as averages over the available oscillators.

Supplementary information

Supplementary Information Description of Additional Supplementary Files Supplementary Data 1 Transparent Peer Review file