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

Domain: physics
Published: 2026-02-27
Editorial outcome: major_revision

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.

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 equations1. 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 dynamics2, automated adaptive model inference3, data-driven discovery of intrinsic lower-dimensional dynamics4,5, machine learning techniques based on reduced order models6, and others7. 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 biology8,9. They play a key role in the study of synchronization phenomena10–12 and have a direct connection to more complex real-world models through the standard phase reduction procedure10,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 sensing16,17, circadian rhythm generators18,19, metachronal waves in cilia carpets20, brain disorders21,22, and other physiological processes related to synchrony23. 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 θ. Each oscillator has a label p_j ∈ R^Np 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

dθ_j / dt = F(θ_j, p_j, W_j), (1)

where

W_j = 1/N * Σ_{k=1}^N Q(p_j, p_k, θ_k) (2)

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 p_j 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

W_j → W[ρ](p_j) = ∫_{R^Np} ∫_{-π}^π Q(p_j, p, θ) ρ(θ, p, t) dθ dp

and write a nonlinear integro-differential continuity equation

∂ρ/∂t + ∂/∂θ (F(θ, p, W[ρ])ρ) = 0 (3)

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 analysis10. 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

Ω_j = ⟨dθ_j / dt⟩,

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

Ω_j = 1/π ∫_{-π}^π ρ_se(θ, p_j, t) F(θ, p_j, W[ρ_se](p_j)) dθ, (4)

where the use of the conditional probability density ρ_se(θ, p_j, t)/g(p_j) is due to the fact that we are considering an oscillator with label p_j. Formula (4) gives an algebraic relationship between the time-averaged observable Ω_j and the system parameters {p_j}. 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 oscillators24. It is famous as a prototype system for chimera states25–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.

Model
We consider a ring of N nonlocally coupled identical phase oscillators

dθ_j / dt = ω - 2π/N * Σ_{k=1}^N G(x_j - x_k) sin(θ_j - θ_k + α). (5)

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 x_j ∈ [−π, π] and the nonlocal interaction between oscillators is determined by a scalar symmetric 2π-periodic coupling function G(x). More precisely, the positions x_j are assumed to be uniformly distributed on the interval [−π, π], although in most numerical examples below we will use a special deterministic choice x_j = −π + 2πj/N.

Pattern
It is well-known24,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 function24

G(x) = κ / (2(1 - e^{-πκ})) * e^{-κ arccos(cos x)},

cosine function25

G(x) = 1/(2π) * (1 + A cos x),

top-hat function28

G(x) = 1/(4πσ) * (1 + 1/σ * arccos(cos x) |arccos(cos x)| ≤ πσ),

and many others27. (Note that due to the periodic boundary conditions in model (5), above we used the expression arccos(cos x), which is equal to |x| in the interval [−π, π], and is defined as 2π - |x| in the rest of the real line.) But later their existence was confirmed experimentally in systems of chemical Belousov-Zhabotinsky oscillators29,30 and in systems of electrochemical oscillators31. In addition, their similarity to dynamic patterns in various biological systems has been established. These include synchronization patterns of elastic cilia32,33, collective states of coupled inner-ear hair cells34, and epileptic seizures35. 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 oscillators36, 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 dθ_j(t)/dt into Eq. (5) and solve the resulting system with respect to the unknown parameters37. 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

ζ_j = ⟨e^{iθ_j(t)}⟩ / ⟨|Z(t)|⟩,

where

Z(t) = 1/N * Σ_{k=1}^N e^{iθ_k(t)} (6)

is the global order parameter of all oscillators and Z(t) 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)

(ω - Ω_j) |ζ_j|^2 = (ω - Ω_j) / (1 + |ζ_j|^2), (7)

Re(ζ_j ξ_j) = |ζ_j|^2 / (1 + |ζ_j|^2), (8)

ξ_j = 2ζ_j / (1 + |ζ_j|^2) for |ζ_j| < 1, (9)

where

ξ_j = e^{iβ} / (ω - Ω_j) * Σ_{k=1}^N G(x_j - x_k) ζ_k e^{i(x_{k+1} - x_{k-1})/2} (10)

and

β = π/2 - α.

(Note that in (10) the notations x_{N+1} - x_N = 2π + x_1 - x_N and x_1 - x_{N-1} = 2π + x_1 - x_{N-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 states38 and breathing chimera states39 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 motion38.

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 x_j, 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

Ω_j(T) = 1/T ∫_0^T dθ_j(t)/dt dt, (11)

ζ_j(T) = 1/T ∫_0^T e^{iθ_j(t)} / |Z(t)| dt, (12)

where

Ω(T) = 1/T ∫_0^T 1/|Z(t)| d|Z(t)|/dt dt, (13)

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

δ_1(T) = 1/N Σ_{j=1}^N |(ω - Ω_j(T)) |ζ_j(T)|^2 - (ω - Ω_j(T)) / (1 + |ζ_j(T)|^2)|,

δ_2(T) = 1/N Σ_{j=1}^N |Re(ζ_j(T) ξ_j(T)) - |ζ_j(T)|^2 / (1 + |ζ_j(T)|^2)|,

δ_3(T) = 1/N_c Σ_{j: |ζ_j(T)| < 1 - 1/√N} |ξ_j(T) - 2ζ_j(T) / (1 + |ζ_j(T)|^2)|,

where N_c is the number of indices j satisfying the inequality |ζ_j(T)| < 1 - 1/√N. 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.

Parameter reconstruction algorithm
Suppose that the statistical equilibrium relations (7)–(9) are satisfied for the observed variables Ω_j and ζ_j in model (5). How can we reconstruct the parameters ω, β and G(x)? It is easy to see that the values of ω and ω - Ω_j can be found from the statistical equilibrium relations (7) written in the form

Ω_j = ω + (Ω_j - ω) η_j with η_j = 2|ζ_j|^2 / (1 + |ζ_j|^2).

Indeed, standard linear regression yields, see Fig. 4a, ω = (S_Ω S_ηη - S_η S_Ωη) / (S_η^2 - S_ηη), (ω - Ω_j) = (S_Ω S_η - S_Ωη) / (S_η^2 - S_ηη),

where

S_η = 1/N Σ η_j, S_Ω = 1/N Σ Ω_j, S_ηη = 1/N Σ η_j^2, S_Ωη = 1/N Σ Ω_j η_j.

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

Re(ξ_j) = Σ_{k=1}^N G(x_j - x_k) Re(e^{iβ} ζ_k) / (ω - Ω_j) |ζ_j|^2

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

G(x) = Σ_{m=0}^M c_m q_m(x) with q_m(x) = cos(mx). (14)

It is not known a priori, so M in Eq. (14) is chosen empirically. According to (8) we can expect that the vector ({c_m}, β) is the minimizer of the functional

J({c_m}, β) = 1/N Σ_{j=1}^N |1/(1 + |ζ_j|^2) - Σ_{m=0}^M c_m Q_{jm}(β)|^2,

where

Q_{jm}(β) = Σ_{k=1}^N q_m(x_j - x_k) Re(e^{iβ} ζ_k) / (ω - Ω_j) |ζ_j|^2.

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

∂_{c_n} J({c_m}, β) = - 2/N Σ_{j=1}^N Q_{jn}(β) (1/(1 + |ζ_j|^2) - Σ_{m=0}^M c_m Q_{jm}(β)) = 0

in the matrix form

A(β)c = b(β), (15)

where

A_{nm}(β) = Σ_{j=1}^N Q_{jn}(β) Q_{jm}(β)

and

b_n(β) = Σ_{j=1}^N 2 Q_{jn}(β) / (1 + |ζ_j|^2).

Then, the solution of (15) reads

c(β) = A^{-1}(β) b(β).

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

J_incoh(β) = 1/N_c Σ_{j: |ζ_j| < 1 - 1/√N} |ξ_j - 1/(1 + |ζ_j|^2) Σ_{m=0}^M c_m(β) Q_{jm}(β)|^2,

and look for its global minimum β_min, 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 |ζ_j| < 1 - 1/√N in our definition of J_incoh(β).) To find the global minimum of J_incoh(β), 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 J'_incoh(β) = 0 using Newton’s method. The obtained global minimizer β_min is interpreted as an approximate value of the phase lag β in Eq. (5), see Fig. 4b. Respectively, formula (14) with c_m = c_m(β_min) 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 Q_{jm}(β) and Q_{jm}(β). All their elements can be computed in O(N^2 M) operations, while the entire matrix A(β) in Eq. (15) can be computed in O(NM^2) operations. Moreover, solving the linear system (15) requires an additional O(M^3) operations. Together, this means that the minimization of the functionals J({c_m}, β) and J_incoh(β) requires O(N^2 M + NM^2 + M^3) operations. Therefore, if M ≪ N, the resulting computational cost is much less than the standard number of operations O(N^3) 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).

For a fixed 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 c_m in formula (14). Note that for the top-hat coupling function, these coefficients can be calculated analytically

c_0 = 1/(2π), c_m = sin(πmσ)/(π^2 mσ) for m = 1, 2, ..., (16)

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.

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 c_m 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 c_1, we can uniquely determine the ranges and decay rates of the corresponding coupling functions, which confirm the reliability and efficiency of the proposed approach.

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 (x_j, Ω_j, ζ_j), j = 1, …, N, only a sufficiently large subset of them is used. Roughly speaking, from N points x_j we can randomly select a smaller subset {x_j : 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({c_m}, β) and J_incoh(β). 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

Z(t) = 1/#{S} Σ_{j∈S} e^{iθ_j(t)}.

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.

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

ζ_j = 1/T ∫_0^T e^{iθ_j(t)} / |Z(t)| dt,

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

Ω_j = lim_{T→∞} 1/T Σ_{k=2}^{N_T} (θ_j(t_k) - θ_j(t_{k-1})).

Here, we need to satisfy the Nyquist criterion that the minimum of the intervals t_k - t_{k-1} is less than the half-period of the corresponding oscillations. But the intervals t_k - t_{k-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 exploit 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 states42 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

Z_2(t) = 1/N Σ_{k=1}^N e^{2iθ_k(t)}

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

ζ_j = ⟨e^{iθ_j(t)} / √|Z_2(t)|⟩,

where the complex square root is calculated in such a way that it varies as a continuous function of time. The value of ζ_j 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

ζ_j = ζ_j / |Z_2|, where Z_2 = 1/N Σ_{j=1}^N ζ_j.

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 models of limit cycle oscillators. The physical phase is defined as a variable that increases uniformly from 0 to 2π during one cycle of its oscillations10,11. It is this phase that appears in the formulation of Kuramoto-type models obtained using the phase reduction method10,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 transformation43. This means that additional effort is usually required to reconstruct the function Φ(θ) along with other system parameters. Fortunately, for stationary coherence-incoherence patterns, such effort 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, the effective frequencies are calculated by

Ω_j = 1/T ∫_0^T dϕ_j(t)/dt dt.

Recalling that the coherent region of the chimera state corresponds to a frequency Ω_c, as in Fig. 2, we can easily verify that Ω_j = Ω_c. 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 Z(t) ≈ Z_0 e^{iΩ_c t}.

Next, we consider the modified local order parameters

ζ_j = ⟨e^{iϕ_j(t)} e^{-iΩ_c t}⟩.

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 ζ_j, 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

ζ_j = ζ_j / Z_b, where Z_b = 1/N Σ_{j=1}^N ζ_j.

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 Φ(θ).

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):

dθ_j / dt = ω - Im(W_j e^{-iθ_j} e^{iα}), (17)

where

W_j = 2π/N Σ_{k=1}^N G(x_j - x_k) e^{iθ_k} (18)

and W_j denotes the complex-conjugate of W_j. 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

∂ρ/∂t + ∂/∂θ [(ω - ∫_{-π}^π ∫_{-π}^π G(x - x_0) sin(θ - θ_0 + α) ρ(θ_0, x_0, t) dθ_0 dx_0) ρ] = 0.

It is known27,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 manifold68,69 consisting of functions of the form

ρ(θ, x, t) = 1/(2π) [1 + Σ_{n=1}^∞ (z_n(x, t) e^{inθ} + z_n(x, t) e^{-inθ})],

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

dz/dt = iωz + 1/2 e^{-iα} Gz - 1/2 e^{iα} z^2 Gz (19)

with the integral operator

(Gz)(x, t) = ∫_{-π}^π G(x - x_0) z(x_0, t) dx_0,

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

z(x, t) = a(x) e^{iΩt} (20)

Inserting this ansatz into Eq. (19) and denoting

w(x) = ∫_{-π}^π G(x - x_0) a(x_0) dx_0, (21)

we obtain

e^{-iβ} w(x) a^2(x) - 2a(x) + e^{iβ} w(x) = 0. (22)

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

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

∫_{-π}^π ∫_{-π}^π ρ(θ, x, t) dθ dx = 2π, ∫_{-π}^π e^{ikθ} ρ(θ, x, t) dθ = a_k(x) e^{ikΩt} for k ∈ N. (23)

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

Z(t) = 1/(2π) ∫_{-π}^π ∫_{-π}^π e^{iθ} ρ(θ, x, t) dθ dx = Z_0 e^{iΩt},

where

Z_0 = 1/(2π) ∫_{-π}^π a(x) dx.

Similarly, using the ergodicity property, we obtain

ζ_j = 1/|Z_0| ∫_{-π}^π e^{iθ} e^{-iΩt} ρ(θ, x_j, t) dθ = a(x_j) / |Z_0|, (24)

and therefore

ζ_j = |Z_0| as N → ∞. (25)

Let us denote a_j = a(x_j) and w_j = w(x_j), then using the trapezoidal rule we write an approximate version of the definition (21)

w_j = 1/(ω - Ω_j) Σ_{k=1}^N G(x_j - x_k) a_k (x_{k+1} - x_{k-1})/2, (26)

where due to the periodicity of the variable x we assume x_{N+1} - x_N = 2π + x_1 - x_N and x_1 - x_{N-1} = 2π + x_1 - x_{N-1}. Multiplying this by e^{iβ} / |Z_0| and using (24), we obtain

ξ_j = w_j e^{iβ} / |Z_0|, (27)

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

ξ_j ζ_j^2 - 2ζ_j + ξ_j = 0. (28)

Proposition. Suppose that ξ_j ∈ C and ζ_j ∈ C satisfy equation (28), and moreover |ζ_j| ≤ 1 for all j ∈ [−π, π]. then

ξ_j = 2ζ_j / (1 + |ζ_j|^2) for |ζ_j| ≠ 1 (29)

and

Re(ξ_j ζ_j) = 1 for |ζ_j| = 1. (30)

Moreover, for all values of |ζ_j| we have

Re(ξ_j / ζ_j) = 1 / (1 + |ζ_j|^2) and Re(ξ_j ζ_j) = 2|ζ_j|^2 / (1 + |ζ_j|^2).

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

ξ_j ζ_j^2 - 2ζ_j + ξ_j = 0,

or equivalently ξ_j = 2ζ_j - ξ_j ζ_j^2. Inserting this into Eq. (28), we obtain

2|ζ_j|^2 ζ_j - ξ_j |ζ_j|^4 - 2ζ_j + ξ_j = 0,

or

ξ_j (1 - |ζ_j|^4) = 2ζ_j (1 - |ζ_j|^2).

If |ζ_j| ≠ 1, this yields (29). On the other hand, if |ζ_j| = 1, then ζ_j = 1/ζ_j, and therefore dividing (28) by ζ_j, we obtain

ξ_j ζ_j - 2 + ξ_j ζ_j = 0

what is equivalent to (30). ■

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

W_j = 2π/N Σ_{k=1}^N G(x_j - x_k) a_k e^{iΩt} = (ω - Ω_j) w_j e^{iΩt} = (ω - Ω_j) ξ_j e^{-iβ} |Z_0| e^{iΩt}.

Inserting this into (17), we obtain

dθ_j / dt = ω - (ω - Ω_j) Im(i ξ_j ζ_j e^{-iΩt} e^{iθ_j}).

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

Ω_j = ⟨dθ_j / dt⟩ = ω - (ω - Ω_j) Im(i ξ_j ζ_j) = ω - (ω - Ω_j) Re(ξ_j ζ_j) = ω - (ω - Ω_j) 2|ζ_j|^2 / (1 + |ζ_j|^2), (31)

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

dϕ_j / dt = dθ_j / dt = Ω_j. (32)

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

ζ_j = ⟨e^{iϕ_j(t)} e^{-iΩ_c t}⟩.

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

e^{iΦ(θ)} = Σ_{k=-∞}^∞ φ_k e^{ikθ} with some φ_k ∈ C,

and therefore

ζ_j = ⟨e^{iΦ(θ_j(t))} e^{-iΩ_c t}⟩ = Σ_{k=-∞}^∞ φ_k ⟨e^{ikθ_j(t)} e^{-iΩ_c t}⟩.

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

⟨e^{ikθ_j(t)} e^{-iΩ_c t}⟩ = ∫_{-π}^π e^{ikθ} e^{-iΩ_c t} ρ(θ, x_j, t) dθ = a_k(x_j) e^{i(k-1)Ω_c t} = a(x_j) δ_{k1} = |Z_0| ζ_j δ_{k1} for k ∈ N,

and similarly

⟨e^{ikθ_j(t)} e^{-iΩ_c t}⟩ = 0 for k = 0, -1, -2, ...

This means

ζ_j = φ_1 |Z_0| ζ_j.

Defining

Z_b = 1/N Σ_{j=1}^N ζ_j (33)

and using (25), we obtain Z_b = φ_1 Z_0 and

ζ_j / Z_b = ζ_j / |Z_0|.

Given that for each coherent oscillator we have |ζ_j| = 1, and for each incoherent oscillator |ζ_j| < 1, we write |φ_1| = max_j |ζ_j|. Thus, we obtain the formulas

ζ_j = ζ_j / Z_b * max_k |ζ_k|, (34)

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

dθ_j / dt = ω_j - q_j/N Σ_{k=1}^N q_k sin(θ_j - θ_k + α), (35)

where the natural frequencies ω_j and coupling strengths q_j > 0 are drawn randomly and independently from distributions H_ω(ω) and H_q(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 density55,70,71

ρ(θ, ω, q, t) = H_ω(ω) H_q(q) / (2π) [1 + Σ_{n=1}^∞ (a_n(ω, q, t) e^{inθ} + a_n(ω, q, t) e^{-inθ})],

where

a(ω, q, t) = h((ω - Ω) / (pq)) e^{iΩt}

and

h(s) = (1 - √1 - s^2) / s for |s| ≤ 1, (1 - i√s^2 - 1) / s for |s| > 1, (36)

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

p = i e^{-iα} ∫_{-∞}^∞ ∫_0^∞ H_ω(ω) H_q(q) h((ω - Ω) / (pq)) q dq. (37)

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

dθ_j / dt = ω_j - Im(W e^{iθ_j} e^{iα}),

where

W(t) = 1/N Σ_{k=1}^N q_k e^{iθ_k(t)},

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

ζ_k = (ω_k - Ω) / (p q_k) h((ω_k - Ω) / (p q_k)), (38)

Ω_k = Ω + p q_k Q((ω_k - Ω) / (p q_k)), (39)

where K = p e^{iα} / |Z_0|, where |Z_0| is replaced by the time-averaged magnitude of the global order parameter Z(t) defined by (6).

Z_0 = ∫_{-∞}^∞ ∫_0^∞ H_ω(ω) H_q(q) h((ω - Ω) / (pq)) dq dω (40)

and

Q(s) = s (1 - √1 - s^2) / s for |s| ≤ 1, s (1 - i√s^2 - 1) / s for |s| > 1.

According to their derivation, relations (38) and (39) are exact for N → ∞ and independent of the distributions H_ω(ω) and H_q(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 model72

dθ_j / dt = ω_j - K/N Σ_{k=1}^N sin(θ_j - θ_k + α), (41)

which describes the dynamics of all-to-all coupled phase oscillators with different frequencies ω_j. In this case H_q(q) = δ(q - K), and therefore the self-consistency equation (37) and formula (40) can be written in the form

p = i K e^{-iα} Z_0, (42)

where

p = p / K and Z_0 = ∫_{-∞}^∞ H_ω(ω) h((ω - Ω) / p) dω.

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

ζ_k = i e^{-iα} h((ω_k - Ω) / p), (43)

Ω_k = Ω + p Q((ω_k - Ω) / p). (44)

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 - 2sh + 1 = 0 and relations (43) hold, the parameter α can be found from an explicitly solvable minimization problem

-α = argmin_{β ∈ [0, π]} Σ_{k=1}^N |ζ_k^2 e^{-2iβ} + 1 - 2ζ_k e^{-iβ} s_k|^2,

and moreover

s_k = (ω_k - Ω) / p = Re(ζ_k^2 e^{2iα} + 1) / (2ζ_k e^{iα}).

Once the ratios s_k are known, the parameters Ω and p can be determined by linear regression from the set of equations (44). Then, ω_k = Ω + p s_k. And finally, from (42) we get the coupling strength K.

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 equations9

dθ_j / dt = ω_0 - K/N Σ_{k=1}^N a_{jk} sin(θ_j - θ_k + α), (45)

where {a_{jk}} denotes an adjacency matrix such that a_{jk} = 1 if there is a link between the jth and kth oscillators, and a_{jk} = 0 otherwise. The number d_j = Σ_{k=1}^N a_{jk} is called the degree of the oscillator j, and d_mean = 1/N Σ_{j=1}^N d_j denotes the mean degree. For densely connected random networks, their dynamics can be described using the annealed approximation9,73

dθ_j / dt = ω_0 - K/(N d_mean) Σ_{k=1}^N d_j d_k sin(θ_j - θ_k + α),

or equivalently

dθ_j / dt = ω_0 - 1/N Σ_{k=1}^N q_j q_k sin(θ_j - θ_k + α), (46)

where q_j = d_j K / d_mean.

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

ζ_k = Z_0 h(s_k), (47)

Ω_k = Ω + (ω_0 - Ω) s_k Q(s_k), (48)

where s_k = (ω_0 - Ω) / (p q_k).

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

p = i e^{-iα} W_0, (49)

where W_0 = ∫_0^∞ H_q(q) h((ω_0 - Ω) / (pq)) q dq. (50)

Next, we will demonstrate that the parameters ω_0, α and q_j in Eq. (46) can be found from the observed values of ζ_k and Ω_k. First, we note that the definitions (36) and (40) imply that Im(Z_0) ≤ 0 and -π ≤ arg Z_0 ≤ 0. Recalling that h(s) satisfies the quadratic equation h^2 - 2sh + 1 = 0 and relations (47), we write

arg Z_0 = ϕ_0 := argmin_{ϕ ∈ [−π, 0]} Σ_{k=1}^N |ζ_k^2 e^{2iϕ} + 1 - 2ζ_k e^{iϕ} s_k|^2,

and s_k = Re(ζ_k^2 e^{2iϕ_0} + 1) / (2ζ_k e^{iϕ_0}).

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

1/N Σ_{k=1}^N p q_k h(s_k) = 1/N Σ_{k=1}^N p q_k h(s_k) = p W_0.

Substituting this into the self-consistency equation (49), we finally get p = √S, α = arg(iS), and q_k = (ω_0 - Ω) / (p s_k).

Note that using the above approach, we can only find the values of q_k = K d_k / d_mean, but we cannot separate them into the coupling strength K and the degree d_k. 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 H_q(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 Z_0 or W_0, 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.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Report on \"Inverse problems for dynamic patterns in coupled oscillator networks: When larger networks are simpler” by\nOmel’chenko\nIn this very well written and original contribution the author employs sophisticated (and, to this referee, rather surprising)\nexplicit analytical expressions for certain thermodynamic variables to determine unknown parameters of the dynamical\nsystem on a network of many thousand oscillators. As an example to illustrate their method the author chose the challenging\nproblem of chimera states which have attracted much attention in the last 10 years, from physicists, mathematicians and\nbiologists. The author corroborates their theory with excellent numerical results, which are clearly presented, illustrate the\nefficacy and accuracy of the proposed method and are reproducible given the matlab code, provided by the author.\nThe methods and general strategy developped here will be applicable to a much wider class of systems (the author actually\nchose a very challenging one), including power grids. I therefore believe that this work will be of interest to a large part of the\nscientific community.\nAll calculations are explicitly given and the numerical code\nThe author's method to infer unknown parameters from a time series has several advantages over other existing\napproaches: robustness to (unbiased) noise, non-reliance on equidistantly sampled data, no reliance on knowledge of\ntemporal derivatives (and instead only using smoothing averaging), no requirement on full observations, and no inversion of\nof NxN matrices (where N is the number of dynamic variables). This makes the newly proposed method stand out against\nprevious methods and makes it a very powerful and widely applicable new methodological framework.\nI only have a few comments the author may wish you take on board.\nOn page 5, the author restricts to symmetric functions G(x). Could the author comment on this choice and also on possible\nways to tackle anti-symmetric G or G which do not have any symmetry (are there any chimeras in these situations\npossible?)?\nI was wondering about the number of Fourier modes required to infer the parameters. It is nicely shown that 10 Fourier\nmodes are sufficient to obtain a very good accuracy for the parameters, despite the coupling function G(x) not being\naccurately represented (as seen in Fig 4, right figure). Does the accuracy grow monotonically with the number of Fourier\nmodes or will the Gibb’s phenomenon lead to a degradation?\nA brief comment on the computational cost of minimizing the function J would be helpful to know.\nThe author lists a few works on other approaches to infer unknown parameters, but missed to mention the recent work\nSmith & Gottwald, Proc Roy Soc A (2025), 481 (2025) 20240813, entitled \"Data assimilation for networks of coupled oscillators:\nInferring unknown model parameters from partial observations”. This method, employing the ensemble Kalman filter and\nBayesian inference allows for noisy partial observations (but also smaller networks for which the thermodynamic limit would\nnot be valid), but requires a constant sampling time and knowledge of the coupling function.\nTypos:\n* comma after (30 missing\n* page 2: 3rd line in Results: of a chimera state\n* several commata are missing in formulae which are followed by the word “where”", "reviewer_2": "The author proposed an inverse approach to infer on system parameters and more importantly to reconstruct the coupling\nfunctions using the chimera pattern data obtained from the Kuramoto-Battogtokh model with nonlocal coupling, based on\nobservation of the average frequency and the local order parameter. In particular, the author used three different coupling\nschemes to validate their method. While this approach might be a promising step towards the inverse problem, in general,\nfor the future, with its analytical foundation and numerical confirmation with different coupling functions, I have major concern\nabout the limitations of this method.\n1. While the conclusion suggests that the method may have broad applicability, the abstract and introduction do not\nadequately convey the underlying motivation and physical framework necessary for inverse analysis. A more explicit\ndiscussion on the advantages and disadvantages of the proposed approach would be valuable. For example, it is not clear\nhow the reconstruction would perform for systems with more complex connectivity structures—such as small-world networks,\nand even globally coupled network—rather than the nonlocal topology considered here. It would also be important to clarify\nwhether this method can handle distributed phase lag parameters (α) or heterogeneity in intrinsic parameters such as the\nnatural frequencies. In particular, could partial synchronization arising in classical Kuramoto oscillators with complex\nconnectivity and Lorentzian-distributed frequencies still allow for recovery of the coupling function and frequency set? If not,\nthese limitations should be explicitly acknowledged in the text.\n2. The current analysis is limited to phase oscillator networks that exhibit a particular chimera pattern. A broader discussion\nof how the inverse method might be adapted to other dynamical regimes—such as clustering states, breathing chimeras, or\nturbulent incoherent states—would strengthen the claim of generality.\n3. The study focuses on a chimera configuration with two coherent/incoherent domains. It would be insightful to explore the\npotential outcomes if multi-headed chimera states were used instead. Would the SER-based reconstruction remain robust?\nSuch an extension could provide valuable insights into the generality and scalability of the proposed approach.\n4. It seems, the current work assumes a coupling strength fixed at unity. It would be useful to discuss whether the algorithm\ncould be extended to cases with variable coupling strengths. If so, would the method still reliably infer other parameters\nunder these conditions, or would additional observables/constraints be needed? Addressing this would clarify the flexibility of the\nmethod for more realistic scenarios.\n5. The reconstruction of coupling functions in nonlinear dynamical systems is a well-established topic. In particular, methods\nfor identifying coefficients in polynomial expansions such as the SINDy framework (Proc. Natl. Acad. Sci. U.S.A. 113, 3932\n(2016)) and other related approaches (Proc. Natl. Acad. Sci. U.S.A. 115, 9300 (2018); Phys. Rev. Lett. 130, 117401 (2023);\nChaos 29, 103116 (2019)) have been extensively studied. The author should clearly explain the specific novelty of their\nmethod compared to these existing techniques. For example, what distinct advantages does the use of statistical equilibrium\nrelations provide such as scalability to very large networks, robustness to measurement noise, or applicability with\npartial/uneven data that go beyond the capabilities of these prior methods?", "reviewer_3": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "Response to Reviewers’ comments on\n“Inverse problems for dynamic patterns\nin coupled oscillator networks:\nWhen larger networks are simpler”\nby Oleh Omel’chenko\nWe thank the Reviewers for their careful reading of our manuscript and for their valuable\ncomments on our results. In the revised version of the manuscript, we have attempted to take\nallofthesecommentsintoaccountbyaddingnewmaterialandexpandingourexplanations. We\nare pleased that both Reviewers recognized the potential of our method and raised questions\nabout how it would work in other situations. We hope that our responses provide a general\nanswer to these questions, and more detail can be obtained from further research. In any case,\ntheaimofourworkistodescribeanewanalyticalmethodbyapplyingittoseveralparadigmatic\nexamples, rather than an exhaustive description of all the results that can be obtained using\nit. We hope that the responses below will be found satisfactory by the Reviewers, and that the\nresubmitted manuscript will be accepted for publication in Nature Communications.\nResponse to Reviewer 1:\nIn this very well written and original contribution the author employs sophisticated (and, to\nthis referee, rather surprising) explicit analytical expressions for certain thermodynamic vari-\nables to determine unknown parameters of the dynamical system on a network of many thousand\noscillators. As an example to illustrate their method the author chose the challenging problem\nof chimera states which have attracted much attention in the last 10 years, from physicists,\nmathematicians and biologists. The author corroborates their theory with excellent numerical\nresults, which are clearly presented, illustrate the efficacy and accuracy of the proposed method\nand are reproducible given the matlab code, provided by the author.\nThe methods and general strategy developped here will be applicable to a much wider class of\nsystems (the author actually chose a very challenging one), including power grids. I therefore\nbelieve that this work will be of interest to a large part of the scientific community.\nAll calculations are explicitly given and the numerical code\nThe author’s method to infer unknown parameters from a time series has several advantages\nover other existing approaches: robustness to (unbiased) noise, non-reliance on equidistantly\nsampled data, no reliance on knowledge of temporal derivatives (and instead only using smooth-\ning averaging), no requirement on full observations, and no inversion of of N × N matrices\n(where N is the number of dynamic variables). This makes the newly proposed method stand out\nagainst previous methods and makes it a very powerful and widely applicable new methodological\nframework. I only have a few comments the author may wish you take on board.\nResponse: We thank the Reviewer for this positive assessment of our manuscript.\nOn page 5, the author restricts to symmetric functions G(x). Could the author comment on\nthis choice and also on possible ways to tackle anti-symmetric G or G which do not have any\nsymmetry (are there any chimeras in these situations possible?)?\nResponse: Our study focuses on stationary coherence-incoherence patterns, i.e., patterns\nthat do not change their position on the network and do not change their shape at the macro-\nscopic level. (We have now added a better explanation of this fact to the last two paragraphs\nof the section “Statistical equilibrium relations (SER)” on page 4.) On the other hand, it is\nwell-known that asymmetric coupling topologies with G(−x) ̸= G(x) typically lead to patterns\nthat drift across the network. Such are, for example, travelling chimera states described in [1].\nIn this case, uniform drift causes the local order parameters of different oscillators (as well as\ntheir effective frequencies) to become the same after averaging over a long time interval. This\nmeansthatallinformationaboutthenetworkheterogeneitiesislostandcannotberecovered. It\nis likely that a modified version of our method can be developed that will allow the calculation\nof local order parameters and effective frequencies in a suitably chosen corotating frame. Then\nthe coupling function G(x) should be sought in the form of a general Fourier series with cosine\nand sine harmonics cos(mx) and sin(mx). However, such a consideration is beyond the scope\nof the presented study, so we leave it for future research.\nI was wondering about the number of Fourier modes required to infer the parameters. It\nis nicely shown that 10 Fourier modes are sufficient to obtain a very good accuracy for the\nparameters, despite the coupling function G(x) not being accurately represented (as seen in\nFig 4, right figure). Does the accuracy grow monotonically with the number of Fourier modes\nor will the Gibb’s phenomenon lead to a degradation?\nResponse: For a given set of local order parameters ζ and effective frequencies Ω ,\nthe\nk k\nL2-accuracy of the reconstruction of the coupling function G(x) increases almost monotonically\nwith the number of Fourier modes M, but only up to some threshold value M , while for\nopt\nM > M the accuracy starts to deteriorate, see Fig. 1 below. Remarkably, the value M is\nopt opt\nnot related to the Gibbs phenomenon, but to the finite accuracy of the statistical equilibrium\nrelations. Since these relations are exact for an infinite number of oscillators N and an infinite\naveraging time T, we know that M → ∞ as N,T → ∞. But a more specific functional\nopt\ndependence of M on N and T is difficult, if not impossible, to obtain.\nopt\n[Figure 1 omitted]\nFigure 1: L2-accuracy of the reconstructed coupling function G (x) relative to the original\nreconstr\ncoupling function G(x). The chimera state for the parameters from Fig. 2 in the manuscript\nwas numerically simulated for system sizes N = 1024, 2048, and 4096, and the corresponding\neffective frequencies Ω and local order parameters ζ were calculated as time-averaged over\nk k\nT = 2000. The obtained data were used in our reconstruction algorithm with different M. In\neach case, the accuracy initially increases with M, but for the number of Fourier harmonics M\ngreater than some M starts to deteriorate. The black curve shows the maximum possible\nopt\naccuracy provided by a truncated Fourier series with all spatial modes cos(mx) where m ≤ M.\nA brief comment on the computational cost of minimizing the function J would be helpful\nto know.\nResponse: We followed the Reviewer’s recommendation and added a paragraph discussing\nthe computational cost of our algorithm at the end of the section “Parameter reconstruction\nalgorithm”, see page 6.\nThe author lists a few works on other approaches to infer unknown parameters, but missed to\nmention the recent work Smith & Gottwald, Proc Roy Soc A (2025), 481 (2025) 20240813, enti-\ntled ”Data assimilation for networks of coupledoscillators: Inferring unknown model parameters\nfrom partial observations”. This method, employing the ensemble Kalman filter and\nBayesian inference allows for noisy partial observations (but also smaller networks for which thethermo-\ndynamic limit would not be valid), but requires a constant sampling time and knowledge of the\ncoupling function.\nResponse: We are very grateful that this recently published work was brought to our\nattention. We have included it in the bibliography and cited it in the section “Discussion”.\nTypos:\n* comma after (30 missing\n* page 2: 3rd line in Results: of a chimera state\n* several commata are missing in formulae which are followed by the word “where”\nResponse: We thank the Reviewer for spotting these typos. All they were corrected.\nResponse to Reviewer 2:\nThe author proposed an inverse approach to infer on system parameters and more impor-\ntantly to reconstruct the coupling functions using the chimera pattern data obtained from the\nKuramoto-Battogtokh model with nonlocal coupling, based on observation of the average fre-\nquency and the local order parameter. In particular, the author used three different coupling\nschemes to validate their method. While this approach might be a promising step towards the\ninverse problem, in general, for the future, with its analytical foundation and numerical con-\nfirmation with different coupling functions, I have major concern about the limitations of this\nmethod.\n1. While the conclusion suggests that the method may have broad applicability, the abstract\nand introduction do not adequately convey the underlying motivation and physical framework\nnecessary for inverse analysis. A more explicit discussion on the advantages and disadvantages\nof the proposed approach would be valuable. For example, it is not clear how the reconstruction\nwould perform for systems with more complex connectivity structures — such as small-world\nnetworks, and even globally coupled network — rather than the nonlocal topology considered\nhere. It would also be important to clarify whether this method can handle distributed phase\nlag parameters (α) or heterogeneity inintrinsic parameters such as the natural frequencies. In\nparticular, could partial synchronization arising in classical Kuramoto oscillators with complex\nconnectivity and Lorentzian-distributed frequencies still allow for recovery of the coupling func-\ntion and frequency set? If not, these limitations should be explicitly acknowledged in the text.\nResponse: We are very grateful to the Reviewer for this extremely important comment.\nAlthough our initial intention was to describe the concept of statistical equilibrium relations in\nthe context of chimera states, we now see that other applications should be mentioned in the\nmanuscriptforabetterunderstandingofourmethod. Accordingly, wehaveaddedanewsection\n“Otherapplications”afterthesection“Methods”. Inthenewsection, weexplainhowstatistical\nequilibrium relations can be obtained for the classical Kuramoto-Sakaguchi model and for the\nannealed approximation of the Kuramoto model on a random network. For the first model,\nwe also show how statistical equilibrium relations can be used to find all oscillator frequencies,\ncoupling strength, and phase lag parameter from observed values of local order parameters and\neffective frequencies. (Importantly, this can be done for arbitrary frequency distributions, not\njust Lorentzians!) For the second model, representing a Kuramoto-type random network, we\nshow that statistical equilibrium relations can be used to recover the relative node degrees and\ntheir distribution, which can potentially help to identify the type of the network. Finally, we\nexplain that in both of the above cases, the parameter reconstruction algorithm can be applied\nto partially collected data, provided that the observed oscillators are uniformly distributed\nacross the network.\nTo reflect the above changes to the text, we have updated the abstract and introduction\naccordingly.\nAs for the case of distributed phase lag parameters, it also has a thermodynamic limit\ndescription and corresponding statistical equilibrium relations. But these relations contain\nother functional dependencies, so their analysis in the context of parameter reconstruction will\nbe qualitatively different from the two examples that we have already added to the manuscript\nin the section “Other applications”. Honestly, we currently do not have a solution to the\nparameter reconstruction problem in this case and prefer to work on this in a separate paper.\nAgain, we cannot apply our method to all interesting problems at once.\n2. The current analysis is limited to phase oscillator networks that exhibit a particular\nchimera pattern. A broader discussion of how the inverse method might be adapted to other\ndynamical regimes — such as clustering states, breathing chimeras, or turbulent incoherent\nstates — would strengthen the claim of generality.\nResponse: Theapproachdevelopedinourworkfocusesonstationarycoherence-incoherence\npatterns and we clearly state this in the text. Clustering states, where the oscillator network\nis divided into several groups with identical dynamical behaviour in each of these groups, are\nnot suitable for parameter reconstruction. And this applies to any parameter reconstruction\nmethod, not just ours. Imagine that you observe a two-cluster state in a network of N ≫ 2\noscillators, then the details about the heterogeneities of the oscillators and the coupling topol-\nogy in each of the clusters are no longer present in the observed data and therefore cannot be\nrecovered! Similarly, turbulent incoherent states are also a poor source of input data for our\nmethod, as we only use time-averaged quantities, and for incoherent states they typically take\nthe same value across the network due to the ergodicity property. As for breathing chimeras,\nthey are non-stationary coherence-incoherence patterns, so the current version of our param-\neter reconstruction algorithm cannot be applied to them. But, using the analytical approach\npresented in [2] one can hope to obtain statistical equilibrium relations for this type of patterns\nas well. However, such a consideration is beyond the scope of the presented study, therefore\nwe leave it for future research. To emphasize this, we have added new text to the sections\n“Statistical equilibrium relations (SER)”, see page 4, and “Discussion”, see page 10.\n3. The study focuses on a chimera configuration with two coherent/incoherent domains. It\nwould be insightful to explore the potential outcomes if multi-headed chimera states were used\ninstead. Would the SER-based reconstruction remain robust? Such an extension could provide\nvaluable insights into the generality and scalability of the proposed approach.\nResponse: In general, the multi-headed chimera states found in the Kuramoto-Battogtokh\nmodel can be analyzed in the same way as the chimera states with two coherent/incoherent\ndomains. The SER-based parameter reconstruction algorithm remains unchanged. The only\nproblem that may arise is that multi-headed chimera states usually have antiphase adjacent\ncoherent regions, and therefore their global order parameter is almost zero. This circumstance\nmakes the definition of the local order parameters ζ given in the Introduction unsuitable. In\nj\nthe new section “Chimera states with zero global order parameter” (see pages 8–9) we explain\nhow this problem can be solved using a Daido order parameter.\n4. It seems, the current work assumes a coupling strength fixed at unity. It would be useful\nto discuss whether the algorithm could be extended to cases with variable coupling strengths.\nIf so, would the method still reliably infer other parameters under these conditions, or would\nadditional observables/constraints be needed? Addressing this would clarify the flexibility of the\nmethod for more realistic scenarios.\nResponse: There is nothing special about the coupling function G(x) being normalized to\nunity in all our examples. It is simply a tribute to the tradition established in the literature\non chimera states. If you examine the rationale for our parameter reconstruction algorithm,\nyou will find that the normalization condition is never used. Therefore, the proposed algorithm\nworks equally well for coupling functions with other normalizations. We added an explanation\nof this fact at the end of the section “Examples” on page 7.\n5. The reconstruction of coupling functions in nonlinear dynamical systems is a well-\nestablished topic. In particular, methods for identifying coefficients in polynomial expansions\nsuch as the SINDy framework (Proc. Natl. Acad. Sci.U.S.A. 113, 3932 (2016)) and other\nrelated approaches (Proc. Natl. Acad. Sci. U.S.A. 115, 9300 (2018); Phys. Rev. Lett. 130,\n117401 (2023); Chaos 29, 103116 (2019)) have been extensively studied. The author should\nclearly explain the specific novelty of their method compared to these existing techniques. For\nexample, what distinct advantages does the use of statistical equilibrium relations provide such\nas scalability to very large networks, robustness to measurement noise, or applicability with\npartial/uneven data that go beyond the capabilities of these prior methods?\nResponse: The SINDy framework is an extremely general method designed to reconstruct\nan arbitrary n-dimensional dynamical system of the form x˙ = f(x) from noisy measurement\ndata [3]. But its generality can also be a source of problems. According to its authors, this\n“method is ill-suited for a large state dimension n, because of the factorial growth of Θ [a\nlibrary vector consisting of candidate nonlinear functions of x] in n and because each of the n\nrow equations in Eq. 4 [a system of linear equations that determines the coefficients at different\nterms of library functions] requires a separate optimization” [3, p. 3933]. Even if we assume\na certain ansatz for f(x) (which is actually done in the present manuscript), using the SINDy\napproach we will obtain an (nm)×(nm) system of linear equations, where m is the number of\nsampling times. Finding a solution to this system usually requires O((nm)3) operations, while\nour SER-based algorithm requires only O(n2) operations, as explained in the new paragraph\nof the section “Parameter reconstruction algorithm” on pages 6–7.\nThe authors of [3] also describe a situation where the curse of dimensionality can be over-\ncome relatively easy: “Fortunately, many high-dimensional systems of interest evolve on a\nlow-dimensional manifold or attractor that is well-approximated using a low-rank basis Ψ” [3,\np. 3933]. Regarding this remark, we note that most partially synchronized states on coupled\noscillator networks, including chimera states, are hyperchaotic trajectories with Kaplan-Yorke\ndimension proportional to the system size n [4]. The low-dimensional nature of the dynam-\nics of such states becomes apparent only in the thermodynamic limit and only when using a\nprobability density representation. This is another feature that is taken into account in our\nSER-based algorithm, but cannot be simply implemented in the SINDy framework.\nThe fact that we do not need to calculate time derivatives and can use partial observations\nis also an advantage of the SER-based algorithm compared to the SINDy approach. But our\nmethod is more focused on a specific (albeit quite broad) class of models and must be tuned to\neach model individually.\nA more specific advantage of our method is explained in the new section “Calculations\nwith protophases” on pages 9–10. We show that the effective frequencies and local order\nparameters required for our parameter reconstruction algorithm can be easily calculated even\nin the case where, instead of the physical phases θ (t), we observe their transformed values\nj\nϕ (t) = Φ(θ (t)) with an unknown transformation function Φ(θ). This situation is typical in\nj j\nreal-world applications, where the physical phase is not measured directly, but is reconstructed\nfrom a protophase obtained from an oscillatory signal [5].\nResponse to Reviewer 3:\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This\nis part of the Nature Communications initiative to facilitate training in peer review and to\nprovide appropriate recognition for Early Career Researchers who co-review manuscripts.\nResponse: All comments from the reports of Reviewer 1 and Reviewer 2 were very useful\nto us and helped to improve the presentation and clarity of our results. Therefore, we would\nlike to thank Reviewer 3 for his/her contribution to one of the above reports.\nResponse to Reviewer 4:\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This\nis part of the Nature Communications initiative to facilitate training in peer review and to\nprovide appropriate recognition for Early Career Researchers who co-review manuscripts.\nResponse: All comments from the reports of Reviewer 1 and Reviewer 2 were very useful\nto us and helped to improve the presentation and clarity of our results. Therefore, we would\nlike to thank Reviewer 4 for his/her contribution to one of the above reports."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The author has addressed all my concerns and I am very happy to recommend publication. Great work!\nJust one suggestion the author may take on board if he sees fit: Fig 1 from the response letter could be included in a supplementary file together with a discussion. But this should be left to the author's decision on how this fits.", "reviewer_2": "The author responded to our questions, hence recommend acceptance for publication.", "reviewer_3": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "Response to Reviewers’ comments on\n“Inverse problems for dynamic patterns\nin coupled oscillator networks:\nWhen larger networks are simpler”\nby Oleh Omel’chenko\nWe thank the Reviewers for all their comments (from this and previous reports), which\ngreatly helped us to improve the quality of our work. We hope that with the revisions made,\nthe manuscript will finally be accepted for publication in Nature Communications.\nResponse to Reviewer 1:\nThe author has addressed all my concerns and I am very happy to recommend publication.\nGreat work!\nJust one suggestion the author may take on board if he sees fit: Fig. 1 from the response\nletter could be included in asupplementary file together with a discussion. But this should be\nleft to the author’s decision on how this fits.\nResponse: We thank the Reviewer for this positive assessment of our work. We followed\nhis/her suggestion and added Fig. 1 from the previous response letter along with its discussion\nto the Supplementary Information file (see Supplementary Note 1). Accordingly, we added two\nnew sentences referring to this note in the main text after Eq. (14).\nResponse to Reviewer 2:\nThe author responded to our questions, hence recommend acceptance for publication.\nResponse: We are glad we were able to answer all of the Reviewer’s questions, which\nprompted us to describe additional aspects of our method.\nResponse to Reviewer 3:\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports.\nThis is part of the NatureCommunications initiative to facilitate training in peer review and to\nprovide appropriate recognition for Early CareerResearchers who co-review manuscripts.\nResponse: All comments from the reports of Reviewer 1 and Reviewer 2 were very useful\nto us and helped to improve the presentation and clarity of our results. Once again, we would\nlike to thank Reviewer 3 for his/her contribution to one of the above reports.\nResponse to Reviewer 4:\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports.\nThis is part of the NatureCommunications initiative to facilitate training in peer review and to\nprovide appropriate recognition for Early CareerResearchers who co-review manuscripts.\nResponse: All comments from the reports of Reviewer 1 and Reviewer 2 were very useful\nto us and helped to improve the presentation and clarity of our results. Once again, we would\nlike to thank Reviewer 4 for his/her contribution to one of the above reports."}]

Decision letter