Reconstructing effective phase connectivity of oscillator networks from observations

We present a novel approach for recovery of the directional connectivity of a small oscillator network by means of the phase dynamics reconstruction from multivariate time series data. The main idea is to use a triplet analysis instead of the traditional pairwise one. Our technique reveals an effective phase connectivity which is generally not equivalent to a structural one. We demonstrate that by comparing the coupling functions from all possible triplets of oscillators, we are able to achieve in the reconstruction a good separation between existing and non-existing connections, and thus reliably reproduce the network structure.


Introduction
Understanding and mapping of the brain connectivity is one of the most challenging problems of neuroscience [1,2,3,4,5,6].The information provided by the matrix of interconnections of the human brain, or "connectome" [7], is essential for both basic and applied neurobiological research.It is believed that numerous severe disorders like Parkinson's and Huntington's diseases, autism, schizophrenia, and dyslexia are related to disorders in brain connectivity (see [1] and references therein).A particular important problem is to reveal the connectivity, and especially the directional connectivity, from multivariate data, e.g. from multichannel EEG or MEG measurements [8].Another, but related, problem is to find out whether two correlated nodes are connected directly or indirectly.In a more general context, the problem of revealing the connectivity can be formulated as a general data analysis problem for oscillator networks: how to recognize the network structure, i.e. the directions and the strengths of the couplings, from the observation.In particular, similar tasks arise also in the analysis of other physiological systems, e.g. in quantification of cardio-respiratory interaction [9,10,11] and in a number of other fields, e.g. in climate physics [12,13,14] and ecology [15].
Various techniques, mostly based on different measures of correlation [16,17], synchronization [18,19,20], and mutual or transfer information [21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36] have been exploited to tackle the connectivity problem.A separate group of methods relies on the coupled oscillators theory; these techniques assume that the nodes of the network are active, self-sustained, oscillators and reveal directional connectivity by reconstructing phase dynamics [37,38,39,20,40,41,42].In this paper we further develop the latter approach, by substituting the pair-wise analysis of a network by an analysis of triplets of nodes.We show that this essentially improves the reliability of the method.Although motivated by the problems of neuroscience, the technique can be exploited to recover the structure of networks of various nature.

Structural, functional, and effective phase connectivity
Before proceeding with the presentation of the technique we discuss different notions of connectivity.In neuroscience one typically differentiates between anatomical connections, functional connections, corresponding to correlations between pairs of nodes which may be not connected anatomically, and effective connections, representing direct or indirect causal influences of one node on the other [43,44].We illustrate this important issue by the following example.Suppose we have three self-sustained oscillators, coupled as shown in Fig. 1a.Here the arrows indicate physically existing, unior bidirectional, connections.For physical systems these connections are implemented by means of resistors, optical fibers, etc.For biological systems they correspond to anatomical connections, e.g.via synapses or via white matter tracts between pairs of brain regions [44].Mathematically, this means that, e.g., for the second node we can write ẋ2 = G 2 (x 2 ) + εH 2 (x 2 , x 1 ) , where the vectors ẋk are state variables of the kth oscillator, functions G and H describe the autonomous dynamics and the coupling, respectively; the coupling parameter ε explicitly quantifies the strength of interaction.Since the dynamics of the second oscillator depends on the state of the first one, we say that there is a structural connection 1 → 2. Obviously, the structural connectivity is directional.
Next, consider two physically uncoupled nodes 1 and 3. Since node 1 acts indirectly on the node 3 via the node 2 and both nodes 1, 3 have a common drive, they can synchronize or, generally, exhibit correlated oscillations.A measure of correlation or synchrony would reflect the non-directional functional connectivity (Fig. 1b).Synchrony is typically quantified by means of the synchronization index, also known as the phase locking value, see e.g., [45,19,46].Since it is computed from the phases of interacting units, see Eq. ( 6) below, it quantifies functional phase connectivity.Now we discuss the notion of effective phase connectivity, crucial for our approach.For this purpose we recall, that if the coupling between the nodes is weak than the dynamics of N coupled systems can be reduced to that of N phases.It means that each (c) If there is a directional structural link k → l, then phase dynamics of node l depends on the node k, i.e. the same link exists in the effective phase connectivity map.However, additional links may appear.So, for this example, phase dynamics of node 3 may depend on that of node 1, i.e. these nodes may be also effectively connected (dotted line).Notice that the phase connectivity is also directional.So, e.g., the 1 → 3 effective phase connectivity is possible in this configuration, while the 3 → 1 connectivity is not.
oscillator, which can be a high-dimensional dynamical system, can be characterized by a single variable, the phase ϕ k [47,48].This variable is introduced in such a way, that it grows uniformly in time, i.e. φk = ω k .
The weakness of coupling means that the motion in the state space takes place on the surface of the smooth invariant N -dimensional torus.The equations of the phase dynamics are: where the new phase coupling functions h k can be obtained from the state space coupling functions H 2 (x 1 , x 2 , . . ., x k ) by means of a perturbative reduction [47] in the form of the power series of the coupling parameter ε: ) If the coupling function h k depends on ϕ l , we say that there exist a directed phase connectivity link l → k.
The crucial issue is that even if the structural connections are pair-wise, i.e.
the phase coupling functions h k can have terms, depending on many phases, not only on the phases of directly coupled nodes [42].So, for the example illustrated in Fig. 1a, the equation of the phase dynamics for the third node has the form Notice, that the indirect coupling 1 → 2 → 3 is reflected by the term ∼ ε 2 ; hence, for very weak coupling the indirect link 1 → 3 is negligible and the effective phase connectivity coincides with the structural one.However, for not very small ε, the indirect coupling 1 → 3 becomes not small either.Notice, that the indirect coupling 3 → 2 → 1 is not possible, since there is no information flow from node 3 to node 2. Thus, generally, the structural and effective phase connectivities of a network differ (Fig. 1c), although for a small coupling they nearly coinside.Below we present a technique for reconstruction of the effective phase connectivity from multichannel data.

Effective connectivity from phase dynamics
Suppose the phase model Eq.(2) of the system of N interacting self-sustained oscillators is known.It is convenient to represent the right hand side of Eq. ( 2) as a Fourier series: The norm of the coupling function h k quantifies the influence of the rest of the network on the oscillator k.In order to quantify the action of a particular oscillator j, we introduce the partial norms F (k) 0,...,l k ,0,...,l j ,0,...
i.e. from the Fourier series we pick up only the terms depending on the phases ϕ k , ϕ j and terms, depending on ϕ j .Notice that in a similar way, picking up the terms, depending on three phases, one can quantify the joint action of two oscillators on the chosen one, m, j → k, etc, see [42].Thus, if the phase model can be reconstructed from data, we can compute partial norms N k←j for all combinations k, j and in this way obtain a description of the network connectivity.This description is quantitative and directional, since generally N k←j = N j←k .

General approach
The first step of the approach is to obtain phases of all oscillators.We assume that we measure the outputs of all nodes of the network and that these outputs are suitable for phase estimation.Hence, we start with N oscillatory time series and then compute N continuous-time protophases θ k (t) [38,39], using e.g., the Hilbert Transform or any other embedding.(The choice of the optimal technique for protophase estimation is beyond the framework of the current study.)We recall, that the true phase of any autonomous self-sustained oscillator shall satisfy Eq. ( 1).The angle variables (protophases) θ k provided by the Hilbert transform or other techniques generally do not have this property and therefore shall be transformed to true phases θ k (t) → ϕ(t), as discussed in the Appendix A.
Next step is to compute numerically the derivatives of all phases (here it is done simply by computing central final difference; for noisy data the application of the Savitzky-Golay filter is recommended) and to fit the coefficients of the model (2), as discussed in [42].
As discussed below, the reconstruction is not possible if oscillators synchronize.Practically, in order to exclude from the consideration the syncronous states, we compute pair-wise and triplet synchronization indices [49] according to and Notice that in Eq. ( 6) n, m are positive integers, while in Eq. ( 7) n, m, p can be both positive and negative.We exclude from the further analysis the states where γ j,k , γ j,k,l > 0.5.

Effect of network size
In the presented approach we face a numerical problem related to the system's dimension.Consider first N = 2, then the coupling functions h 1,2 are functions of two variables ϕ 1,2 and in order to recover them by means of fitting we need that the data points cover the square 0, 2π, 0, 2π in the ϕ 1 , ϕ 2 plane.If two oscillators synchronize (i.e. the index ( 6) is close to one), then ϕ 1 , ϕ 2 are functionally dependent, the points do not cover the square but fall on a line, and the function of two variables cannot be recovered.
Let now N = 3 and suppose that oscillators are far from synchrony.For successful data reconstruction we need enough points to cover the cube 0 < ϕ j ≤ 2π.This is feasible, however, the amount of data required to fill the hyper-cube in dimension N grows rapidly with the dimension, making the reconstruction of high-dimensional coupling functions practically impossible already for N > 3. Thus, in dimensions larger than three only partial phase dynamics reconstruction is possible.

Partial phase dynamics
The simplest and most typical approach is to perform a pair-wise analysis of the network, reconstructing the phase dynamics for all pairs of nodes.In this way one does not reconstruct the full Eq.( 2), but fits the pair-wise coupling functions h kj using where all time series except for ϕ k (t), ϕ j (t) are ignored.The norms P k←j = ||h kj || are then used for the quantification of the network connectivity.As we show in details in Section 5 below, this estimation may yield spurious effective phase connections.Indeed, suppose one has a network like in Fig. 1a and determines phase connectivity according to (8).Because ϕ 1 is correlated with ϕ 2 , and ϕ 2 drives ϕ 3 , consideration of the pair ϕ 1 , ϕ 3 will give non-zero coupling from ϕ 1 to ϕ 3 .Only full phase dynamics according to Eq. ( 2) would reveal absence of direct coupling between ϕ 1 and ϕ 3 and presence of indirect coupling 1 → 2 → 3. Making use of this observation, we suggest to perform partial triplet analysis for networks with more than three oscillators as well.In doing this, we reconstruct the threedimensional coupling functions for all possible triplet configurations.Namely, from the time series of three phases ϕ j , ϕ k , ϕ l (all indices j, k, l are different) we first reconstruct the coupling terms h jkl , h klj , h ljk , ignoring all other phases.From these functions we compute, according to Eq. ( 5), partial norms Tj←k (l) for all binary connections within this triplet (all these norms for a triplet correspond to different permutations of symbols j, k, l).Next, we repeat this procedure for all triplets within the network.As a result, we obtain N − 2 estimates of Tj←k for the connection j ← k, since the systems j and k belong to N − 2 different triplets.We suggest to take the minimal value of these estimates as the final triplet-based measure of the binary effective phase connectivity: To support this approach, let us assume that the motif in Fig. 1a is a part of a larger network.From the triplets including oscillators 1, 3, and some disconnected oscillator n one would obtain a strong binary term 3 ← 1.However, in the triplet {123} this term will be small, as here the analysis recognizes that in fact oscillator 3 is driven by oscillator 2, and not by oscillator 1.This triplet will correctly yield a small value for the link 3 ← 1.

Three oscillators
To demonstrate the method, we start with the case of three coupled units, where the triplet analysis gives a full description of the phase dynamics, and compare the results with the partial pair-wise analysis.We introduce our approach using the example of three coupled van der Pol oscillators: Here the coefficients σ ij are either zero or one; they determine the structure of the coupling network in the three-oscillator system, whereas the parameter ε = 0.2 determines the coupling strength.The other parameters are: µ = 0.5, 75483.The choice of the frequencies is motivated by an attempt to avoid synchronous states, where the reconstruction is not possible.We generate three-channel data by solving Eq. ( 10) with time step 0.05s and use 10 5 data points per channel for the phase dynamics reconstruction.The protophases θ 1,2,3 are computed via the Hilbert transform and then transformed to genuine phases ϕ 1,2,3 via the phase transformation, see Appendix A. Next, we estimate coupling functions, approximated by a Fourier series of the order 5, cf.Eq. ( 4), and compute partial norms N k←j according to Eq. (5).For comparison, we also perform pairwise analysis, i.e. we compute the coupling function for the pair k, j neglecting the third oscillator to obtain P k←j .
We consider two examples of three-oscillator networks, illustrated in Fig. 2. The corresponding results are given in Tables 1(a,b), which represent the reconstructed coupling matrices.First we discuss the results for the chain of oscillators (Fig. 2a), see Table 1(a).We emphasize several important issues.(i) The reconstructed measures of the coupling strength, i.e. the norms N k←j , P k←j , corresponding to the existing structural connections are almost two orders of magnitude larger than the norms for unconnected oscillators, with an exception for the link 1 ← 3, discussed next.(ii) While oscillator 3 is not connected to oscillator 1 directly, there exist an indirect causal link mediated by oscillator 2. Therefore the exceptionally large value N 1←3 = 0.018 captures some real flow of information.Contrary, all values of the norms which correspond neither to a direct nor to an indirect causal link, practically vanish.This fact illustrates that although the effective phase connectivity generally differs from the structural one, this difference is not due to artifacts, but reflects the properties of the network.(iii) As expected, the obtained norms are asymmetric, what shows that our analysis captures essentially more than simple correlation.(iv) While the results from the triplet and the pairwise analyses, i.e.N k←j and P k←j practically coincide for existing structural links, they differ for some absent structural connections.Here, the values from the triplet analysis are notably smaller than those from the pairwise one, cf.N 2←1 = 0.002 and P 2←1 = 0.009.This happens because the oscillator 3 is not included into the pairwise model and its effect on the correlation between oscillators 1 and 2 is erroneously assigned to the link from 2 to 1.This essential difference between the pairwise and the triple analysis will be illustrated in detail below.
In the second example (Fig. 2b) two oscillators are driven by the third one.The results in Table 1(b) corroborate those from the first example.Indeed, the reconstructed norms for existing connections are again essentially larger than those for the non-existing connections.Moreover, here the analysis yields a perfect coincidence of the effective phase connectivity and of the given structural one.Since in this configuration there are no indirect causal links, this fact substantiates the conjecture that differences between effective phase and structural connectivities are due to indirect coupling.The results from the triplet and the pairwise analysis are again identical in case of existing structural links and likewise differ when no direct structural coupling is present and the dynamics is not autonomous, as can be seen for links 1 → 3 and 3 ← 1.We see that here the pairwise analysis again mistakenly takes correlations due to the common drive for causal interaction.Hence, for this link the pairwise analysis fails, while the triplet analysis yields reasonable results.

Four oscillators
With the next example of a network (Fig. 3) we want to address the obvious problem: for N > 3 each oscillator enters more than one triplet configuration.So, e.g. the link 1 ← 2 can be quantified by the analysis of the triplet configurations {123} and {124}, and we have to choose between two results.To discuss this problem, we use the example with ω 1 = 1, ω 2 = 1.3247, ω 3 = 1.75483, and ω 4 = 1.5333.All other parameters of the system, of the numerical simulation, and of the data processing are the same as in the previous examples.The results are given in Table 2.The values in the last row are due to numerical errors; indeed, since oscillator 4 is autonomous, these values shall be zero.However, they are much less than the values corresponding to structural connections.The analysis of other entries supports our conjecture that the value of T k←j is overestimated, if an oscillator which drives the unit k is excluded from the triplet.The reason is that the correlations between the oscillators k, j due to an external (with respect to the considered triplet) drive are spuriously explained as a coupling between k, j.Consider, e.g., the values of T 1←4 : in the triplet configuration {134} which includes both units driving the first oscillator, T 1←4 ≈ 0.016, i.e. about two times smaller than in the configuration {124}.Similar observation can be done for the link 2 ← 4. For a comparison we also present the values obtained from the pair-wise analysis; for the absent links they are overestimated as well.
The results support the suggested strategy to analyze networks of oscillators by covering them with all possible triplets.In a network of N oscillators each partial norm T k←j is then computed from N − 2 triplets, and, hence, N − 2 different estimates of its value are available.Since improper triplets overestimate T k←j , we choose the minimal over all triplets value min{T k←j } for the estimate of the directional coupling between oscillators k, j.

Random oscillator networks
In this section we report on the statistical analysis of the quality of reconstruction of effective phase connectivity in random networks of N = 5 and N = 9 oscillators, using the method presented above.The equations of the model read For each run, the random frequencies ω k are taken from the uniform distribution 0.5 < ω < 1.5.The random asymmetric connection matrix σ kl is composed of zeros and ones, with a fixed number of two incoming connections for each oscillator for N = 5, and with four incoming connections for N = 9.Parameter Θ governing effective phase shift of the coupling is taken from a uniform distribution 0 ≤ Θ < 2π.The protophase estimates are obtained as φ k = arctan( ẋk /ω k x).From the phases we calculated the maximum over binary and triplet synchronization indices according to Eqs. (6,7) and analyzed only cases with the maximal index less than 0.5; for computation of the indices we used |n|, |m|, |p| ≤ 5.
The results for N = 5 and for three different values of ε are presented in Fig. 4.Here we show the pairwise coupling norms P vs the coupling norms T , based on the triplet analysis, for the same data.Norms corresponding to the existing couplings (σ kl = 1) are marked with green, while norms corresponding to non-existing connections (σ kl = 0) are shown with red; 10 6 points are used for the analysis.Summarizing these results, we conclude: (i) For the existing connections, the triplet analysis practically does not change the value of the partial norm, obtained from the pairwise analysis: all green points lie very close to the diagonal (dashed line), so that P ≈ T .(ii) For non-existing connections, the triplet analysis significantly reduces the values of the coupling norms.For small coupling strength ε = 0.02 this reduction is rather strong, with typical factor about 0.3.For larger coupling the reduction is less pronounced.(ii) The values of the coupling from the pair-wise analysis do not allow to distinguish between existing and non-existing connections unambiguously, as the distributions of the obtained norms overlap (panels (d,e,f), cf.similar observation in [42]).However, for the norms obtained from the triplet analysis, the existing and non-existing couplings can be clearly separated, at least for the small coupling strength ε = 0.02.The separation remains non-perfect for larger coupling strengths, but nevertheless it is definitively better than in case of the pair-wise analysis.
In our approach we relied on Eq. ( 3) and on the statement that the indirect connections are reflected by the Fourier terms of the coupling function, proportional to ε 2 , while the corresponding terms for the structurally existing links appear already in the first approximation in ε.This is confirmed by Fig. 5, where we present the re-scaled  (d,e,f) show distributions for the norms P (red for non-existing and blue for existing connections, also marked as ∓2) and T (green for non-existing and magenta for existing connections, also marked as ∓3).For existing connections the distributions for P, T practically coincide, while for non-existing connections they are essentially different: the distributions for the norms T obtained from the triplet analysis demonstrate very good separation between existing and non-existing connections.
distributions of norms T .The overlap of re-scaled distributions firmly supports that the effective indirect phase connectivity is not the artifact of the method, but a real coupling appearing in the second order in the coupling strength.
The results for nine-oscillator networks, are shown in Figs.6,7, for two values of the coupling strength; here 10 5 points are used for the analysis.The results are generally similar to those for five coupled oscillators.Notice the shift of maxima of curves A,C in Fig. 7 with respect to those of curves B,D; this reflects the fact that our approach always overestimates the norms.Since the coupling in this example is rather weak (because of a larger number of incoming links compared to five-oscillator networks), the norms proportional to ∼ ε 2 are of the order of numerical precision.However, as follows from Figs. 6c,d, the separation between existing and non-existing links is still very good.the performance of these methods is a serious computational and programming task, and requires a separate study; see, e.g., [50,51,52,53,54].Several algorithms are implemented in software packages [55,56,57], see also http://hermes.ctb.upm.es.

Discussion
Here we discuss briefly two other methods.The first one is based on Granger causality and is implemented as a user-friendly toolbox tailored to evaluate the causal connectivity underlying neural data [57].With this toolbox we estimated the coupling structure for the four coupled van der Pol oscillators, Eqs.(11).Comparison with our method, presented in details in Appendix B, shows that both methods basically correctly estimate structural links.On the one hand, our method provided more homogeneous results (all existing links have approximately the same level of coupling, while the Granger causality approach gives couplings differing up to 75%).On the other hand, the phase model based results indicate weak indirect coupling for not directly linked oscillators, while the values of the Granger connectivity in this cases correspond exactly to the structural coupling.This comparison shows that quality of our approach of revealing the connectivity of oscillatory networks is comparable to the Granger causality approach, with an advantage that its results are widely insensitive to the choice of parameters.
Next, we discuss a relation to the technique [58], based on computation of synchronization indices, cf.[59].Since the indices are a symmetric measure, this technique yields a non-directed measure of phase interdependence, i.e. it quantifies functional phase connectivity.The method also differentiates between existing and non-existing links.Reliable performance requires that synchronization indices are rather high, so that the system shall be close to synchrony, and is based on the statistics of the deviations from this synchrony (ideal synchrony yields identity matrix which contains no information).On the contrary, performance of our technique is better if the network is far from synchrony; thus, in this respect the methods are complimentary.However, the method [58] requires the 1 : 1 synchrony, i.e. closeness of the frequencies of all nodes; our technique does not have this limitation.

Conclusions
Discussing the pros and cons of our approach we emphasize, that it explicitly exploits the assumption that all nodes of the network are active, self-sustained oscillators, so that the phase dynamics description is meaningful.This is a disadvantage of the approach, compared to the information theory based methods which are free of any assumptions regarding the type of the dynamics.On the other hand, for this particular case our techniques yields an appropriate description which admits a clear physical interpretation.We clearly interpret the output of our algorithms and the recovered network structure, by theoretical argumentation and numerical demonstration, as revealing the effective phase connectivity which coincides with the structural one in case of very weak coupling and shows additional links when the coupling is not too weak.Noteworthy, the additional links are not artifacts: most of them quantify the indirect influence between nodes appearing in the higher orders of phase reduction, and correspond to a causal information flow.Although we cannot unambiguously separate the links which are physically meaningful from those which are due to computational errors, we suggest as a rule of thumb a practical approach based on the triplet analysis.Indeed, the computed norms are always positively biased and can never be zero.However, as shown in the examples of Fig. 2a and Fig. 3, see also [42], the norms corresponding to the real information flow are always at least several times larger than the noise level.Figure 5 demonstrates, that the width of the distribution for the norms of non-existing structural connections is about two orders of magnitude.Hence, we can reasonably assume that norms that are, say, five times larger than the minimal one, correspond to indirect connections.The discrimination can be improved using the information on the direct connections (which can be recovered reliably, as discussed above): if there is a path between nodes j, k via node m, then the indirect connection j → k is very likely to exist.
We stress that although the method is based on the phase description Eq. ( 2), it is not restricted to the case of very weak coupling.Indeed, analytical derivation of Eq. ( 2) requires weakness of coupling, but numerical reconstruction is possible as long as the signals can be considered as nearly periodic or weakly chaotic.Certainly, application of the method implies that the outputs of all nodes, suitable for phase estimation, are available.
An important issue is that Eq. ( 2) is valid also for transient processes.Hence, the techniques does not imply stationarity of the data.So, e.g., if the network is repeatedly stimulated (cf.[60]), the pieces of data between the stimuli can be used for the phase dynamics reconstruction and, therefore, characterization of the network connectivity.Moreover, repeated perturbation can be a useful tool in case when the network (or some nodes) are close to synchrony, so that the reconstruction without perturbation fails.Perturbing the system and observing its relaxation to the synchronous state one can obtain enough data for successful application of the technique.
As an issue for further extension of the technique we mention accounting of terms, depending on the phases of three nodes.These terms appear in the second approximation in ε even for a purely pair-wise coupling and can be of the order of ε, if the coupling is nonlinear [42].These terms describe the joint action of two oscillators on the third one.Quantification of this action is straightforward, but representation of results for a network of more than three oscillators is not trivial.
In summary, we have shown that the triplet analysis of an oscillator network yields an essential improvement in the quantification of the connectivity compared to the conventional pair-wise analysis.The described technique provides effective phase connectivity and allows one a reliable differentiation between structurally existing and non-existing links.

Appendix A. Protophase to phase transformation
The true phase ϕ of an autonomous oscillator is introduced according to Eq. ( 1), i.e. it grows uniformly with time.However, an angle variable, or a protophase θ, obtained from a scalar time series by means of a two-dimensional embedding, e.g.via the Hilbert transform, generally does not have this property and obeys θ = ω + g(θ) .
If the oscillator is coupled, e.g. to another one, then its phase dynamics is described by where the coupling function h -the goal of our analysis -is small while the function g is generally not, cf.Eq. (2).Hence, the reconstruction of the function h is hampered by the function g and the latter shall be eliminated by means of the transformation θ → ϕ(θ).For a noise-free system the transformation is easily obtained from Eq. (1) using the chain rule: For noisy systems we have to average dt dθ (t) over the trajectory to obtain the transformation function σ(θ): Up to a factor 2π, this function is the probability density of θ and, hence, can be easily obtained from the time series Θ(t k ), k = 1, . . ., M , see [38,39].The final transformation reads: where S n = M −1 M k=1 exp(−inΘ k ) are coefficients of the Fourier expansion of σ(θ).Practically, one has to truncate the Fourier series and compute some finite number K of terms in Eq. (A.1).An efficient algorithm for determination of the optimal K was suggested by C. Tenreiro [61].A Matlab code for the protophase-to-phase transformation with implementation of the Tenreiro optimization can be downloaded from www.stat.physik.uni-potsdam.de/∼mros/damoco.
We illustrate the importance of the transformation by the following example.We consider two uncoupled Hindmarsh-Rose neuronal oscillators: where I 1 = 5, I 2 = 5.1.For the chosen parameter values both systems exhibit periodic spiking.Next, we compute the protophases via Hilbert transform of x 1,2 .Finally, using Eq. ( 6) we compute the synchronization index twice, from the protophases and from the true phases.Since the systems are uncoupled, the true value is zero.The index, obtained from the protophases yields an obviously spurious value ρ ≈ 0.13, while the estimation using the true phases after the transformation above yields a reasonable result, ρ ≈ 0.02.

Appendix B. Comparison with Granger causality method
Here we describe details of causal connectivity evaluation, based on Granger causality concept, with the help of the toolbox [57].Namely, we reconstructed the coupling pattern of four coupled van der Pol oscillators, see Eqs. (11), using the numerically obtained solutions x k (t) as an input for the toolbox.We followed the steps described in the manual, i.e. the signals have been detrended, demeaned, and standardized.Next, we have checked the statistical properties of the time series, namely whether the data are covariance stationary.We found, that neither the data nor their derivatives meet this requirement due to pronounced periodicity, what means that the Granger causality approach is possibly not well-suited for our problem.Thus, the results should be interpreted carefully.Evaluation of the Granger connectivity from data requires selection of a time lag parameter, L. The algorithm for optimization of L did not converge for our data.Therefore, we tried different parameter values and found reasonable connectivity patterns for time lags L = 2 and L = 3.In both cases all values indicating causal influence where judged as significant results, see Table B1.
Although both results indicated a similar pattern of connectivity, the absolute values strongly varied with the time lag.Since the results for L = 2 are closer to the true structural connectivity of the system, we compared them to our results.To simplify the comparison of different magnitudes, we normalized the maximal values of the respective connectivity matrices to 1.The values corresponding to structural coupling of the oscillators are shown with bold font in Table B2.Ideally, this bold values should all be 1 while all others should be 0. The values corresponding to structural coupling obtained from the phase model (the entries marked by PM) have less than ∼ 10% deviance, while the results based on Granger causality (GC) differ up to 75%.On the other hand,

Figure 1 .
Figure 1.Illustration of the notions of connectivity.(a) Nodes 1, 2 and 2, 3 are physically (structurally) connected.This connection can be bi-or unidirectional, as shown by arrows.(b) Structural connection causes correlated (or synchronous)dynamics, i.e. functional connection (dotted lines).Since measures of correlation and synchronization are symmetric, functional connectivity is not directional.Notice that although nodes 1, 3 are not physically connected, they may be correlated due to the common drive from the node 2 and indirect action of node 1, what yields an additional link if compared with (a).(c) If there is a directional structural link k → l, then phase dynamics of node l depends on the node k, i.e. the same link exists in the effective phase connectivity map.However, additional links may appear.So, for this example, phase dynamics of node 3 may depend on that of node 1, i.e. these nodes may be also effectively connected (dotted line).Notice that the phase connectivity is also directional.So, e.g., the 1 → 3 effective phase connectivity is possible in this configuration, while the 3 → 1 connectivity is not.

Figure 2 .
Figure 2. Two examples of three-oscillator networks.

log 10 PFigure 4 .
Figure 4. Results for random five-oscillator networks, for ε = 0.02 (a,d), ε = 0.05 (b,e), and ε = 0.1 (c,f).In (a,b,c) the points, corresponding to existing and nonexisting connections are shown with green and red, respectively; dashed line is the identity line.Panels (d,e,f) show distributions for the norms P (red for non-existing and blue for existing connections, also marked as ∓2) and T (green for non-existing and magenta for existing connections, also marked as ∓3).For existing connections the distributions for P, T practically coincide, while for non-existing connections they are essentially different: the distributions for the norms T obtained from the triplet analysis demonstrate very good separation between existing and non-existing connections.

log 10 (Figure 5 .log 10 PFigure 6 .
Figure5.Scaling of coupling terms for random five-oscillator networks.Here the distributions of T from Fig.4d,e,f are re-scaled: the norms of existing links are divided by ε, whereas the norms of non-existing connections are divided by ε 2 .The result confirms our conjuncture that terms, describing indirect phase connectivity appear in the second order of the phase approximation, i.e. are ∼ ε 2 .Moreover, it confirms that most of the non-existing connections revealed by the technique are not artifacts, but reflect causal information flow, mediated by the indirect driving.Curves A,B correspond to the rescaled distributions marked as ∓3 in Fig.4d, while C,D and E,F are the corresponding re-scaled distributions from Fig.4eand Fig.4f, respectively.

Table 1 .
(a) Reconstruction of the coupling for the example in Fig.2a.An entry in the k-th row and j-th column represents N k←j , P k←j , which are measures of the coupling strength.With bold font we show existing structural connections.Hence, the corresponding entries should be much larger than the others; for the example under consideration this is obviously true.Notice that since the norms are by definition positive, we always obtain overestimated values for non-existing connections.(b)Sameas (a), but for the network configuration shown in Fig.2b.

Table 2 .
Reconstruction of the coupling for the example in Fig.3.An entry in the k-th row and j-th column represents two values for the partial norm, Tk←j , obtained from two different triplets which include nodes k, j; the corresponding triplets are denoted as {kjn}, the values obtained from the pair-wise analysis are denoted as {kj}.Existing structural connections are shown with bold font.The minimal coupling norms T k←j are shown in boxes.

Table B1 .
Results of application of the connectivity estimation based on Granger causality, for two time lags L = 2 and L = 3 (shown in {}).Osc 4 0.000 {2} 0.000 {2} 0.000 {2} 0.000 {3} 0.000 {3} 0.002 {3} the phase model based results indicate weak indirect coupling for not directly linked oscillators, while the values of the Granger connectivity in this cases correspond exactly to the structural coupling.This comparison shows that our approach is better suited for reconstructing the connectivity of oscillatory networks, because its results are widely insensitive to the choice of parameters.