Coupling functions in networks of oscillators

Networks of interacting oscillators abound in nature, and one of the prevailing challenges in science is how to characterize and reconstruct them from measured data. We present a method of reconstruction based on dynamical Bayesian inference that is capable of detecting the effective phase connectivity within networks of time-evolving coupled phase oscillators subject to noise. It not only reconstructs pairwise, but also encompasses couplings of higher degree, including triplets and quadruplets of interacting oscillators. Thus inference of a multivariate network enables one to reconstruct the coupling functions that specify possible causal interactions, together with the functional mechanisms that underlie them. The characteristic features of the method are illustrated by the analysis of a numerically generated example: a network of noisy phase oscillators with time-dependent coupling parameters. To demonstrate its potential, the method is also applied to neuronal coupling functions from single- and multi-channel electroencephalograph recordings. The cross-frequency δ, α to α coupling function, and the θ, α, γ to γ triplet are computed, and their coupling strengths, forms of coupling function, and predominant coupling components, are analysed. The results demonstrate the applicability of the method to multivariate networks of oscillators, quite generally.


Introduction
The networks found in nature [1] range from large-scale climatic interactions [2], through medium-scale synchronously-firing ensembles of neurons in the brain [3], to small-scale coupled molecular systems [4]. The nodes of the networks and their links (edges) can be static, or they may consist of dynamical systems and processes. There can also be a clustered subnetwork exerting a common influence over the rest of the network. The measured data for such networks may involve global observables, e.g. mean fields, or individual measurements of the dynamical nodes and their connections. The main challenge to be faced is how to characterize and reconstruct networks based on these kinds of data.
In what follows, we will focus on networks of oscillatory dynamical systems, a widespread class of networks that is particularly important for physiological processes like the brain-cardiovascular interactions [5][6][7][8]. Numerous methods exist for the inference of network couplings between oscillators [9][10][11][12][13][14][15][16][17][18][19], and it has also been shown that different forms of coupling can coexist between interacting systems [20]. From the statistics of the coupled signals, e.g. correlation and (bi-) coherence measures, one can find the functional connectivity [21] but such methods provide no information about causality or about the form of the coupling functions. We show, however, that cross-frequency coupling functions and their associated causality can be inferred from real data, thus yielding the effective connectivity [21]. Our approach will be based on a coupled-phase-oscillator model [9,22] and utilizes the recently proposed method of dynamical Bayesian inference [23][24][25][26][27]. Building on earlier work in this area [28][29][30], we will extend the method to encompass the inference of the coupling functions that prescribe the nature of the links (edges) between the oscillating nodes of a network.
In section 2 we develop the method itself. We start with a physically motivated introduction to coupling functions, describe how the interactions in a network can be decomposed into separate couplings, and outline the way in which Bayesian inference can be used in the analysis to reconstruct both pairwise connections and connections of higher order, even in the case of noisy, time-variable, dynamics. We also define some quantitative measures (coupling strength, similarity of form, and main coupling component) that can be used to characterize the coupling relations. Examples of the application of the method are presented in section 3 considering, first, data from a numerically generated system where all the answers are known in advance. Applications to brain dynamics data in the form of single-channel and multi-channel electroencephalograph (EEG) recordings are then described. Finally, we discuss the results and draw conclusions in section 4.

Coupling functions
Before going into details of the method and the analysis, we first focus our attention on the meaning and interpretation of coupling functions. In doing so, we use an elementary example to present the basic physics underlying a coupling function. We considered two oscillators that are unidirectionally phase-coupled: The goal is to describe the influence of the coupling function ϕ ϕ q ( , ) 2 2 through which the first oscillator affects the second one. From the expression for φ 2 in equation (1) one can appreciate the fundamental role of the coupling function: ϕ ϕ q ( , 2 is added to the frequency ω 2 . Thus changes in the magnitude of ϕ ϕ q ( , ) 2 1 2 will contribute to the overall change of the frequency of the second oscillator. Hence, depending on the value of ϕ ϕ q ( , ) 2 1 2 , the second oscillator will either accelerate or decelerate. This description of the coupling function is illustrated schematically in figure 1. Because in real situations one measures the amplitude state signals, we explain how the amplitude signals (figures 1(a) and (d)) are affected depending on the specific phase coupling function (figures 1(b) and (c)). In all plots, time is scaled relative to the period T 1 of the amplitude signal originating from the first oscillator x t ( ) 1 (e.g. presented on a π π × 2 2 grid (figure 1(b)) resembles a shifted cosine wave, which changes only along the ϕ 1 -axis: this kind of representation is used extensively in the discussion that follows below. Because all the changes occur along the ϕ 1 -axis, and for easier comparison, we present in figure 1(c) a ϕ 2 -averaged projection of ϕ ϕ q ( , )  decreases, x t ( ) 2 decelerates. Of course, the coupling functions of real systems can be much more complex than the simple example presented ( ϕ π + cos ( 2.5) 1 ), but the basic principle and physical significance of the coupling function are the same. A similar description can be elaborated for amplitude state coupling functions (e.g. q x x ( , 2 )-with the difference that the change in coupling function values will increase or decrease the amplitude of the driven oscillator.
Furthermore, the coupling functions determine the possibility of qualitative transitions between the oscillations, e.g. routes into and out of phase synchronization. This has important implications for many physiological interactions, including e.g. the phase-transitions of the cardiorespiratory interaction [31,32]. The coupling functions decomposition can enable a description of the functional contribution from each subsystem within a coupling relationship. It is this that has led to much recent progress towards the extraction and reconstruction of coupling functions between interacting oscillatory processes, with applications to cardiorespiratory interactions [23,33,34], chemistry [35][36][37], and communications [38]. Thus, the coupling function amounts to much more than just a new way of investigating correlations: it opens up a whole new perspective on the mechanisms underlying the functionality of a network.

Coupling decomposition
The properties of a network of N coupled periodic oscillators subject to noise can be investigated by focusing on its phase dynamics [22,39,40]. To do so, a system of N stochastic differential equations with time varying parameters can be built as where the frequency of each oscillator φ i is defined as the sum of its natural frequency ω i , a function q i of the phases ϕ … N 1, , of the oscillators influencing it, including its self-dynamics, and a stochastic part represented by ξ i , which is modelled as Gaussian white noise such that ξ ξ τ In this configuration it is easy to appreciate how the natural frequency ω i of each oscillator is deterministically influenced in time by the additive coupling function q i .
In the case of such an N-network model, the direct coupling from one oscillator to another, discussed in figure 1, is just one of the several possible combinations of couplings. In equation (3) the deterministic part q i of equation (2) is decomposed into the sum of the partial contributions of different orders of coupling [33,41], i.e. ′ q i represents the coupling from one oscillator, ″ q i the coupling from two oscillators, and so on  -Let us call the case in which the source and driven oscillator of the coupling are the same self-coupling (figure 2(a)). The self-coupling ϕ q ( ) A A has no physical relevance [42], and has the same morphology as that presented in figure 1, but the sinusoidal wave propagates along the ϕ A dimension as no other phase is involved.
-If the driven oscillator is different from the source, we refer to the coupling as direct (figure 2(b)): this is exactly the case discussed in figure 1, and ϕ q ( ) B A consists of a wave propagating along ϕ A while ϕ B remains constant [33,41].
• When the coupling comes from two oscillators, the form of the coupling shows a wave running in the diagonal direction (figures 2(c) and (d) at the bottom), generated from the presence of a sinusoidal function in both the dimension plotted. One can distinguish the cases when the driven oscillator is included among the sources, or not: in which the driven oscillator is involved in the coupling source common from two (figure 2(c)).
-The case ϕ ϕ q ( , ) C A B that does not include the driven oscillator in the coupling dynamics is referred to as a direct from two coupling (figure 2(d)).
• Even if higher orders of coupling cannot be plotted because they have more than two dimension in their domain, the nomenclature for classification can intuitively be extended, so that we call ϕ ϕ ϕ q ( , , ) common from three (figure 2(e)), and ϕ ϕ ϕ q ( , , ) C direct from three (figure 2 (f)), and so on. • At all levels of coupling, the description common is applied if and only if the driven oscillator has a selfcoupling.
• As notation in what follows, a single-source coupling from ϕ 1 to ϕ X will be indicated by ϕ ϕ → X

1
. For the overall coupling from ϕ ϕ … , , N 1 to ϕ X , the notation will be ϕ ϕ ; and to indicate the coupling components originating from a specific subset of sources among the overall coupling, the notation The deterministic periodic part of the differential equation (2) can be decomposed by Fourier approximation for each oscillator into a sum of base functions Decomposing the sum in equation (4) in the same way as in equation (3), in order to isolate the different orders of the network coupling, one obtains: elements), and so on. Therefore, the overall number ν of elements for the c matrix characterizing the Λ-order coupled network is given by

Dynamical Bayesian inference
The network model of N coupled phase oscillators, equation (5), is to be inferred by a dynamical Bayesian approach [23,27].
If the sampling frequency is high enough, i.e. the sampling step h is small enough relative to the dynamics, the phase dynamics described by equation (2) can be well-approximated from the time series  using the Euler midpoint discretization ϕ . Because the noise is treated as white and Gaussian, i.e. statistically independent, the likelihood at each time is considered as a product over l of the probability of observing ϕ + i l , 1 . The likelihood function is computed through the stochastic integral of the noise term over time, as where H is the Cholesky decomposition of the noise matrix D, andz i is a vector of normally distributed random variables. The joint probability density of the phase dynamics process in respect of ϕ ϕ − is calculated using the joint probability density of z i by imposing , where ξ ϕ J is the Jacobian term of the transformation of variables that can be calculated from the base functions Φ i k , . From this, the minus log-likelihood function where summation over the repeated indices k is implicit, and the dot index in ϕ · is substituted with the relevant index.
Assuming that the prior probability of parameters  is a multivariate normal distribution, and taking into account the quadratic form of the log-likelihood (7), the posterior probability will also be a multivariate normal distribution. With this particular distribution for the parameters c, with mean c, and covariance matrix Σ Ξ ≡ − prior prior 1 , the stationary point of S is calculated recursively [27] using the equations This inference technique is applied to the information provided by a stream of sequential blocks coming from the time-series. The current distribution (8) is computed, based on the evaluation of the previous block of data, i.e. informative priors are used. At each iteration information is propagated between the data windows, and the current prior depends on the previous posterior; because the first initial prior cannot contain any information, it is set to a flat normal distribution with Ξ = 0 prior and = c 0 prior . To handle networks characterized by time-variable interacting dynamics, explicit information propagation is used between consecutive blocks of data [23]: the covariance matrix of the next prior is computed by convolution of the current posterior with the current diffusion matrix Σ l diff . This measure describes how much the parameters can change:

Quantitative measures
One can utilize the inferred c to quantify certain characteristics of the coupling relations. They can be used either as quantitative measures, or to compare different coupling mechanisms.

Coupling strength
The coupling strength quantifies the coupling amplitude. It is defined as the Euclidean norm of the inferred parameters corresponding to the Fourier components of the coupling to the oscillator ϕ i from the combination of oscillators σ: part of the vector, and it is composed of × K 2 2 2 elements; the strength of each coupling from three oscillators is indexed into the ‴ c i ( ) part of the vector, and is composed by × K 2 3 3 elements, and so on.

Similarity of form
By calculating the correlation of two coupling functions, one can calculate the similarity of their forms, irrespective of their amplitudes [34]. In this way one can quantify the form of the coupling function-which is their unique characteristic. The similarity index is defined as where〈 × 〉 denotes averaging over the π π × 2 2 phase grid, and q i are the standard deviations = − 〈 〉 q q q i i i . The correlation ρ can quantify the similarity between two inferred coupling functions, or quantify how much one inferred coupling function is similar in functional form to some predefined characteristic coupling function.

Main coupling component
The concept of correlation has been extended to extract information about the predominant functional source of the coupling. The correlation with analytically predefined waves qp, can give a measure of the prevalence of a particular wave-propagating component within a coupling function. To take account of the possibility of the main wave being phase-shifted into a coupling function lying within the interval π − 0 2 , we compute the correlation between the coupling function and a series of M waves phase-shifted by an increment π m M 2 . To detect the similarity to narrower or broader waves, we repeat the procedure for waves with natural frequencies of ω × n . The maximum of these correlation values is chosen for each directional comparison: For the coupling from two oscillators (e.g. Figure 2(b)) we phase-shift bi-variate waves and check both the diagonals.

Numerical example
In order to illustrate and validate the inference procedure, we now apply the method to a numerically simulated system. The test network presented in this section is composed of five phase oscillators, coupled as illustrated in figure 3(a). The first oscillator X 1 oscillates without being influenced by any other, and with a time-varying natural frequency (dashed circle), X 2 is driven by the direct from two coupling + → X X X 1 3 2 (light green link), X 3 is driven by the direct from three coupling + + → X X X X 4 (light blue link), andX 5 is driven by the common from two coupling + → X X X 3 5 5 (dark blue link). The system of stochastic differential equations associated with the network is: where ξ i represents addictive white Gaussian noise, the length of the time-series is 2000 s and the sampling step is set as h = 0.01 (see [43] for details of the computation schemes). All the coupling coefficients σ c i: and the natural frequency ω 1 are built as time-varying and are defined by: The other frequencies are set to be constant, i.e. ω = 11.5 2 , ω = 21.5 3 , ω = 28 4 , ω = 42 5 . The time window for the dynamical Bayesian inference is set equal to 50 s, which has been chosen as a compromise between the need to include in each window enough information about the dynamics of the system (i.e. each 50 s window includes 5000 samples) and the resolution needed to follow the time variability of the parameters (i.e. 40 windows of 50 s).
In order to validate the results, a set of 100 surrogates of the original network was generated by shifting each phase time-series ϕ i by a random increment. This technique allows us to destroy any phase-to-phase correlation within the system, while still preserving the statistical properties of the time-series [44]. Figure 3(b), on the left, shows the coupling strength of the links effectively present in the net (green and blue bars), compared with the corresponding average coupling strengths plus two standard deviations (+2SD) computed from the set of surrogates (overlapped grey bars): it can be seen that the surrogates only attain values that are significantly lower, showing that the coupling links are correctly identified among the possible combinations checked. Figure 3(b), on the right, shows-on a different scale-the coupling strength inferred for some links that do not exist in the net (purple bars), plotted against the corresponding average +2SD strengths detected from the set of surrogates (overlapped grey bars): it can be seen that, in this case, the surrogates reach significantly higher values than the calculated couplings, implying that a coupling is (correctly) inferred to be below the significance level in cases when it is actually absent. Figure 3(c) presents a time-plot of the actual time-varying parameters (thin solid lines), as compared with their values inferred from the time-series (dashed lines). It can be seen that the method is not only able to detect the average value of a coupling (as indicated by the bars in figure 3(b)), but is also able to follow the time evolution of the parameters.

Application to brain dynamics
The ability of the method presented to infer time-varying and noisy dynamics makes it particularly suitable for application to EEG time series. We demonstrate this through an extension of the recently-introduced neuronal cross-frequency coupling functions [41]. In order to test the method on real data, two examples with quite different features have been analysed. First we analyzed a single frontal-probe recording from a subject under general anaesthesia. The second set of data was a 19-channel EEG recording from a young subject with autistic spectrum disorder (ASD), which allowed us to apply the methods to a distinctively different case. In performing the analysis, we concentrate in each case on the physical mechanisms linking the coupled systems; the physiological and clinical implications go beyond the scope of the present paper.

Preprocessing
To treat the system as a network of coupled phase oscillators the cognitive bands, i.e. the δ (0.8-4 Hz), θ (4-7.5 Hz), α (7.5-14 Hz), β (14-22 Hz) and γ (22-80 Hz) intervals, were filtered out from the signals by application of finite-impulsive-response and Butterworth filters. The phase was then extracted from each filtered time series using the Hilbert transform [45]. During this preprocessing procedure, particular care must be taken to minimise overlap between the filtered spectra [46]: overlaps of consecutive frequency intervals would result in overestimation of the corresponding phase-to-phase coupling.

Coupling relations of interest
First, we will analyse and discuss the δ, α α → coupling function. It has been found that the δ-waves, which typify deep sleep in adults [47] and still appear in the waking EEG of young subjects [48], can influence α-activity, which is related to the processing of information [47,49].
Secondly, we will study couplings in the γ-interval, which is associated with attention, memory and sensory processing [50]. The θ γ → and α γ → couplings have been already investigated pairwise, separately: the former coupling was related to memory tasks, while the latter was being found to be predominant during attentionrelated activities [11,14,51]. We will also apply our network approach to elucidate the dynamics of the θ α γ γ → , , triplet coupling. Figure 4(a) shows the δ α α → , coupling functions from the subject with a single EEG channel recording. In calculating the coupling function, all the partial contributions involved in δ α α → ,

Application to a single-channel EEG
were considered, namely the α α → self-coupling, the direct coupling δ α → , and the common coupling δ α α + → . The shape of the function is similar to the form presented in figure 2(b), implying that the coupling is dominated by the direct component. In particular, the predominance of tridimensional waves propagating along the ϕ δ dimension, together with the almost constant level along the ϕ α -axis, reveals that it is the δ that leads the α oscillations. Moreover, the position of the ridge and the hollow in the coupling function implies that the δ oscillations are accelerating the α when their phase is between π and π 2 , and decelerating them when they are between 0 and π (see figure 1 for details). Figure 4(b) shows the coupling that θ and α exerts on γ oscillations. In this case, the partial contributions included in the form are the direct couplings θ γ → and α γ → , and the direct coupling from the joint pair α θ γ + → . The shape of the coupling function indicates that the predominant role is played by the direct influence of θ over the γ oscillations.

Application to multi-channel EEG
The method was also applied to a set of multichannel EEG measurements. The recordings were derived from 19 probes distributed over the head of the subject according to the standard 10-20 configuration [52] illustrated in figure 5(a). The two channels defined by the 10-20 protocol as references, positioned on the earlobes, are not included in this analysis. The coupling strength was calculated for δ α α → , and θ α γ γ → , , . Unlike the form of the coupling function, the coupling strength representation does not take into account all the arguments of the coupling components. Hence the partial components can either be separately visualized or summed together. Figure 5(b) shows the strength of the δ α α → , coupling, which includes the self-coupling α α → , the direct coupling δ α → , and the common coupling δ α α + → . It is evident that the coupling strength is not equally distributed among the probes. Higher δ α α → , coupling is evident on the left hemisphere, with the frontal regions having highest coupling values. Figure 5(c) shows the strength of the θ α γ γ → , , coupling, which includes the self-coupling γ γ → , the direct coupling α γ → and θ γ → , the common coupling α γ γ + → and θ γ γ + → , the direct coupling α θ γ + → and finally the common coupling θ α γ γ + + → . The θ α γ γ → , , triplet coupling strength is mostly distributed among the outer probes. Figure 6(a) shows how the form of the δ α α → , coupling functions varies in relation to their spatial locations on the head. It can be seen in this figure that the tridimensional waves propagate mostly in the δ dimension. This tendency can be better observed in figure 6(b), which shows the averaged coupling function. Its form again depends predominantly on the direct δ oscillation, but it is now similar to two periods of a sine-wave, with two minima and two maxima along the ϕ δ -axis. This reveals that, within one δ cycle, there are two epochs of acceleration and deceleration of the α oscillations.
By computation of the correlation as described by equation (11), a numerical value can be obtained to quantify the visual evidence about the similarity of the coupling components. Figure 7 shows the values of the

Discussion and conclusions
The method presented allows one to study a network of coupled oscillators in greater depth than has been possible hitherto. Up to now, the concept of coupling had been limited mainly to the statistical properties of the time series in approaches that detected the existence of connections, and in some cases [9,10,[16][17][18] were also able to quantify coupling strengths. The present technique, on the other hand, allows one to discover the underlying mechanism of coupling as well, detecting the functional features of the interactions, i.e. from how many, and which, oscillators they are generated. The form of the coupling functions represents a new dimension in the characterization of a coupling. The underlying phase mechanism is also revealed, as the form of the coupling function determines how the intensity of the coupling depends on the relative phases of the two oscillators.
The method itself is based on the dynamical Bayesian framework [23][24][25][26][27], which infers well the timeevolving deterministic dynamics and separates out the noise perturbations. This makes the method especially suitable for inference of interactions among (oscillatory) biomedical systems that possess these characteristics. The method exploits a phase dynamics model and reveals the underlying causal relationships, which places it within the effective coupling class of methods [21]. It is mainly applicable to small-scale networks consisting of a few oscillators. Large numbers of oscillators will increase considerably the dimension and number of coefficients and base functions (equation (6)), thereby making the computations difficult in practice. In the case of networks consisting of a relatively large number of interacting oscillators, a compromise to obtain a reasonable number of coefficients can be set by adjustment of parameters such as the required order for the Fourier approximation (K in equation (6)) or the required complexity for the network interactions (Λ in equation (6)), and additional computing power (e.g. high performance clusters) might then be needed.
The network approach investigates all of the possible interactions simultaneously, correctly associating each coupling with the specific functional link that generated it. Figure 8 shows a comparison between the averaged δ α α → , coupling function computed (a) with the pairwise approach and (b) with the network approach. Clearly, the higher part of the function is computed consistently by the two methods, but the details at lower values appear to be slightly different. This demonstrates an important advantage of our network-based method over the pairwise-based approaches. That is, the approach proposed here and in [28,41,53] can distinguish the network-specific coupling contributions, i.e. triplets, quadruplets and higher-than-pairwise couplings in general, that are intrinsically present in the multivariate multidimensional interactions.
The feasibility of the method was demonstrated on a numerical example where the correct answers were known in advance. The method passed this test convincingly, and was shown to work with considerable precision. The specific application to brain dynamics is particularly promising, because of the intrinsic oscillatory nature of cognitive waves, and their underlying couplings. The ability of the method to decompose noise [23,27] makes it especially suitable for this kind of signal. The quantitative indices inferred allow one to detect differences in coupling among cognitive waves in different states, or having different functional structures, for example in pathological states. Inclusion of lower frequency oscillations in the network, e.g. those corresponding to cardiac and respiratory activity [54], could bring additional insight into the interactions between different physiological networks. The wide prevalence of networks of noisy oscillators makes the method applicable to a broad area of science and technology.