Inferring connectivity of an oscillatory network via the phase dynamics reconstruction

1 Introduction

1.1 The connectivity problem

Inferring network connectivity from observation represents one of the most challenging data analysis problems. This task finds practical applications in many fields and particularly in life sciences. For example, determining brain area connectivity is essential for studying normal and pathological brain function. This brief review summarizes one of the existing approaches to the connectivity problem. This approach is not general: it applies to the case of oscillatory networks when each node represents an active, self-sustained oscillator. However, as we argue below, the approach is natural for this case, and the results admit a clear interpretation. The technique relies on the phase dynamics reconstruction from observations; it assumes that the outputs of all network nodes are available.

Before proceeding with the approach’s description and discussion, we briefly formulate how we understand the connectivity. We start by formulating an ensemble of N coupled dynamical systems, the couplings are not necesseraly pairwise, but can include many-body interactions. Suppose the dynamics of the kth node are described by ${\dot{\mathbf{x}}}_k={\mathbf{F}}_k({\mathbf{x}}_k)+\epsilon {\mathbf{H}}_k({\mathbf{x}}_k,{\mathbf{x}}_1,\dots ,{\mathbf{x}}_N)$, where the vector x_k consists of the state variables at the node k (below we focus on the case of weak coupling, thus we introduce a parameter ɛ having the dimension of frequency). Notice that even a more generic setup could be considered, where couplings themselves are described by some differential equations. Vector functions F_k and H_k determine the autonomous and the interaction dynamics of the node k. Since the numbering of oscillators is arbitrary, for the sake of brevity it is convenient to write, without loss of generality, the dynamics equation for the first unit, i.e., to set k = 1 and omit one index. Thus, we write

$\dot{\mathbf{x}}=\mathbf{F}\left(\mathbf{x}\right)+\epsilon \mathbf{H}\left(\mathbf{x},{\mathbf{x}}_2,\dots ,{\mathbf{x}}_N\right).(1)$

This equation and its phase approximation version written below help to quantify all incoming links 1 ← j, j = 2, … , N. If the coupling term H depends explicitly¹ on x_j, we say that there exists a structural link j → 1. In other words, a structural link means some physical connection, e.g., resistive coupling of electronic circuits, synaptic or gap-junction neuronal coupling, connection of brain areas via fibers, etc. Such connections are represented by an explicit dependence in the dynamical equations, which is generally directional (i.e., existing connection j → 1 does not imply the link 1 → j and vice versa; if both links exist, their strengths generally differ). Naturally, node 1 can have many incoming links. In this case, the function H depends on many arguments. Below, we mainly concentrate on the typical case of pair-wisely coupled networks when

$\mathbf{H}={\displaystyle \sum _{j\ne 1}}{\mathbf{H}}_j\left(\mathbf{x},{\mathbf{x}}_j\right),(2)$

but in Section 2.3.4, also many-body (triplet) couplings will be discussed shortly.

Even if the nodes k, j are not structurally connected, they may exhibit correlated dynamics due to a common drive or indirect coupling via a third node, etc. Different measures partially discussed below quantify the degree of correlation or, in other words, the functional connectivity. To illustrate the difference between structural and functional connectivities, consider a simple motif of three oscillators, where the second unit drives the other two, 1 ← 2 → 3. Obviously, units 1 and 3 are not structurally linked but may be correlated, i.e., functionally connected.

Most important for us is the notion of effective phase connectivity because it is precisely what the approach discussed below yields. A characterization of oscillations with their phases (to be discussed in more details in Section 1.2 below) replaces dynamical variables at each node x_k with the phases φ_k. Generally, the phase dynamics equation of the first node reads

$\dot{\phi }=\omega +\epsilon h\left(\phi ,{\phi }_2,\dots ,{\phi }_N\right),(3)$

where ω is the natural frequency, and h is the (phase) coupling function. If h depends on φ_j, we say that there is an effective phase connection j → 1. We discuss the relation between functions H and h and the relation between structural and effective phase connectivities in the next section after recalling the basic results of the phase reduction theory. We emphasize that in the context of the connectivity problem, the word “effective” sometimes means the data-based estimation of the true connectivity. We use the term to describe the connectivity as the true phase dynamics give it; as we argue below, it generally differs from structural connectivity.

1.2 Basic facts about phase reduction

As already emphasized, the reviewed approach applies to networks of self-sustained units. Although natural oscillators are inevitably noisy or weakly chaotic, many analysis techniques rely on the phase reduction theory, which is mostly easily formulated for coupled limit-cycle oscillators (and then extended under some approximations to other cases). In this Section, we briefly present the main results of this theory. Thus, here we assume that isolated systems possess stable limit cycles and, hence, exhibit periodic oscillations. It is well-known that for such systems the phase variable can be introduced obeying

${\dot{\phi }}_k={\omega }_k.(4)$

It is important that the phase is defined both for the limit cycle and its basin of attraction. The description of the dynamics in terms of the phases means an immediate reduction of the problem dimensionality: each limit-cycle oscillator, two- or multidimensional, is now quantified by only one variable φ_k. The dynamics of interacting units obeys Eq. 3 and the phase space of the system is the N-dimensional torus, unless the coupling becomes too strong and destroys the torus (Afraimovich and Shilnikov, 1983). We underline that derivation of functions h from given H remains an unsolved problem in the general case. The well-known theory provides a recipe for the derivation in the first-order approximation in the interaction strength, and even in this case the analytical solution requires knowledge of phase sensitivity functions [see detailed discussion by Nakao (2016)]. However, for network reconstruction from data, this derivation is not needed. We only need the general properties of functions h.

Let us assume that the coupling strength ɛ is small compared to the frequency ɛ ≪ ω, and express the function h in powers of ɛ. Eq. 3 becomes

$\dot{\phi }=\omega +\epsilon Q_1\left(\phi ,{\phi }_2,\dots ,{\phi }_N\right)+{\epsilon }^2{Q}_2\left(\phi ,{\phi }_2,\dots ,{\phi }_N\right)+\cdots .(5)$

In the limit of weak coupling, one keeps only the first-order term Q₁, and the general theory says, that for the pair-wise coupling, see Eq. 2, this term reads

$Q_1\left(\phi ,{\phi }_2,\dots ,{\phi }_N\right)={\displaystyle \sum _{j>1}}{q}_j\left(\phi ,{\phi }_j\right).(6)$

It means that in the weak-coupling limit, the pair-wise coupling on the level of full equations results in the pair-wise coupling in the phase model. This property is not preserved if the high-order terms, i.e., terms ∼ɛ², ∼ɛ³, … are not negligible. So, the analysis of the analytically solvable case of three oscillators coupled pair-wisely in a chain (i.e., only structural couplings 1 ↔ 2 and 2 ↔ 3 are present) yields the second-order term Q₂(φ, φ₂, φ₃) in the equation for the first unit dependent on all three phases (Gengel et al., 2021), though the first and third oscillators are not structurally linked: they interact only through the second oscillator. This theoretical result confirms earlier numerical observations (Kralemann et al., 2011; Kralemann et al., 2014) demonstrating the difference between the structural and effective phase connectivity. We expect that third-order terms Q₃ provide dependence on the phase of the lth unit interacting with the first oscillator indirectly through two mediators, e.g., in the following configuration, 1 ↔ n ↔ m ↔ l, and so on.

A remark is in order. Consider for a moment two oscillators. The first-order phase approximation yields a general function of two phases, Q₁ = Q₁(φ₁, φ₂). If the coupling is weak ɛ ≪ ω, then one typically performs another approximation step, by averaging the phase equation over the oscillation period 2π/ω and obtaining the averaged coupling function (describing slow variations of the phases) that depends on the slow phase difference only, ${\tilde{Q}}_1={\tilde{Q}}_1({\phi }_1-{\phi }_2)$ (the Kuramoto-Daido form), see, e.g., Pikovsky et al. (2001). This step essentially facilitates the theoretical treatment of the phase dynamics, but using the Kuramoto-Daido form instead of the general dependence may be an unnecessary simplification for the connectivity inference from data. Another useful representation of the phase dynamics is the Winfree form. If the term H₂ can be written as H₂(x₁, x₂) = sH₂(x₁, x₂) = sh₂(φ₁, φ₂), where s is a constant unit vector and we substituted the variables x_1,2 on the limit cycles as functions of φ_1,2, then one obtains

$Q_1\left({\phi }_1,{\phi }_2\right)=Z\left({\phi }_1\right)h_2\left({\phi }_1,{\phi }_2\right),(7)$

where Z(φ₁) is the phase sensitivity function for a perturbation in direction s, or phase response curve (PRC). In many cases one assumes that the driving term does not contain the driven coordinates, i.e., h₂(φ₂) contains the driving phase only. In case of more than two units, the sensitivity functions generally differ for different inputs.

All real-world oscillators are noisy, and in the simplest case of a state-independent noise we shall add a random term to the right-hand side of Eqs 1–5. However, the main properties of the phase dynamics remain, also for the case of weakly chaotic oscillators.

In summary, the main idea of the approach reviewed in this paper is to estimate phases from observed time-series data and reconstruct Eq. 3 for each network node. These equations provide effective phase connectivity that is close but not identical to the structural one. Two factors render efficient reconstruction: 1) phase representation is low-dimensional (there is only one equation per unit), and 2) the function h on the right-hand side has a relatively simple form (it can be represented as a multiple Fourier series). This inference is possible if the data are suitable for phase estimation; we discuss this case in the next section. The situation differs if the data look like a sequence of spikes. We address this problem in Section 3.

2 Case I: smooth oscillatory dynamics

2.1 From time series to phases

2.1.1 From time series to a protophase

The first step in the phase dynamics analysis is the phase inference from the observed signals. Generally, observables are functions of the state variables x. For periodic dynamics (i.e., without coupling) these observables are periodic functions of time, but under coupling, these observables are nearly periodic signals with amplitude and phase modulations. An extraction of the phase is simple, if at least two scalar observables $y_k^{(1)}(t),y_k^{(2)}(t)$ for a system k are available. In this case, a two-dimensional plot of a trajectory on the plane (y⁽¹⁾, y⁽²⁾) is a nearly closed curve, and one can characterize the phase by a coordinate along the curve. In the simplest case when the curve is close to a circle, one can take ${\theta }_k=\text{arg}[y_k^{(1)}-⟨y_k^{(1)}⟩+\mathrm{i}(y_k^{(2)}-⟨y_k^{(2)}⟩)]$ as such a coordinate. If the curve is more complex, e.g., it has several loops, one can define this coordinate as a normalized length along the curve,

$\theta =2\pi \frac{L\left(t\right)-L\left(t_i\right)}{L\left(t_{i+1}\right)-L\left(t_i\right)}$

(here t_i are time instants at which the ith revolution completes; each revolution corresponds to one oscillation). This approach ensures that the obtained coordinate θ grows monotonically with time. A complex waveform appears, e.g., in the analysis of an electrocardiogram, where one heartbeat cycle consists of several “waves.” For an advanced algorithm developed for the ECG analysis, see Kralemann et al. (2013a).

The problem is less trivial if only one scalar observable y_k(t) is available. To perform the embedding, one needs to construct from it the second observable. The method of choice in the literature is obtaining the second time series by virtue of the Hilbert Transform (King, 2009; Feldman, 2011): ${\tilde{y}}_k(t)=\widehat{H}[y_k(t)]$. The variable θ based on the embedding $(y_k(t),{\tilde{y}}_k(t))$ is often called “Hilbert phase” (quite often it is additionally assumed “by default” that θ is an argument of a complex observable $y_k(t)+\mathrm{i}\widehat{H}[y_k(t)]$). This method works reasonably well for a nearly harmonic signal with a weak and slow modulation of the phase and small amplitude modulation. If the phase modulation is not slow and small, the HT method is inaccurate but still used in many applications because of its simplicity. Recently, an improvement of the HT method has been suggested (Gengel and Pikovsky, 2019; Gengel and Pikovsky, 2021; Gengel and Pikovsky, 2022). The procedure implies an iterative approach: one takes the variable θ⁽¹⁾ based on the embedding $(y_k(t),\widehat{H}[y_k(t)])$ as a first approximation (thus the superscript 1), and uses it as a new time variable to perform a new embedding $(y_k({\theta }^{(1)}),\widehat{H}[y_k({\theta }^{(1)})])$, where now in the computation of HT one considers the signal as a function of θ⁽¹⁾ instead of a function of time. This step begets a new, improved coordinate θ⁽²⁾, etc. Gengel and Pikovsky (2019), Gengel and Pikovsky (2021), and Gengel and Pikovsky (2022) demonstrated that for a signal with a pure phase modulation, this iterative procedure converges, allowing for a perfect inference of a phase. However, if amplitude modulation is also present, the inference is inaccurate. Developing a reliable method for phase and amplitude demodulation of a signal remains a challenging problem; see also Matsuki et al. (2023).

The variable θ discussed above is monotonic in time and grows by 2π at each oscillation, but it is not the phase because the genuine phase of an oscillator has a special property—it grows, according to Eq. 4, uniformly in time, while the variable θ obeys

$\dot{\theta }=\omega +f\left(\theta \right),(8)$

where the function f(θ) depends on the form and selection of the observables y and on the method of its extraction. To emphasize this difference, variable θ is called protophase.

2.1.2 From a protophase to the phase

The transformation from a monotonic protophase to the true phase is a necessary and straightforward step in the phase reduction theory². Writing $\frac{\mathrm{d}\phi }{\mathrm{d}t}=\omega =\frac{\mathrm{d}\phi }{\mathrm{d}\theta }\frac{\mathrm{d}\theta }{\mathrm{d}t}$, one obtains the desired transformation

$\frac{\mathrm{d}\phi }{\mathrm{d}\theta }=\omega {\left[\frac{\mathrm{d}\theta }{\mathrm{d}t}\right]}^{-1}=\sigma \left(\theta \right),\text{or}\phi ={\int }_0^{\theta }\sigma \left({\theta }^{\prime }\right)\mathrm{d}{\theta }^{\prime }.(9)$

Finding the transformation function σ(θ) in data analysis is non-trivial: we must estimate it from a noisy trajectory. Hence, we must compute an average of dt/dθ over the trajectory. Kralemann et al. (2007) and Kralemann et al. (2008) demonstrate an efficient algorithm that avoids numerical differentiation. Let the discretely sampled protophases be ${\widehat{\theta }}_j$, j = 1, … , N_p, where N_p is the number of points. Then the transformation reads:

$\phi =\theta +2{\displaystyle \sum _{n=1}^M}\text{Im}\left[S_n\left(e^{\mathrm{i}n\theta }-1\right)/n\right],\text{where}{S}_n=N_p^{-1}{\displaystyle \sum _{j=1}^{N_p}}{e}^{\mathrm{i}n{\widehat{\theta }}_j}.(10)$

Here M should not be chosen too high to avoid overfitting. Numerical tests by Kralemann et al. (2008) demonstrate that the protophase-to-phase transformation (Eq. 10) provides a nearly homogeneously growing phase. It removes the deviations from ωt with the characteristic time 2π/ω, but preserves low-frequency features like phase diffusion. It is important that the transformation (Eq. 10) is invertible and therefore not a “filtering” or “smoothing.”

2.2 Undirected connectivity

The main qualitative effect of weak interaction between self-oscillatory systems is a possibility of their synchronization manifested by frequency and phase locking. So, for two periodic oscillators, which, being uncoupled, have close but different frequencies, synchronization means that the frequencies become precisely equal if the interaction is switched on and its strength exceeds some threshold. Equality of frequencies means boundness of the phase difference ψ = φ₁ − φ₂, i.e., |ψ| < const. Generally, if the frequencies of uncoupled systems fulfill nω₁ ≈ mω₂, where n, m are some positive integers, the interaction may result in synchronization, which can be defined according to frequency and phase locking conditions, nΩ₁ = mΩ₂ and |ψ_n,m| = |nφ₁ − mφ₂| < const, respectively; here, Ω_1,2 denote frequencies of coupled systems (observed frequencies). In a noisy environment, ψ_n,m remains constant for long time intervals, but can exhibit relatively rapid ± 2π jumps, called phase slips. Then, the locking condition |ψ| < const does not hold, but the distribution of ψ_n,m mod 2π remains non-uniform. Note that phase slips also appear in a noise-free case when the oscillators are slightly outside the border of synchrony. A measure of the non-uniformity of the distribution of the phase differences quantifies thus phase coherence. In other words, it quantifies how far the system is from the asynchronous case of two uncoupled oscillators where the distribution of ψ_n,m mod 2π is uniform.

The corresponding measure can be the entropy of the distribution or the amplitude of its first Fourier mode. The latter is known as the phase locking value or synchronization index and is most popular because it has no parameters (Rodriguez et al., 1999; Mormann et al., 2000; Rosenblum et al., 2001). For the general case of n: m locking, where n, m are some positive integers, this measure is

${\gamma }_{n,m}=|⟨\exp \left[\mathrm{i}\left(n{\phi }_1-m{\phi }_2\right)\right]⟩|,(11)$

where ⟨⋅⟩ means averaging over time. Obviously, 0 ≤ γ_n,m ≤ 1; this index quantifies the constancy of the phase difference. Since the ratio between the frequencies of the uncoupled systems is unknown, though it can be roughly estimated, one tries different combinations of n, m and picks the values which maximize γ_n,m.

For a network, computing the index γ_n,m, we can assign its value to the link between these nodes. This procedure provides an undirected network of functionally connected units but does not yield reliable information about the underlying physical connections. To illustrate this issue, let us consider a motif of three units coupled in the following configuration: 1 ↔ 2 ↔ 3. It is known that the 1st and 3rd oscillators can synchronize, while the central unit 2 remains asynchronous. This state is denoted as remote synchrony, see Bergner et al. (2012); its appearance can be explained by the second-order phase reduction that yields the coupling term ∼ɛ² dependent on φ₁ − φ₃ (Kumar and Rosenblum, 2021). For this state, quantifying functional interrelations through the index γ will yield a strong link 1 ↔ 3 and no links for the pairs (1, 2) and (2, 3).

With Figure 1 we illustrate the importance of the protophase-to-phase transformation for computation of γ. See Kralemann et al. (2008) for analytical examples demonstrating that calculating the index from protophases can yield either an under- or overestimated value. Here, we generate the signals by simulating the dynamics of two coupled van der Pol oscillators:

${\ddot{x}}_{\mathrm{1,2}}-3\left(1-x_{\mathrm{1,2}}^2\right){\dot{x}}_{\mathrm{1,2}}+{\left(1\pm \mathrm{\Delta }\omega \right)}^2{x}_{\mathrm{1,2}}=0.1\left({\dot{x}}_{\mathrm{2,1}}-{\dot{x}}_{\mathrm{1,2}}\right).(12)$

FIGURE 1

FIGURE 1. Synchronization index as a function of frequency mismatch in a system of two coupled van der Pol oscillators, see Eq. 12. The solid black line shows the index computed from the true phases; this computation requires knowledge of the systems equations and serves as a benchmark. The vertical dashed line marks the border of the synchronization domain. We see that computing γ from protophases 1) yields the protophase-dependent results, and 2) the index value can be essentially smaller than one for synchronous states. Next, we see that the protophase-to-phase transformation provides an essential improvement.

Taking x_1,2 as observables, we compute the Hilbert protophases ${\theta }_{\mathrm{1,2}}^{(H)}$ and length protophases ${\theta }_{\mathrm{1,2}}^{(L)}$. Then we compute γ_1,1 from protophases and transformed protophases. For the benchmark, we take the index calculated from true phases³. The results demonstrate that the index γ depends on the chosen protophase, while the protophase-to-phase transformation eliminates this ambiguity and crucially improves the estimate. Finally, we mention that the index is large outside of the synchronization domain, which means that it is not a measure of synchronization but of the interaction strength reflected in phase correl ations.

The definition of the synchronization index can be easily extended to quantify triplet interactions. Consider triplet locking defined via the condition

$|n{\phi }_1+m{\phi }_2+l{\phi }_3|<\text{const},n{\omega }_1+m{\omega }_2+l{\omega }_3=0,(13)$

where integers n, m, l can be both positive and negative, while the conditions of the pairwise locking are not satisfied for any pair of units (Kralemann et al., 2013b). The corresponding index is

${\gamma }_{n,m,l}=|⟨\exp \left[\mathrm{i}\left(n{\phi }_1+m{\phi }_2+l{\phi }_3\right)\right]⟩|.(14)$

We emphasize that many-body interaction when coupling terms depend on more than two phases naturally appears in oscillatory networks, see Pikovsky and Rosenblum (2022) for a review and Nijholt et al. (2022) for an experimental illustration.

A natural question is: what is the advantage of index γ_n,m over the usual correlation directly computed from the time series without phase estimation? The reason is twofold. First, the synchronization index distinguishes the locking states (or, more precisely, states close to locking) of different orders. Second, if the oscillators are noisy or weakly chaotic, the amplitudes may strongly fluctuate while the phases can lock. Thus, a measure exploiting phases is more sensitive to interaction.

2.3 Directed connectivity

Determination of the true structural connectivity of a network requires knowledge of the state space Eq. 1. One can compute the norm of the corresponding coupling function H as a measure of the connection’s strength. So the effect of the jth node on the first one is quantified by ‖H(x, x_j)‖. However, such quantification is not complete since the effect of the forcing depends not only on the force’s strength but also on the system’s sensitivity to this influence. Anyway, reconstruction of Eq. 1 from scalar data is hardly feasible. In contradistinction, inference of the phase dynamics Eq. 3 is possible; some limitations are discussed below. In the rest of this section, we treat separately the cases of two, three, and more than three units.

2.3.1 Two oscillators

Suppose we observe two oscillators and register two observables y_1,2. Suppose also that these observables are suitable for phase estimation. Then, we perform a two-step procedure for each observable: first, we compute protophases exploiting an appropriate technique and then execute the protophase-to-phase transformation. The next step is to estimate phase derivatives (the instantaneous frequencies). To this end, we first unwrap each phase to make it a monotonically growing function of time and then apply a Savitzky–Golay filter, i.e., for each time point, perform a local polynomial fit in a running window, and obtain the derivative at this point from the fitted curve. Local fitting is a smoothing filter, making this approach especially useful for noisy data. Having phases and their derivatives, we find the right-hand sides of coupled equations

${\dot{\phi }}_{\mathrm{1,2}}={\omega }_{\mathrm{1,2}}+Q^{\left(\mathrm{1,2}\right)}\left({\phi }_{\mathrm{1,2}},{\phi }_{\mathrm{2,1}}\right)$

(here, it is convenient to write the index of the first system explicitly)⁴. There are at least two ways to achieve this. The first one is to use the fact that the coupling function is 2π-periodic in both arguments, represent it by a double Fourier series, and find the coefficient of this series by the least mean squares fit. The second one is to use the kernel density estimation to obtain the desired function on a rectangular grid. Note that practically we compute the r.h.s. of phase equations and then represent it as a sum of the constant non-zero term (frequency ω) and zero-mean function⁵.

Now, we introduce an index d_1→2 quantifying the directionality of coupling (Rosenblum and Pikovsky, 2001). It is defined so that d_1→2 = 1 for unidirectional coupling from the first to the second unit, d_1→2 = −1 for the opposite case, and −1 ≤ d_1→2 ≤ 1 for bidirectional coupling. The directionality index is

$d_{1\to 2}=\frac{c_2-c_1}{c_1+c_2},$

where c_1,2 = ‖Q^(1,2)‖/ω_1,2 and $\Vert Q^{(\mathrm{1,2})}{\Vert }^2=\sum |{\mathcal{F}}_{n,m}^{(\mathrm{1,2})}{|}^2$, where ${\mathcal{F}}_{n,m}^{(\mathrm{1,2})}$ are the Fourier coefficients of the corresponding zero-mean function Q^(1,2). Coefficients c_1,2 provide a relative dimensionless measure of the external action on a unit. For effects of noise and data length on the efficiency of this technique, see Smirnov and Bezruchko (2003). Comparison of this approach with the partial directed coherence can be found in Smirnov et al. (2007).

We conclude the discussion of the two-units case with two remarks. 1) Inference of a function of two variables φ₁, φ₂ requires that the given points scatter over the square domain (toroidal surface) 0 ≤ φ_1,2 < 2π. In other words, the systems shall not synchronize. In the opposite case, the trajectory is just a curve on the surface of the torus, and the fitting procedure fails. 2) A widely used approach is reconstructing the coupling function as a function of the phase difference, i.e., in the Kuramoto-Daido form. Such inference requires less numerical effort than recovering the general function of two phases, but the price is that one has to choose the n: m ratio beforehand.

2.3.2 Three oscillators

For a motif of three oscillators, we first calculate three protophases, then perform three times the transformation to phases, and infer the coupling functions exploiting the least mean squares fit⁶. To be exact, by fitting, one finds the coefficients of the Fourier representation, e.g., for the first oscillator:

${\dot{\phi }}_1={\displaystyle \sum _{l_1,l_2,l_3}}{\mathcal{F}}_{l_1,l_2,l_3}\exp \left[\mathrm{i}\left(l_1{\phi }_1+l_2{\phi }_2+l_3{\phi }_3\right)\right].$

Then we set $\omega ={\mathcal{F}}_{\mathrm{0,0,0}}$ and assign the sum of all terms except for the constant one to Q⁽¹⁾. To quantify a link, e.g., a pairwise action from 2 → 1, we compute a partial norm

${\mathcal{P}}_{1\leftarrow 2}^2={\displaystyle \sum _{l_1,l_2\ne 0}}|{\mathcal{F}}_{l_1,l_2,0}{|}^2.(15)$

Quantification of the joint action of units 2 and 3 on the first one is given by

${\mathcal{T}}_{1\leftarrow \mathrm{2,3}}^2={\displaystyle \sum _{l_1,l_2\ne 0,l_3\ne 0}}|{\mathcal{F}}_{l_1,l_2,l_3}{|}^2.(16)$

We emphasize that one can perform a pairwise analysis, e.g., reconstructing the coupling function Q⁽¹⁾ from phases φ₁, φ₂ as described in the previous section. In this way, one ignores the third oscillator. The norm of the reconstructed function ‖Q⁽¹⁾‖ generally differs from ${\mathcal{P}}_{1\leftarrow 2}$; indeed, the former one captures both direct and indirect causal effect from 2 to 1, while the latter is less sensitive to indirect effects.

2.3.3 More than three oscillators

An extension to the case of N > 3 units seems straightforward. However, by reconstructing a high-dimensional coupling function, we face the curse of dimensionality. For a successful inference, we must fill the N-dimensional hypercube, which requires large data sets. A possible solution is to use triplet analysis (Kralemann et al., 2014). Suppose we want to quantify a link from unit k to unit j. To this goal, we reconstruct the equation

${\dot{\phi }}_j={\omega }_j+Q^{\left(j,k,m\right)}\left({\phi }_j,{\phi }_k,{\phi }_m\right)$

for all m ≠ j, k, i.e., for all possible triplets. Then, from each triplet we obtain partial norm

${\tilde{\mathcal{T}}}_{j\leftarrow k}^2\left(m\right)={\displaystyle \sum _{l_j,l_k\ne 0}}{\left|{\mathcal{F}}_{l_j,l_k,0}^{\left(j\right)}\right|}^2(17)$

and take

${\mathcal{T}}_{j\leftarrow k}^2=\underset{m}{\min }{\tilde{\mathcal{T}}}_{j\leftarrow k}^2\left(m\right)(18)$

as the final triplet-based measure of the link strength. Figure 2 illustrates the idea.

FIGURE 2

FIGURE 2. Triplet analysis of networks with N > 3 oscillators. Suppose we quantify the link 1 → 3 absent in the given configuration. However, if we estimate the link strength using Eq. 17 exploiting the triplet (3, 1, 5), shown by the orange color in (A), we obtain a non-zero value, indicated by a dashed line. The explanation is that units 1 and 3 are connected through unit 2, and this connection is not captured by the triplet used. We obtain a much better estimation using the triplet (3, 1, 2) (B). Notice that the resulting measure is always positive. Therefore, taking minimum over all triplets, see Eq. 18, yields the best estimation.

Numerical experiments (Kralemann et al., 2014) with small networks (N = 5 and N = 9) demonstrate that the triplet analysis provides a very good separation between existing and non-existing connections. Moreover, the results confirm that most of the non-existing links revealed by the technique are not artifacts but reflect causal information flow mediated by indirect driving. This driving is due to terms emerging in the high-order phase reduction.

The technique was thoroughly tested by Rings and Lehnertz (2016) for networks of weakly chaotic Rössler oscillators, also in the presence of noise. The main result is that, for networks with high mean degrees and a larger number of nodes (N ≫ 10), the triplet analysis does not appear to provide information about directional couplings other than the one obtained with the pairwise analysis.

2.3.4 Reconstruction of large networks

The generalization of the presented techniques to the case of large networks is straightforward. So, Pikovsky (2018) has treated three basic models of coupling in large networks (below, we assume that the initial preprocessing, including protophase extraction and a transformation to the phase, is completed):

${\dot{\phi }}_k={\omega }_k+{\displaystyle \sum _{j=1;j\ne k}^N}{\mathrm{\Gamma }}_{kj}\left({\phi }_j-{\phi }_k\right),(19)$

${\dot{\phi }}_k={\omega }_k+{\displaystyle \sum _{j=1;j\ne k}^N}{Q}_{kj}\left({\phi }_j,{\phi }_k\right),(20)$

${\dot{\phi }}_k={\omega }_k+{\displaystyle \sum _{j=1;j\ne k}^N}{\displaystyle \sum _{l=1,l\ne k,l>j}^N}{G}_{\mathit{kjl}}\left({\phi }_j,{\phi }_k,{\phi }_l\right).(21)$

Here Eq. 19 describes a so-called Kuramoto-Daido network with arbitrary and generally different pairwise coupling functions depending on the phase differences; Eq. 20 describes a network where pairwise coupling functions are general 2π-periodic functions of the phases; Eq. 21 describes a hypernetwork with triplet interactions.

As the first step, we calculate from the time series of phases their time derivatives ${\dot{\phi }}_k$ [the most common approach here is to use Savitzky-Golay filtering, see Ahnert and Abel (2007)]. In the next step, we represent the coupling functions as sums of elementary sin and cos functions (in the case of the Kuramoto-Daido coupling) or as sums of products of these functions. For example, the coupling function Q_kj in Eq. 20 is represented as

$Q_{kj}\left({\phi }_j,{\phi }_k\right)={\displaystyle \sum _{m=1}^{2M}}{\displaystyle \sum _{l=1}^{2M+1}}{q}_{\mathit{kjml}}{f}_m\left({\phi }_j\right)f_l\left({\phi }_k\right),(22)$

where

$f_{2n-1}\left(x\right)=\cos nx,f_{2n}\left(x\right)=\sin nx,f_{2M+1}=1,n=1,\dots ,M.$

Parameter M defines the number of harmonics in the test functions; it should be adjusted according to complexity of the coupling function. Exploiting this representation, one performs minimization of the squared error for the discrete time series φ_k(n)

${\displaystyle \sum _n}{\left({\dot{\phi }}_k-{\tilde{\omega }}_k-{\displaystyle \sum _j}{Q}_{kj}\left({\phi }_j,{\phi }_k\right)\right)}^2,$

where Q_kj(φ_j, φ_k) is substituted using Eq. 22, and then finds unknown coefficients q and frequencies ${\tilde{\omega }}_k$ by virtue of singular value decomposition. We illustrate the performance of the approach for the Winfree-type model with random coupling functions (containing p = 3 harmonics) and a fully connected random network of 16 oscillators in Figure 3.

FIGURE 3

FIGURE 3. Results of the reconstruction of the coupling constants for one oscillator in the random network (Eq. 20) using 5,000 observation points. Norms of the reconstructed coupling functions are plotted vs. the true ones. Square, plus, and circle symbols correspond to M = 1, M = 2, and M = p = 3, respectively. Dashed lines help to see validity of a linear relation between the reconstructed and true norms. For the proper number M = 3 of explored harmonics the reconstruction is very good.

Similar approaches to reconstruction of large oscillatory networks have been discussed by Casadiego and Timme (2018), Panaggio et al. (2019), Tokuda et al. (2019), Novaes et al. (2021), and Rings et al. (2022). We mention that there is still a problem of discriminating small couplings from nonexisting ones. Generally, one obtains false positive and false negative interactions. See Cecchini et al. (2021) for an extended analysis of relations of these errors to the network structure.

3 Case II: pulse-like oscillatory dynamics

Suppose we observe neuronal spiking, or the data we measure looks like spiking, so reducing the signals to point processes is appropriate. If it is reasonable to assume the self-sustained activity, then the proper model is that of pulse-coupled integrate-and-fire units⁷. In this one-dimensional model, the phase of each unit grows uniformly, ${\dot{\phi }}_k={\omega }_k$. When the phase φ_k attains 2π, the unit issues a spike, and its phase is instantaneously reset to zero.

3.1 Undirected connectivity

Suppose we deal with two units and know the times at which they fire. Without loss of generality, we say that phase is zero when a unit spikes, and it grows to 2π within an inter-spike interval. The simplest way to quantify the degree of synchrony for two units is to approximate their phases linearly between the spikes and compute the synchronization index. For advanced techniques focused on analysis of spike train data, see Kreuz et al., 2007, Kreuz, 2011, and Satuvuori et al., 2017.

3.2 Directed connectivity

Here, we discuss the connectivity inference for the case when a network of pulse-coupled oscillators can reasonably model the system, and the available data represent sequences of spikes (Cestnik and Rosenblum, 2017). We assume that the coupling is sufficiently weak to justify the phase description. Next, we suppose that a unit’s phase response curve (PRC) is the same for all incoming connections, though PRCs of different units generally differ. We describe individual units by the integrate-and-fire model. It is, without interaction, the phase of each unit grows uniformly, ${\dot{\phi }}_k={\omega }_k$. When the phase φ_k attains 2π, the unit issues a spike, and its phase is instantaneously reset to zero. At this moment, the unit k sends a stimulus through all outgoing links and thus kicks all units connected to it through these links. When unit j receives a kick from unit k, its phase is instantaneously reset according to its PRC Z_j(φ):

${\phi }_j\to {\phi }_j+{\epsilon }_{jk}{Z}_j\left({\phi }_j\right),$

where coefficient ɛ_jk quantifies the strength of the link k → j. Generally, ɛ_jk ≠ ɛ_kj.

In our approach, we choose one oscillator (let it be the first one) and consider all its incoming connections. To simplify the notations, we omit the index 1, so for the incoming links we rename ɛ_1k by ɛ_k, with k = 2, 3, … , N. Using an iterative procedure described below, we infer this oscillator’s frequency ω, PRC Z, and ɛ_k. Next, we perform the same analysis for all other units.

We denote the inter-spike intervals for the chosen (first) unit as T_m, where m = 1, 2, … , M and M + 1 is the number of spikes in the pulse train of the first unit. We write the following equations for the phase evolution within each interval T_m:

$\omega T_m+{\displaystyle \sum _{k=2}^N}{\epsilon }_k{\displaystyle \sum _{l=1}^{n_m\left(k\right)}}Z\left({\phi }_m^{\left(k,l\right)}\right)=2\pi .(23)$

Here ${\phi }_m^{(k,l)}$ is the phase of the first unit at the instant when it receives the lth spike from the unit k, within the inter-spike interval number m; n_m(k) is the number of incoming stimuli from the unit k within interval T_m. For sufficiently large M, these equations can be solved by iterations to provide unknown ω, ɛ_k, Z, see Cestnik and Rosenblum (2017).

The key idea for solving Eq. 23 is as follows. Suppose for a moment that we know phases ${\phi }_m^{(i,l)}$ and coupling coefficients ɛ_k. Representing Z(ϕ) as a finite Fourier series, we obtain the overdetermined linear system for unknown Fourier coefficients (provided M > 2N_F + 1, where N_F is the order of the Fourier expansion). Suppose, vice versa, that we know phases and PRC; then we obtain a linear system to find coupling coefficients ɛ_k. These observations suggest a process to solve the problem: we start with some initial estimates for ${\phi }_m^{(i,l)},{\epsilon }_k$ and obtain the first estimates for Z, ω. Next, we use the first estimates for Z, ω to obtain second estimates for ${\phi }_m^{(i,l)},{\epsilon }_k$, etc. Since the coupling coefficients and PRC enter the phase dynamics equation as a product, we obtain these quantities up to a constant factor. The numerical tests show that the iterative procedure converges for the random initial distribution of ɛ_k or equal values ɛ_k = ɛ. The initial phase is taken as linearly growing within the inter-spike interval. The next iterations are piece-wise linear: the phase grows linearly between the incoming spikes and changes instantaneously when these spikes arrive. Suppose we compute the phase at the end of an inter-spike interval for some iteration; we denote this phase as ψ_m,i, where m, i label the inter-spike intervals and the iteration, respectively. ψ_m,i is calculated according to the left-hand side of Eq. 23, using ω, ɛ_k, Z from the previous iteration. Since the latter values are not exact, generally ψ_m,i ≠ 2π. The outcome of this fact is twofold. First, we rescale the phase by 2π/ψ_k,i. Second, we use the standard deviation σ_i = std(ψ_m,i − 2π) as a measure of convergence: if σ_i decreases with iterations, then the results are reliable. Furthermore, one can reconstruct the same network many times, starting from different initial values of the coupling strength and checking the convergence. If different initial values yield close results, then the latter can be trusted.

The tests demonstrate that the technique is robust and capable of dealing with relatively short data (several hundreds of spikes suffice). If the coupling is not weak enough, the estimation of the network connectivity remains correct, but the obtained PRC becomes amplitude-dependent. The method works well unless the nodes synchronize. It may also fail in case of sparse networks where one can expect purely periodic nodes. Indeed, the iterative procedure requires some variability in the drive. However, noise in realistic networks enhances the reconstruction.

4 Discussion

In this mini-review, we presented one line of a broad field of research aiming at network inference from observations. The presented approach explicitly relies on the coupled oscillation theory. It is, therefore, suitable only for the case where all network nodes are self-sustained oscillators, either noisy periodic or weakly chaotic. An additional crucial requirement is that each node is observed and the observables are appropriate for phase estimation. It is not easy to verify these conditions in practice, which is a weak point of the presented approach. However, if the assumptions used are correct, the approach allows for a straightforward interpretation of the results and yields a reasonable inference of the structural connectivity. As important applications, we mention analysis of the interaction between cardiovascular and respiratory systems in humans (Kralemann et al., 2013a) and between cardiac, respiratory, and delta brain rhythms in anesthetized rats (Musizza et al., 2007), as well as use of coupling functions to reveal microvascular impairment in treated hypertension (Ticcinelli et al., 2017). A difficult problem is the nonstationarity of real-world data, see Stankovski et al. (2012) for an estimation of time-evolving coupling. Finally, we mention the assessment of the statistical significance of the inference. The typical approach implies comparing the obtained results with those for the surrogate data, see, e.g., Andrzejak et al. (2003). A discussion of the corresponding techniques may be the subject of a separate review; here, we only mention a simple procedure exploited by Kralemann et al. (2013a) for the bivariate case. To analyze the significance of the coupling function estimation, they compared this function with those obtained for cardiac and respiratory time series taken from different subjects using a similarity index, which quantifies the similarity of forms of two two-dimensional functions.

Another popular approach exploits information-theoretical measures, see, e.g., Hlaváčková-Schindler et al. (2007), Faes et al. (2011), Xiong et al. (2017), and Krakovská et al., 2018. Multiple techniques from this group do not assume that one deals with self-sustained systems and have, therefore, a broad field of applications, e.g., they can be applied to networks consisting of both active and passive units. However, the interpretation of results is complicated. The detected information flows can be interpreted as functional connectivity, which is generally directed, because the information-theoretical measures are asymmetric, contrary to cross-correlations. This functional connectivity generally differs from genuine physical connections (i.e., the inference procedure can yield links that are not structural according to our definition). Finally, we mention a hybrid approach where one evaluates the directionality of interaction by applying the information theory approach to phases (Paluš and Stefanovska, 2003; Vlachos et al., 2022). For a review of techniques originating from different approaches, see Lehnertz (2011), Rings and Lehnertz (2016), Siggiridou et al. (2019), Moraffah et al. (2021), Shojaie and Fox (2022), and Runge et al. (2023). Unfortunately, a detailed comparison of all existing techniques is missing because it requires much effort, publically available codes, and well-designed test models.

Author contributions

MR: Writing–original draft, Writing–review and editing. AP: Writing–original draft, Writing–review and editing.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. They were funded by the eutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Projektnummer 491466077.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The authors declared that they were an editorial board member of Frontiers, at the time of submission. This had no impact on the peer review process and the final decision.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

¹Explicit dependence means that there exists at least one partial derivative ∂H_l/∂x_jm that does not vanish identically, where indices l, m label components of H and x_j, respectively.

²Knowledge of the protophase suffices if only the frequency of the signal is required, since $⟨\dot{\theta }⟩=⟨\dot{\phi }⟩$, where ⟨⋅⟩ denotes averaging over time.

³See Gengel et al. (2021) for the algorithm providing coupled oscillator phases.

⁴We mention that the phase coupling functions Q(1,2) generally contain high-order terms in the coupling strength, cf. Eq. 5.

⁵This representation assumes that Q has zero mean, which is generally false. The coupling function can have a constant term that scales with the coupling strength. If several observations for different though unknown coupling values are available, it is possible to infer the correct value of the natural frequency, see Kralemann et al. (2007) and Kralemann et al. (2008). In the case of only one observation, assigning the constant term of the r.h.s. to ω seems to be the only option.

⁶Kernel density estimation in three dimensions is computationally inefficient and difficult to interpret.

⁷Integrate-and-fire system is the simplest representation of a relaxation oscillator.

References

Afraimovich, V. S., and Shilnikov, L. P. (1983). Invariant two-dimensional tori, their destroying and stochasticity Methods of Qualitative Theory of Differential Equations (Gorki). 3–28. In Russian; English translation Amer. Math. Soc. Transl. Ser. 2 (149), 201–212. 1991.

Google Scholar

Ahnert, K., and Abel, M. (2007). Numerical differentiation of experimental data: local versus global methods. Comput. Phys. Commun. 177, 764–774. doi:10.1016/j.cpc.2007.03.009

CrossRef Full Text | Google Scholar

Andrzejak, R. G., Kraskov, A., Stögbauer, H., Mormann, F., and Kreuz, T. (2003). Bivariate surrogate techniques: necessity, strengths, and caveats. Phys. Rev. E 68, 066202. doi:10.1103/PhysRevE.68.066202

PubMed Abstract | CrossRef Full Text | Google Scholar

Bergner, A., Frasca, M., Sciuto, G., Buscarino, A., Ngamga, E. J., Fortuna, L., et al. (2012). Remote synchronization in star networks. Phys. Rev. E 85, 026208. doi:10.1103/PhysRevE.85.026208

PubMed Abstract | CrossRef Full Text | Google Scholar

Casadiego, J., and Timme, M. (2018). Accelerated reference frames (arfs) reveal networks from time series data. New J. Phys. 20, 113031. doi:10.1088/1367-2630/aaebb8

CrossRef Full Text | Google Scholar

Cecchini, G., Cestnik, R., and Pikovsky, A. (2021). Impact of local network characteristics on network reconstruction. Phys. Rev. E 103, 022305. doi:10.1103/PhysRevE.103.022305

PubMed Abstract | CrossRef Full Text | Google Scholar

Cestnik, R., and Rosenblum, M. (2017). Reconstructing networks of pulse-coupled oscillators from spike trains. Phys. Rev. E 96, 012209. doi:10.1103/PhysRevE.96.012209

PubMed Abstract | CrossRef Full Text | Google Scholar

Faes, L., Nollo, G., and Porta, A. (2011). Information-based detection of nonlinear Granger causality in multivariate processes via a nonuniform embedding technique. Phys. Rev. E 83, 051112. doi:10.1103/PhysRevE.83.051112

PubMed Abstract | CrossRef Full Text | Google Scholar

Feldman, M. (2011). Hilbert transform applications in mechanical vibration. John Wiley & Sons.

Google Scholar

Gengel, E., and Pikovsky, A. (2019). Phase demodulation with iterative Hilbert transform embeddings. Signal Process. 165, 115–127. doi:10.1016/j.sigpro.2019.07.005

CrossRef Full Text | Google Scholar

Gengel, E., and Pikovsky, A. (2021). “Phase reconstruction with iterated Hilbert transforms,” in Physics of biological oscillators. Editors A. Stefanovska,, and P. McClintock (Springer), 191–208.

CrossRef Full Text | Google Scholar

Gengel, E., and Pikovsky, A. (2022). Phase reconstruction from oscillatory data with iterated Hilbert transform embeddings – benefits and limitations. Phys. D. Nonlinear Phenom. 429, 133070. doi:10.1016/j.physd.2021.133070

CrossRef Full Text | Google Scholar

Gengel, E., Teichmann, E., Rosenblum, M., and Pikovsky, A. (2021). High-order phase reduction for coupled oscillators. J. Phys. Complex. 2, 015005. doi:10.1088/2632-072X/abbed2

CrossRef Full Text | Google Scholar

Hlaváčková-Schindler, K., Paluš, M., Vejmelka, M., and Bhattacharya, J. (2007). Causality detection based on information-theoretic approaches in time series analysis. Phys. Rep. 441, 1–46. doi:10.1016/j.physrep.2006.12.004

CrossRef Full Text | Google Scholar

King, F. W. (2009). Hilbert Transforms. encyclopedia of mathematics and its applications, 2. Cambridge University Press, 2. doi:10.1017/CBO9780511735271

CrossRef Full Text | Google Scholar

Krakovská, A., Jakubík, J., Chvosteková, M., Coufal, D., Jajcay, N., and Paluš, M. (2018). Comparison of six methods for the detection of causality in a bivariate time series. Phys. Rev. E 97, 042207. doi:10.1103/PhysRevE.97.042207

PubMed Abstract | CrossRef Full Text | Google Scholar

Kralemann, B., Cimponeriu, L., Rosenblum, M., Pikovsky, A., and Mrowka, R. (2007). Uncovering interaction of coupled oscillators from data. Phys. Rev. E 76, 055201. doi:10.1103/PhysRevE.76.055201

PubMed Abstract | CrossRef Full Text | Google Scholar

Kralemann, B., Cimponeriu, L., Rosenblum, M., Pikovsky, A., and Mrowka, R. (2008). Phase dynamics of coupled oscillators reconstructed from data. Phys. Rev. E 77, 066205. doi:10.1103/PhysRevE.77.066205

PubMed Abstract | CrossRef Full Text | Google Scholar

Kralemann, B., Frühwirth, M., Pikovsky, A., Rosenblum, M., Kenner, T., Schaefer, J., et al. (2013a). In vivo cardiac phase response curve elucidates human respiratory heart rate variability. Nat. Commun. 4, 2418. doi:10.1038/ncomms3418

PubMed Abstract | CrossRef Full Text | Google Scholar

Kralemann, B., Pikovsky, A., and Rosenblum, M. (2011). Reconstructing phase dynamics of oscillator networks. Chaos 21, 025104. doi:10.1063/1.3597647

PubMed Abstract | CrossRef Full Text | Google Scholar

Kralemann, B., Pikovsky, A., and Rosenblum, M. (2013b). Detecting triplet locking by triplet synchronization indices. Phys. Rev. E 87, 052904. doi:10.1103/PhysRevE.87.052904

PubMed Abstract | CrossRef Full Text | Google Scholar

Kralemann, B., Pikovsky, A., and Rosenblum, M. (2014). Reconstructing effective phase connectivity of oscillator networks from observations. New J. Phys. 16, 085013. doi:10.1088/1367-2630/16/8/085013

CrossRef Full Text | Google Scholar

Kreuz, T. (2011). Measures of neuronal signal synchrony. Scholarpedia 6, 11922. doi:10.4249/scholarpedia.11922

CrossRef Full Text | Google Scholar

Kreuz, T., Haas, J., Morelli, A., Abarbanel, H., and Politi, A. (2007). Measuring spike train synchrony. J. Neurosci. Methods 165, 151–161. doi:10.1016/j.jneumeth.2007.05.031

PubMed Abstract | CrossRef Full Text | Google Scholar

Kumar, M., and Rosenblum, M. (2021). Two mechanisms of remote synchronization in a chain of stuart-landau oscillators. Phys. Rev. E 104, 054202. doi:10.1103/PhysRevE.104.054202

PubMed Abstract | CrossRef Full Text | Google Scholar

Lehnertz, K. (2011). Assessing directed interactions from neurophysiological signals — an overview. Physiol. Meas. 32, 1715–1724. doi:10.1088/0967-3334/32/11/R01

PubMed Abstract | CrossRef Full Text | Google Scholar

Matsuki, A., Kori, H., and Kobayashi, R. (2023). An extended Hilbert transform method for reconstructing the phase from an oscillatory signal. Sci. Rep. 13, 3535. doi:10.1038/s41598-023-30405-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Moraffah, R., Sheth, P., Karami, M., Bhattacharya, A., Wang, Q., Tahir, A., et al. (2021). Causal inference for time series analysis: problems, methods and evaluation. Knowl. Inf. Syst. 63, 3041–3085. doi:10.1007/s10115-021-01621-0

CrossRef Full Text | Google Scholar

Mormann, F., Lehnertz, K., David, P., and Elger, C. E. (2000). Mean phase coherence as a measure for phase synchronization and its application to the EEG of epilepsy patients. Phys. D. 144, 358–369. doi:10.1016/s0167-2789(00)00087-7

CrossRef Full Text | Google Scholar

Musizza, B., Stefanovska, A., McClintock, P., Palus, M., Petrovcic, J., Ribaric, S., et al. (2007). Interactions between cardiac, respiratory and EEG-delta oscillations in rats during anaesthesia. J. Physiol. 580 (Pt 1), 315–326. doi:10.1113/jphysiol.2006.126748

PubMed Abstract | CrossRef Full Text | Google Scholar

Nakao, H. (2016). Phase reduction approach to synchronisation of nonlinear oscillators. Contemp. Phys. 57, 188–214. doi:10.1080/00107514.2015.1094987

CrossRef Full Text | Google Scholar

Nijholt, E., Ocampo-Espindola, J., Eroglu, D., Kiss, I., and Pereira, T. (2022). Emergent hypernetworks in weakly coupled oscillators. Nat. Commun. 3, 4849. doi:10.1038/s41467-022-32282-4

CrossRef Full Text | Google Scholar

Novaes, M., dos Santos, E. R., and Pereira, T. (2021). Recovering sparse networks: basis adaptation and stability under extensions. Phys. D. Nonlinear Phenom. 424, 132895. doi:10.1016/j.physd.2021.132895

CrossRef Full Text | Google Scholar

Paluš, M., and Stefanovska, A. (2003). Direction of coupling from phases of interacting oscillators: an information-theoretic approach. Phys. Rev. E 67, 055201. doi:10.1103/PhysRevE.67.055201

CrossRef Full Text | Google Scholar

Panaggio, M. J., Ciocanel, M.-V., Lazarus, L., Topaz, C. M., and Xu, B. (2019). Model reconstruction from temporal data for coupled oscillator networks. Chaos: An Interdiscip. J. Nonlinear Sci. 29, 103116. doi:10.1063/1.5120784

PubMed Abstract | CrossRef Full Text | Google Scholar

Pikovsky, A. (2018). Reconstruction of a random phase dynamics network from observations. Phys. Lett. A 382, 147–152. doi:10.1016/j.physleta.2017.11.012

CrossRef Full Text | Google Scholar

Pikovsky, A., and Rosenblum, M. (2022). “Non-pairwise interaction in oscillatory ensembles: from theory to data analysis,” in Higher-order systems. Editors F. Battiston,, and G. Petri (Springer), 181–195.

CrossRef Full Text | Google Scholar

Pikovsky, A., Rosenblum, M., and Kurths, J. (2001). Synchronization. A universal concept in nonlinear sciences. Cambridge: Cambridge University Press.

Google Scholar

Rings, T., Bröhl, T., and Lehnertz, K. (2022). Network structure from a characterization of interactions in complex systems. Sci. Rep. 12, 11742. doi:10.1038/s41598-022-14397-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Rings, T., and Lehnertz, K. (2016). Distinguishing between direct and indirect directional couplings in large oscillator networks: partial or non-partial phase analyses? Chaos 26, 093106. doi:10.1063/1.4962295

PubMed Abstract | CrossRef Full Text | Google Scholar

Rodriguez, E., George, N., Lachaux, J.-P., Martinerie, J., Renault, B., and Varela, F. J. (1999). Perception’s shadow: long distance synchronization of human brain activity. Nature 397, 430–433. doi:10.1038/17120

PubMed Abstract | CrossRef Full Text | Google Scholar

Rosenblum, M. G., and Pikovsky, A. S. (2001). Detecting direction of coupling in interacting oscillators. Phys. Rev. E 64, 045202. doi:10.1103/PhysRevE.64.045202

PubMed Abstract | CrossRef Full Text | Google Scholar

Rosenblum, M. G., Pikovsky, A. S., Kurths, J., Schäfer, C., and Tass, P. A. (2001). “Phase synchronization: from theory to data analysis,” in Neuro-informatics and neural modeling of handbook of biological physics. Editors F. Moss,, and S. Gielen (Elsevier), 4, 279–321.

CrossRef Full Text | Google Scholar

Runge, J., Gerhardus, A., Varando, G., Eyring, V., and Camps-Valls, G. (2023). Causal inference for time series. Nat. Rev. Earth Environ. 4, 487–505. doi:10.1038/s43017-023-00431-y

CrossRef Full Text | Google Scholar

Satuvuori, E., Mulansky, M., Bozanic, N., Malvestio, I., Zeldenrust, F., Lenk, K., et al. (2017). Measures of spike train synchrony for data with multiple time scales. J. Neurosci. Methods 287, 25–38. doi:10.1016/j.jneumeth.2017.05.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Shojaie, A., and Fox, E. B. (2022). Granger causality: a review and recent advances. Annu. Rev. Statistics Its Appl. 9, 289–319. doi:10.1146/annurev-statistics-040120-010930

CrossRef Full Text | Google Scholar

Siggiridou, E., Koutlis, C., Tsimpiris, A., and Kugiumtzis, D. (2019). Evaluation of Granger causality measures for constructing networks from multivariate time series. Entropy 21, 1080. doi:10.3390/e21111080

CrossRef Full Text | Google Scholar

Smirnov, D., Schelter, B., Winterhalder, M., and Timmer, J. (2007). Revealing direction of coupling between neuronal oscillators from time series: phase dynamics modeling versus partial directed coherence. Chao: An Interdiscip. J. Nonlinear Sci. 17, 013111. doi:10.1063/1.2430639

PubMed Abstract | CrossRef Full Text | Google Scholar

Smirnov, D. A., and Bezruchko, B. P. (2003). Estimation of interaction strength and direction from short and noisy time series. Phys. Rev. E. 68, 046209. doi:10.1103/PhysRevE.68.046209

PubMed Abstract | CrossRef Full Text | Google Scholar

Stankovski, T., Duggento, A., McClintock, P. V. E., and Stefanovska, A. (2012). Inference of time-evolving coupled dynamical systems in the presence of noise. Phys. Rev. Lett. 109, 024101. doi:10.1103/PhysRevLett.109.024101

PubMed Abstract | CrossRef Full Text | Google Scholar

Ticcinelli, V., Stankovski, T., Iatsenko, D., Bernjak, A., Bradbury, A. E., Gallagher, A. R., et al. (2017). Coherence and coupling functions reveal microvascular impairment in treated hypertension. Front. Physiology 8, 749. doi:10.3389/fphys.2017.00749

PubMed Abstract | CrossRef Full Text | Google Scholar

Tokuda, I. T., Levnajic, Z., and Ishimura, K. (2019). A practical method for estimating coupling functions in complex dynamical systems. Philosophical Trans. R. Soc. A 377, 20190015. doi:10.1098/rsta.2019.0015

PubMed Abstract | CrossRef Full Text | Google Scholar

Vlachos, I., Kugiumtzis, D., and Paluš, M. (2022). Phase-based causality analysis with partial mutual information from mixed embedding. Chaos: An Interdiscip. J. Nonlinear Sci. 32, 053111. doi:10.1063/5.0087910

CrossRef Full Text | Google Scholar

Xiong, W., Faes, L., and Ivanov, P. C. (2017). Entropy measures, entropy estimators, and their performance in quantifying complex dynamics: effects of artifacts, nonstationarity, and long-range correlations. Phys. Rev. E 95, 062114. doi:10.1103/PhysRevE.95.062114

PubMed Abstract | CrossRef Full Text | Google Scholar