We computed the wavelet transform (WT) (Kaiser, 1994; Bračič and Stefanovska, 1998; Stefanovska et al., 1999) in order to evaluate the power content within the 0.8–40 Hz range, converting the signals s(t) to the time-frequency domain: (1)WT(ω,t)=∫0∞ψ(ω(u-t))s(u)ωdu, where ω denotes angular frequency, t is time, and ψ(u) = 1/(2π) (e(i2πf0u)-e(2πf0)2/2)e-u2/2 (with ∫ ψ(t)dt = 0) with central frequency f₀ = 1. The power within each frequency interval was assessed by averaging the spectra over the corresponding frequency ranges.

Amplitude dynamics in living systems is often multidimensional, which can create complications in analysis. In contrast, the phase dynamics of a periodic process in such systems is describable in terms of a single-dimensional observable, which is usually much easier to detect and extract from data. It is well known that brain activity carries the signatures of several distinct neural oscillations that manifest themselves within characteristic frequency intervals (Buzsáki and Draguhn, 2004). The signals extracted from these intervals are periodic, enabling the underlying oscillatory processes and their interactions to be studied effectively through phase dynamics (Kuramoto, 1984), and leading to extraction of phase-to-phase cross-frequency couplings (Jensen and Colgin, 2007; Jirsa and Müller, 2013). The cross-frequency phase couplings coexist in a multivariate and multidimensional space, so we will consider a network model of N coupled phase oscillators, each described by (2)ϕ˙i(t)=ωi(t)+qi(ϕi,ϕj,ϕk,…,ϕN,t)+ξi(t)=ωi(t)+∑nqi(1)(ϕn,t)+∑nmqi(2)(ϕn,ϕm,t)+∑nmlqi(3)(ϕn,ϕm,ϕl,t)+…+ξi(t), for all l, m, n, …, where ϕ˙i(t) is the time derivative of the phase (i.e., the instantaneous frequency), ω_i(t) is the natural frequency and the external stochastic dynamics ξ_i(t) is treated as Gaussian white noise 〈ξ_i(t)ξ_j(τ)ξ_k(τ)…〉 = δ(t − τ)D_ijk…, where D is the matrix of noise diffusion and D_{i, j, k..} gives the noise strength for the particular i, j, k… element. Although we will discuss the inference of neural coupling functions from phase dynamics, the method that we will describe is in principle also applicable to their inference from amplitude dynamics (Stankovski et al., 2014b).

The coupling functions qi(κ) describe the dynamics in terms of the phases of κ interacting oscillators. As can be seen from Equation (2), the coupling functions q_i act in such a way as to modify the natural frequency ω_i(t): in physical terms, a positive coupling coefficient will accelerate the oscillation in question (by increasing its instantaneous frequency ϕ˙i(t)), whilst a negative coupling coefficient will decelerate it (by decreasing ϕ˙i(t)). Thus a coupling function is able to describe in detail, within a single cycle, how one oscillator is accelerated or decelerated as a result of the influence from the other oscillators. This carries important implications for the interpretation of the mechanisms underlying the coupling functions, as will be discussed below. Each function qi(κ) is periodic, for κ ≥ 2 on the κ-dimensional torus, and can be decomposed into a sum of κ-dimensional Fourier series of trigonometric functions. In practice it is assumed that the dynamics can be well-described by a finite number K of Fourier terms (Kralemann et al., 2011; Duggento et al., 2012): ϕ˙i=∑k=-KKck(i)Φi,k(ϕ1,ϕ2,…,ϕN)+ξi=∑k=-KKck(i)exp[i(k1ϕ1+k2ϕ2+…+kNϕN)]+ξi, where i = 1, …, N, Φ_{i, 0} = 1 so that c0(l)=ωl, and the rest of Φ_{i, k}, scaled by ck(i), are the k most important Fourier components. Such Fourier series Φ_{i, k}(ϕ₁, ϕ₂, …, ϕ_N) act as base functions for the dynamical inference method.

Our aim is to reconstruct a dynamical model describing the interactions through the analysis of data, so that the model can then be used for extraction of the coupling functions. Our approach is based on the method of dynamical inference, often referred to as dynamical modeling or dynamical filtering (Kalman, 1960; Sanjeev Arulampalam et al., 2002; Friston et al., 2003; Voss et al., 2004; von Toussaint, 2011).

Note that inference of cross-frequency couplings from the statistics of the coupled signals, e.g., through correlation, (bi-)coherence and Granger causality measures (Geweke, 1982; Baccala and Sameshima, 2001; Kamiński et al., 2001), yields the functional connectivity but it provides no information about the mechanisms of causality. These latter methods are designed to infer statistical effects rather than dynamical mechanisms (Barrett and Barnett, 2013). In what follows, however, with the aid of dynamical inference we discuss how the mechanisms of the associated causality can be inferred from data, thus yielding an effective connectivity (Friston, 2011).

In particular, coupling functions represent one type of dynamical mechanism and their inference yields the effective connectivity. More specifically, the form of the coupling function defines the functional law under which some input of the interactions (i.e., the mutual influence between the oscillations) is translated into an appropriate output. This is related, not only to the quantitative parameters of the net coupling strength i.e., net information flow, but also to how this information is functionally structured to give an effective mechanism. For example, as we will see below, the interactions can be such that the form of the coupling function varies in time (see e.g., Section 3.3.2 and Stankovski et al., 2012). This dynamical change can cause a qualitative transition (like synchronization), irrespectively of the value and the variations of the net coupling strength. This is an example of a case where functional connectivity methods (e.g., Granger causality) will detect only the net coupling strength and not the possible reason for a qualitative transition, unlike coupling functions analysis which can do so (see below).

A number of different techniques are available for estimating a model from data, based on different procedures and theories, and resulting in slightly different properties and characteristics. They include e.g., least-squares and kernel smoothing fits (Rosenblum and Pikovsky, 2001; Kralemann et al., 2013), dynamical Bayesian inference (Smelyanskiy et al., 2005; Stankovski et al., 2012), maximum likelihood (multiple-shooting) methods (Voss et al., 2004; Tokuda et al., 2007), and dynamic causal modeling (Friston et al., 2003).

In what follows we use the dynamical Bayesian inference technique (Smelyanskiy et al., 2005; Stankovski et al., 2012). Briefly, the method applies Bayesian probability theory to the multidimensional time-series to infer the dynamical model in terms of stochastic differential equations. Assuming a normal multivariate distribution for the prior of the scale parameters, by the use of the model base functions, the method constructs a log-likelihood function which also ensures that the posterior probability is normally distributed. Evaluation of the current distribution relies on the evaluation of the previous block of data in the sequence, i.e., informative priors are used and the current prior depends on the previous posterior. For the first time window, in the absence of an earlier block, we set the initial prior to a flat (zero) distribution; which might effect the precision with which the initial coupling function is inferred for that window. To account for the time-variability of the interacting dynamics, the covariance matrix of the next prior is the convolution of the current posterior with the current diffusion matrix which describes how much the parameters can change. Further details of the method can be found in the Supplementary Material, in Smelyanskiy et al. (2005), Stankovski et al. (2012), Duggento et al. (2012), and Stankovski et al. (2014a) and in the references therein.

Using the inferred parameters we can calculate the coupling quantities and characteristics. The coupling functions q_i(ϕ_i, ϕ_j, ϕ_k, …, ϕ_N) acting on the oscillator from each of the i phases are evaluated on a 2π × 2π × … × 2π grid by selecting the relevant base functions, i.e., Fourier components scaled by the corresponding inferred coupling parameters. The coupling strength is calculated as the Euclidean norm ||qi||=〈qiqi〉1/2 of the inferred parameters for a particular coupling, and therefore carries the same unit of measure as the natural frequency (Hz). The correlation ρi(qi,qj)=〈q~iq~j〉/(||q~i|| ||q~j||), of two coupling functions where q~i are the deviations from the mean, q~i=qi-〈qi〉, gives the similarity of their forms, irrespectively of their amplitudes (Kralemann et al., 2013). Here, we propose a further extension of this index. By calculating the correlation of a coupling function q with a sequence of numerically-generated forms Q having specific shape features, taken from a bank, one can determine which of those features is dominant in q. The numerical set simulates the shape of a direct coupling from the slower oscillation to the faster, phase-shifted by an angle ϑ along the 2π axes. Thus, the numerical form Q_ϑ generating the highest ρ carries dual information: the extent of the similarity (described by ρ itself) and the corresponding phase, given by ϑ. See Figure 1 and the animation video 1 in the Supplementary Material. A natural way of presenting this information is by plotting it on the complex plane to provide a polar representation of the similarity index Pq=ρqeiϑ.

In neuroscience, the cross-frequency analyses reported to date have mostly focused on the net coupling. In contrast, coupling functions enable one to study the functional dependences of the distinct contributions from the individual oscillations. This procedure, referred to as coupling decomposition, separates a net pairwise coupling q_i(ϕ_i, ϕ_j) into its partial self-coupling q̄i(ϕi), direct-coupling q̄i(ϕj), and common (or indirect) q̄i(ϕi,ϕj) coupling components (Iatsenko et al., 2013; Stankovski et al., 2016). The inference of both net and partial coupling has been validated numerically (Stankovski et al., 2015). The direct-coupling q̄i(ϕj) describes the influence of the direct unidirectional driving exerted by one oscillator on the other. Arguably, this is the most observed interaction in physiology, often linked to modulation mechanisms; it dominates in a number of the coupling functions discussed below. Similarly, for a triplet coupling function q_i(ϕ_i, ϕ_j, ϕ_k) one can decompose the self, direct, and common components depending on either one or two phase variables. Additionally, one can have the direct component q̄i(ϕj,ϕk) from two phase variables exerting a joint influence, and the common component between all three phases q̄i(ϕi,ϕj,ϕk). Generalization to higher κ-dimensional couplings is implicit. These couplings in a κ-dimensional network could reflect a joint functional influence from a cluster subnetwork.

The multichannel EEG recordings analyzed in this work were downloaded from the Neurophysiological Biomarker Toolbox (NBT) dataset (O'Gorman et al., 2013; Poil et al., 2013). The signals were recorded for a group of 16 subjects (of which 10 were female, median age 27 years, range 21–48) in the resting state for 8 min, with a sampling frequency of 200 Hz. During the first 4 min, subjects were asked to keep their eyes open, and in the following 4 min to keep them closed. Signals from 19 EEG electrodes corresponding to the international 10–20 system were selected from the dataset for the analysis.

The cross-frequency intervals were extracted by a standard (FIR and no-phase-shift) filtering procedure. The boundaries for the conventional frequency intervals were: delta δ = 0.8–4 Hz, theta θ = 4–7.5 Hz, alpha α = 7.5–14 Hz, beta β = 14–22 Hz, and gamma γ = 22–40 Hz. Special care was taken to minimize cardiac components and powerline interference (Lehnertz et al., 2014; Iatsenko et al., 2015). The phases of the filtered δ and α were estimated by use of the Hilbert transform, followed by the protophase-phase transformation (Kralemann et al., 2008).

The extensive changes that the simple closing of the eyes triggers in the brain caught the attention of the very first electroencephalographers (Berger, 1933). It is now known that exclusion of visual input from the central system causes the power of brain activity to increase instantaneously across all the conventional frequency ranges (Barry et al., 2007). The most striking change occurs within the α rhythm, and it has its strongest effect on the occipital part of the scalp, over the visual cortex area. It has been argued that, with eyes open, the desynchronization of α, resulting in a lower power, might occur in order to give way to a more sophisticated pattern of information processing (Klimesch, 1999).

The δ-to-α interaction reflects how δ activity, associated with deep dreamless sleep (Feinberg et al., 1987), influences the α oscillations related to information processing (Pfurtscheller and Lopes da Silva, 1999). Other findings have also suggested that the δ-to-α coupling is mostly located within the frontal regions, that it is stronger during the eyes-closed resting state (Deco et al., 2010; Jirsa and Müller, 2013), and that a strong δ-to-α link exists during non-REM sleep (Bashan et al., 2012).

Cross-frequency interactions are usually mediated by the slower oscillations modulating the faster ones (Brunel and Wang, 2003; Lakatos et al., 2005; Händel and Haarmeier, 2009). In particular, task-based studies suggest that slow oscillations, which are extended across the scalp, modulate the spatial extent of the faster oscillations, which are more localized (Palva et al., 2005; Isler et al., 2008; Canolty and Knight, 2010).

In the light of this, and because of the crucial role that the α oscillation (Klimesch et al., 2007; Eidelman-Rothman et al., 2016) plays in the eyes open (EO) and eyes closed (EC) states, we focused on the analysis of δ-to-α coupling functions. In doing so, we are able to assess, quantify, and describe in detail the functional mechanisms that define the interaction in question.

Moreover, the multichannel recordings allowed us to investigate couplings between δ and α oscillations extracted from different probes, and hence to create connectivity maps illustrating how the δ-to-α modulation differs in the EO and EC states. The coupling strength was first quantified. Note that, in earlier work (Musizza et al., 2007; Jirsa and Müller, 2013; Lehnertz et al., 2014) the use of the terms “coupling causality” and “directionality” refers to the net coupling strength.

When applying non-linear analysis techniques, one should bear in mind that the linear properties of the signals, like autocorrelation or spectral features, are likely to affect the measure. To discriminate the genuine results from the ones happened by chance, one can apply surrogate testing (Theiler et al., 1992; Schreiber and Schmitz, 1996, 2000; Paluš and Hoyer, 1998; Kreuz et al., 2004). The idea behind this technique is to apply the non-linear method in question to independent time series that have the same, or as close as possible, statistical properties as the original time signals, while randomizing the expressions of the non-linear property being measured. This procedure allows one to define a threshold beneath which any result is considered spurious.

In practice, when inferring couplings even from very weakly-coupled (or completely uncoupled) systems, the methods always detect some non-zero values of apparent coupling strength. Surrogate testing can then be used to establish the “zero-level” of apparent coupling corresponding to uncoupled signals. In order not to bias the threshold with effects due to inter-subject or inter-probe variability, we applied the surrogate techniques to the same signals for which the coupling was to be measured, and we therefore define different thresholds for different subjects, pairs of probes and states.

We generated the necessary surrogates by use of the phase-shuffling (PS) method (Schreiber and Schmitz, 2000; Jirsa and Müller, 2013). This acts on the time evolution of the phase of an oscillation, wrapped between 0 and 2π, by randomizing the sequence of full phase-periods that it contains. With this technique, the linear structures of the signals are preserved but the nonlinear properties are changed. Non-stationarities appearing within each period of the oscillations are preserved. The method was applied for each subject, state, and pair of probes, thereby providing pairs of surrogate phases (δ and α). These pairs were used as input for the Bayesian inference to compute the surrogate coupling. The significance thresholds, calculated independently for each subject and combination of probes, were then set as the mean+2 standard deviations of the resultant distributions.

The surrogate populations were tested for normality with the Shapiro-Wilk test, with the null hypothesis that the data come from a normal distribution of unknown mean and variance. The test rejected the null hypothesis at the 5% significance level in only 3% of the surrogates, and we therefore accepted the assumption of a normal distribution. Hence, we could test the coupling from the original signal by comparison with the significance threshold.

The non-parametric Wilcoxon paired test was used to determine the significance of differences between the EO and EC distributions for each frequency within the power spectra, for the averaged power within each frequency interval, for the coupling strength and for the similarity of coupling functions.

Figure 2A shows the difference in spectral power between the EO (blue) and EC (red) states, for spectra averaged across all the probes. The shaded significance area coincides closely with the α band, indicating an increase of power in that interval for EC compared to EO. This increase was independent of scalp location. Figure 2B shows the statistical distributions of the averaged power within each frequency band, with a pairwise probe-by-probe statistical approach. The statistical analysis confirmed the increase of amplitude across all the frequency intervals when comparing EC with EO.

To investigate the quantitative statistics of each group of subjects we calculated the average values of the significant coupling strengths, with the corresponding surrogates' value subtracted, and the moduli of the polar similarity indices for the coupling functions of all the links for each subject. Then we compared statistically the distributions of these values for the two groups of subjects. To present the differences between the distributions visually, we use standard boxplots which refer to the descriptive statistics (median, quartiles, maximum and minimum).

The results in Figure 7 show that there were statistically significant differences for both the coupling strengths and the similarity of coupling functions between the EO and EC states. Figure 7A shows that the coupling strengths detected for the EC is significantly higher than the EO. Similarly, Figure 7B shows that the similarity index for δ-to-α coupling functions for the EC were significantly higher than for the EO. The latter also means that there was larger variability of the coupling functions for the EO state, compared to the EC state. Overall, the similarity of coupling functions for the EO and EC states was not very high (in the interval of [0,1]), indicating that there is relatively high variability of coupling functions for both of the resting states.

Much has already been done, mostly through fMRI and EEG analysis, to demonstrate the existence of resting state interactions, including the formation and dissolution of resting state functional network configurations around a stable anatomical connectivity (Berger, 1933; Barry et al., 2007; Deco et al., 2010; Jirsa and Müller, 2013; O'Gorman et al., 2013). Our application of coupling functions to the resting state revealed the underlying mechanisms of interaction and has identified a number of differences between the EO and the EC states.

As there were more significant couplings in the EC than in the EO state (Figure 3), it is obvious that there will be more coupling links for the EC than for the EO state when presented spatially (Figure 4). What is interesting is that, for EC, different subjects seem to have a preferential pattern of directions, with the δ oscillation from the anterior temporal lobes (probes T3 and T4) acting as “hubs,” influencing the phase of α in both the frontal and occipital directions. The occipital probes O1 and O2 are the most susceptible to the difference in α power (Figure 2) as they are placed over the visual area of the cortex. In EO, they act as a starting point for δ modulating long range connections toward the frontal cortex, which existed in five subjects, and then disappeared in EC. In contrast, for EC these probes receive the influence in their α rhythm from temporal and central probes.

The δ-to-α coupling functions had a specific shape, showing that the coupling is predominantly like a direct sine wave due to the δ influence, which accelerates and decelerates the α oscillations. Importantly, the form was similar for the EO and EC states (Figure 5), with distinctive variations and shifts along the δ oscillation. This similarity implies that the same underlying interaction mechanism exists in the EO and EC states, and that the difference between these two resting states corresponds to increasing and decreasing some of the connection strengths (or to switching them on-off).

Because we reconstructed the form of the coupling functions, we were able to observe what they look like for both individual and averaged connections and subjects. Even though we found relatively similar forms of function, we also observed a certain degree of variability, both inter-subject variability (Figure 5) and time variability (Figure 6) of the form. These should be taken into account when average values are used, for example in making multi-subject statistics.

Finally, for the comparison of the EO and EC states (Barry et al., 2007) our analysis confirmed that the spectral power of the α oscillations in EC is significantly larger than that of EO (Klimesch, 1999). It also showed that there are a larger number of real (i.e., validated by surrogate testing) δ-to-α couplings for the EC state (Jirsa and Müller, 2013), that the form of the coupling functions was similar for EO and EC, that the coupling functions were somewhat less variable for EC than for EO, and that this dominance of the EC state in the interactions was confirmed also by the quantitative boxplot statistics for the whole groups of subjects.

The assessment of neural coupling functions through the phase dynamics of interacting neural oscillations enables us to study their acceleration/deceleration, i.e., timing and coordination. The generalization to amplitude coupling functions is implicit. In such cases, one should be able to determine a plausible state model in relation to the dimensionality of the signals. Amplitude neural coupling functions can reveal the mechanism through which the strength and power of one neural oscillation are affected by the influence of the other oscillations.

Earlier effective connectivity methods for the inference of neural interactions have in principle contained coupling functions within their models of the interacting dynamical systems. The question we address here, in addition to presenting an efficient Bayesian method for determination of coupling functions, is that of how to assess the neural coupling functions. We have shown how to unify a functional unit which can be quantified and compared with other such units, and whose evolution can be followed in time. The key characteristic that distinguishes this assessment is the form of the neural coupling functions. A unified and effective coupling function analysis can provide insights that go far beyond just knowing that neural interactions exist.

The pairwise investigation can further be generalized to higher degrees of network complexity (Kralemann et al., 2014; Stankovski et al., 2015). One might, for example, study the coupling functions between the brain and other physiological oscillations, forming a physiological network (Musizza et al., 2007; Stefanovska, 2007; Bashan et al., 2012; Stankovski et al., 2016). The brain is a heavily connected network (Park and Friston, 2013) and coupling functions could be applied to reveal the functional mechanisms operative at different levels and sublevels of the interactions. In network topology with nodes and edges (Albert and Barabási, 2002) this would mean that, not only could the existence, strength and direction of the edge be studied, but also the underlying functional mechanism giving rise to the edge. The multivariate coupling function assessment can then be linked to hypergraphs (Karypis and Kumar, 2000; Zass and Shashua, 2008), though it was argued recently that, for larger networks (N > 10), there is no significant benefit from using multivariate inference of coupling (functions) and partialization (Rings and Lehnertz, 2016).

The time-varying form of the coupling functions (Figure 6) can be a cause of self-organization transitions, like the emergence of network clustering, or of the systems going into-and-out-of synchronization (Stefanovska et al., 2000; Varela et al., 2001), even for an invariant net coupling strength (Stankovski, 2017). More importantly, having detected and characterized a neural coupling function, one can then use this knowledge to detect, or even to predict, the onset of phase synchronization (Kiss et al., 2005). In such cases, the key feature is the known form of the coupling function which, depending on parameters like frequency, coupling strength, or polar similarity index, can predict the synchronization transition. This could have important implications for the prediction of epileptic seizures (Lehnertz and Elger, 1998; Fell et al., 2001) which occur or disappear as synchronous activity in the brain.

In conclusion, coupling functions bring a novel perspective to neuroscience that is unique in that it provides access to the functional form of a coupling. The polar similarity index that we have introduced allows one to describe the form in quantitative detail. The comparisons of δ-to-α phase-to-phase coupling functions in the EO and EC resting states demonstrate how neural coupling functions can be reconstructed from spatially distributed sources, and what benefits and possibilities are introduced by their assessment. We have confirmed the previous result that the direct coupling is stronger during EC, and we have shown for the first time that the coupling function is significantly less variable in that state. The EO/EC states were taken as an example on which to base a discussion of methodological issues and, in so doing, to point to the wider implications and possibilities of the method itself. One may hope to gain new insights into the neuronal mechanisms underlying certain diseases from studies of coupling functions. In principle, the method can equally be applied to the time series created by any pair of coupled oscillatory processes.

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.