Time-varying coupling functions: dynamical inference and cause of synchronization transitions

TL;DR

This paper investigates how time-varying coupling functions, independent of coupling strength, can induce synchronization transitions in biological systems, demonstrated through EEG analysis and phase oscillator models.

Contribution

It introduces a method to detect synchronization caused solely by changes in coupling functions using dynamical Bayesian inference.

Findings

Time-varying coupling functions are present in biological interactions.

Synchronization transitions can occur without changes in coupling strength.

Dynamical Bayesian inference effectively detects time-varying coupling functions.

Abstract

Interactions in nature can be described by their coupling strength, direction of coupling and coupling function. The coupling strength and directionality are relatively well understood and studied, at least for two interacting systems, however there can be a complexity in the interactions uniquely dependent on the coupling functions. Such a special case is studied here { synchronization transition occurs only due to the time-variability of the coupling functions, while the net coupling strength is constant throughout the observation time. To motivate the investigation, an example is used to present an analysis of cross-frequency coupling functions between delta and alpha brainwaves extracted from the electroencephalography (EEG) recording of a healthy human subject in a freerunning resting state. The results indicate that time-varying coupling functions are a reality for biological interactions. A model of phase oscillators is used to demonstrate and detect the synchronization transition caused by the varying coupling functions, during an invariant coupling strength. The ability to detect this phenomenon is discussed with the method of dynamical Bayesian inference, which was able to infer the time-varying coupling functions. The form of the coupling function acts as an additional dimension for the interactions and it should be taken into account when detecting biological or other interactions from data.

Equations

$$\begin{array}{rl}{\dot{\varphi }}_1& ={\omega }_1+{\epsilon }_1(t)q_1({\varphi }_1,{\varphi }_2)=\\ & ={\omega }_1+c_1(t)\sin ({\varphi }_2-{\varphi }_1)+c_2(t)\cos ({\varphi }_2-{\varphi }_1)\\ & \\ {\dot{\varphi }}_2& ={\omega }_2+{\epsilon }_2(t)q_2({\varphi }_1,{\varphi }_2)=\\ & ={\omega }_2+c_3(t)\sin ({\varphi }_1-{\varphi }_2),\end{array}$$

$${\dot{\varphi }}_i={\omega }_i+q_i({\varphi }_i,{\varphi }_j)+{\xi }_i(t),$$

$${\dot{\varphi }}_i=\sum _{k=-K}^K{\tilde{c}}_k^{(i)}(t){\Phi }_{i,k}({\varphi }_i,{\varphi }_j)+{\xi }_i(t),$$

$$\begin{array}{rl}\mathbf{D}& =\frac{h}{L}{\left({\dot{\varphi }}_n-{\mathbf{c}}_k{\mathbf{\Phi }}_k({\varphi }_{\cdot ,n}^{\ast })\right)}^T\left({\dot{\varphi }}_n-{\mathbf{c}}_k{\mathbf{\Phi }}_k({\varphi }_{\cdot ,n}^{\ast })\right),\\ {\mathbf{r}}_w& =({\mathbf{\Xi }}_{\text{prior}}{)}_{kw}{\mathbf{c}}_w+h{\mathbf{\Phi }}_k({\varphi }_{\cdot ,n}^{\ast })({\mathbf{D}}^{-1}){\dot{\varphi }}_n+\\ & -\frac{h}{2}\frac{\partial {\mathbf{\Phi }}_k({\varphi }_{\cdot ,n})}{\partial \varphi },\\ {\mathbf{\Xi }}_{kw}& =({\mathbf{\Xi }}_{\text{prior}}{)}_{kw}+h{\mathbf{\Phi }}_k({\varphi }_{\cdot ,n}^{\ast })({\mathbf{D}}^{-1}){\mathbf{\Phi }}_w({\varphi }_{\cdot ,n}^{\ast }),\\ \stackrel{~}{}{\mathbf{c}}_k& =({\mathbf{\Xi }}^{-1}{)}_{kw}{\mathbf{r}}_w,\end{array}$$

Full text

I Introduction

Interacting dynamical systems abound in nature, with examples ranging from physics, biology, climate and social sciences Pikovsky et al. (2001); Strogatz (2003); Haken (1983). Many of them are open systems and have connections with other systems and the environment. The connectivity presents a link between two dynamical systems, which can be structural, functional or effective Friston (2011). Very often these interactions are assessed successfully to a relatively great extend when they are considered isolated, however in some cases there is additional complexity due to external influences and the existing time-variability. From the aspect of mathematical models, such systems and their interactions are studied as non-autonomous dynamical systems Kloeden and Rasmussen (2011); Stankovski (2013); Suprunenko et al. (2013); Clemson and Stefanovska (2014); Lancaster et al. (2015). The time-variability can have different effects on the interactions, including for example, changes in frequency, emergence or disappearance of connectivity, transitions to or out of qualitative states, or time-varying form of the coupling functions.

Coupling function describes in great detail the physical rule of how the interactions occur and manifest. For example, the reconstructed coupling function of Belousov-Zhabotinsky chemical interactions revealed how there could be higher-harmonics and bi-stability of the synchronization state Miyazaki and Kinoshita (2006), the knowledge of the coupling function of one pairwise interaction was used to predict the synchronization and clustering of a network of electrochemical oscillators Kiss et al. (2005), and the form of the cardiorespiratory coupling function was linked to respiratory sinus arrythmia (RSA), a known mechanism in physiology Iatsenko et al. (2013); Kralemann et al. (2013). The coupling function as a whole can be described in terms of its strength and form. It is the functional form that has provided a new dimension and perspective, probing directly the functional mechanisms of the interactions. In this way the coupling function can determine the possibility of qualitative transitions between states of the systems e.g. routes into and out of synchronization. Decomposition of a coupling function can also facilitate a description of the functional contributions from each separate subsystem within the coupling relationship. Different methods for coupling function detection have been applied widely in chemistry Kiss et al. (2007); Miyazaki and Kinoshita (2006); Tokuda et al. (2007); Kiss et al. (2005); Blaha et al. (2011), in cardiorespiratory physiology Kralemann et al. (2013); Stankovski et al. (2012); Iatsenko et al. (2013), in neuroscience Stankovski et al. (2015); Stankovski et al. (2016a); Wilting et al. (2015), in mechanical interactions Kralemann et al. (2008), in social sciences Ranganathan et al. (2014) and in secure communications Stankovski et al. (2014a). The study of coupling function is a very active and expanding field of research Stankovski et al. (2016b).

Recently, it was found that in interacting biological oscillating systems, not only the frequency and the coupling strength, but also the form of the coupling functions can be a time-varying process Stankovski et al. (2012). This was demonstrated in the case of cardiorespiratory interactions and showed that not only the net parametric and quantitative properties, but also the mechanisms of the interactions can be time-varying. These varying coupling functions can change in time the physical rules for the interactions, which can cause transitions of physical effects and phenomena, as the most important outcome of the interactions. For example, the time-varying cardiorespiratory coupling function was shown to induce transitions in and out of synchronization, and between different synchronization ratios Stankovski et al. (2012). Therefore, understanding the effects of the time-varying coupling functions on the interactions is of great interest, especially in understanding, detecting and interpreting the interactions of open (biological) systems. In this paper, we further investigate the effects of the time-varying coupling functions on the interactions which can be revealed uniquely by the coupling function analysis.

The time-varying coupling functions are especially important for the detection and inference of interactions from data. Namely, they can be represented by a series of sub-coupling components and their time-variability can introduce a certain complexity in the interactions. Being able to detect and correctly interpret the coupling will depend greatly on the nature of the methods used. For example, some methods assess the amount of information, the net coupling strength and directionality, while other methods perform dynamical inference of a model with differential equations and can describe the underlying mechanisms and coupling functions. Hence, one of the main aspects of the paper will be to discuss a method for analysis of the interactions, in the light of the complexity introduced by the time-varying coupling functions.

II Biological interactions as motivation – the case of neural coupling functions

Biological systems are of great importance in this field as they are a classical example of systems that are not isolated but interact with other systems in the body and can be affected by other forces in the local environment. Therefore, the correct analysis of such systems has great implications for the determination and treatment of various physiological states and diseases.

The electrophysiology of the brain evaluated through electroencephalography (EEG) presents an important characteristic for investigating different neural states and diseases, however there could be a certain complexity due to time-varying coupling functions when one tries to analyze the neural interactions associated with this electrophysiology. Here we present the neural cross-frequency coupling functions Stankovski et al. (2015) as an example of biological systems whose coupling functions are time-varying. For this purpose, the EEG signal of a human subject in a resting state from the public PhysioNet database is used Goldberger et al. (2000); phy ; Schalk et al. (2004). The EEG signals were recorded with 64 electrodes according to the standard international 10-10 system. Only one EEG signal recorded at the frontal Fp1 electrode of the 10-10 system, recorded during eyes-closed resting state of a healthy subject, was analysed here. After extracting the phases for the delta and alpha brainwave oscillations of the filtered delta and alpha brainwave signals from the EEG signal, the neural cross-frequency coupling functions and the coupling strength of two components are reconstructed by dynamical Bayesian inference (for details about the method see Sec. IV.1 below). The two sub-coupling components were chosen arbitrarily for better presentation and in accordance with the other numerical examples that follow.

The delta-to-alpha neural cross-frequency coupling has been found to be generally higher in the eyes-closed than in the eyes-open condition of the resting state, mostly located within the frontal (e.g. at Fp1) and the parieto-occipital regions, and these regions were connected through larger-scale coupling with a different coupling direction Jirsa and Müller (2013); the delta-alpha coupling was significantly increased with more similar forms of the coupling functions from awake to deep general anaesthesia, and this effect was higher for anaesthesia induced with sevoflurane than with the propofol anaesthetic Stankovski et al. (2016a); and a strong link between delta and alpha brain activity was found during non-REM sleep, although alpha waves were greatly diminished and delta waves were dominant Bashan et al. (2012). Fig. 1 (a) shows that two of the delta-alpha sub-coupling components have significant time-variability between the allowed $[-1,1]$ values. This leads to the time-variability of the form of the coupling functions – as can be seen by comparing the three delta-alpha coupling functions in Fig. 1 (b). The fact that the sub-coupling components are time-varying, causing the form of the coupling functions to vary as well, points to the possibility of qualitative transitions of the interactions. Therefore, revealing this complexity can be of crucial importance for the correct analysis and interpretation of the neural interactions.

III The Characteristic problem of time-varying coupling functions– direct approach

The sub-coupling components define the form of the net function i.e. the coupling function is evaluated as a group sum of all the functional components of some set of decomposition functions. For example, for the phase dynamics of interacting oscillators, because of the periodicity, one can decompose the coupling functions into Fourier series. Therefore, the variations of some of the sub-coupling components will also define how the form of a net coupling function will vary. The latter effect, may or may not change the coupling strength, to a greater or lesser extent, however it will change the form of the coupling function and thus the mechanism underlying the interactions, which in turn can cause qualitative transitions. The mechanism is defined by the function that gives the rule through which the input values are translated into output values Friston (2011); Barrett and Barnett (2013).

To further explain and emphasize the complexity and the effects from the varying sub-coupling components, a simple numerical example is used. A system of two interacting phase oscillators Kuramoto (1984); Petkoski and Stefanovska (2012) is considered:

[TABLE]

where $\omega_{1}$,$\omega_{2}$ are the natural frequencies, while $c_{1}$,$c_{2}$,$c_{3}$ are the strengths of the sub-coupling components that scale the appropriate coupling function components. The numerical values of the parameters were set to $\omega_{1}=2.28$, $\omega_{2}=1.2$ and $c_{3}(t)=0.1$. Here, and throughout the manuscript, the focus is on the coupling in only one direction – the coupling from the second to the first phase oscillator. In this coupling direction the strengths of the sub-coupling components $c_{1}(t)$ and $c_{2}(t)$ are set to be time-varying i.e. $c_{1}(t)=\sqrt{t/T}$ and $c_{2}(t)=\sqrt{(T-t)/T}$, where the time changes $t=0\rightarrow T$ and $T$ is the total time of observation. In this way $c_{1}(t)$ varies in the interval $[0,1]$, and $c_{2}(t)$ varies in the interval $[1,0]$. The total net coupling strength, or the net information flow in one coupling direction, is evaluated as the Euclidian norm of all the coupling components Kralemann et al. (2008), which in this case is $\varepsilon_{1}(t)=\sqrt{c_{1}^{2}(t)+c_{2}^{2}(t)}.$

The main purpose of the example is to present constant net coupling strength, but time-varying coupling functions. To do this $c_{1}(t)$ and $c_{2}(t)$ were varied in the interval $[0,1]$ simultaneously, but reducing one and increasing the other, in such a way as to give constant net coupling $\varepsilon_{1}(t)$ for all time. Fig. 2 (a) shows the time-variation of the coupling parameters, with $c_{1}(t)$ and $c_{2}(t)$ varying, while $\varepsilon_{1}(t)$ is being constant. The important part is that even though the net coupling is constant, the coupling function also varies due to the variations of the sub-coupling components.. The latter can be observed by comparing the different forms of the coupling functions at the beginning Fig. 2 (b) and the end of the observation time Fig. 2 (c).

The variation of the coupling function changes the mechanism and the physical rule under which the interactions are manifesting, which can lead to qualitative transitions into or out of certain physical effects and phenomena of the interactions. The latter could include transitions to synchronization, amplitude or oscillator death, clustering in networks or the emergence of chimera states. With the current example we present the transition to synchronization. Fig. 2 (d),(axis right) shows that the phase difference of the interacting system (1) is not bounded at the beginning, and then at approximately $t\sim 2500s$ there is a transition to bounded phase difference, as in the case of synchronized systems. The same can be verified by the synchronization detection index Tass et al. (1998), which quantifies, in the interval $[0,1]$, the synchronization (or phase coherence) from the phase difference. The method first divides each phase interval into $N$ bins ($N=10$ used here, and window size 50s). Then, for each bin $l$, it calculates the dependence of $\phi_{2}(t_{j})$ on $\phi_{1}(t_{j})$, such that $\phi_{1}(t_{j})$ belongs to this bin $l$, and $M_{l}$ is the number of points in the bin. The average over all bins leads to the final index: $I_{Sync}(t_{j})=1/N\sum_{l=1}^{N}|M^{-1}_{l}\sum e^{i\phi_{2}(t_{j})}|.$ In this way, $I_{Sync}$ measures the conditional probability for $\phi_{2}$ to have a certain value provided $\phi_{1}$ is in a particular bin. If the values are greater than $\sim 0.95$ (surrogate level resulting from the mean plus two standard deviations of synchronization indexes from phase random shuffling realizations Schreiber and Schmitz (2000); Stankovski et al. (2014b), shown with dashed lines in Fig. 2 (d)), synchronization is detected. Fig. 2 (d) (axis left) shows consistently that at approximately $t\sim 2500s$ there is a transition to synchronization. Therefore, although the net coupling strength was constant between the systems, due to the change of the coupling functions there is a transition to synchronization.

IV Interaction analysis – inverse approach

The main objective here is to be able to detect or infer from data the coupling relations between the systems. This constitutes the so-called inverse approach, starting from the data and attempting to learn the nature of the connectivity between dynamical systems. The methods for coupling inference have great importance as they allow different applications in various fields.

Many efficient methods for coupling detection have been designed Smirnov and Bezruchko (2009); Andrzejak et al. (2006); Paluš and Stefanovska (2003); Bahraminasab et al. (2008); Jamšek et al. (2010); Hlaváčkováá-Schindler et al. (2007); Faes et al. (2011); Stankovski et al. (2012); Kralemann et al. (2008); Kiss et al. (2005); Tokuda et al. (2007), and they are based around two main aspects. The first aspect is that the assessment of the strength of the interaction and its predominant direction can be used to establish if certain interactions exist at all. In this way, one can determine whether some apparent interactions are in fact genuine, and whether the systems under study are really connected or not. The second aspect is that one can infer the underlying interacting systems, including the appropriate coupling functions, and learn about how an interaction occurs. The former aspect determines the existence of coupling and relates to functional connectivity, while the latter describes the mechanisms and is part of the effective connectivity methods Friston (2011). Below, an effective connectivity method which infers coupling functions inherently is presented in light of the interaction complexity, as discussed in the previous sections.

IV.1 Dynamical Bayesian Inference

Dynamical inference of coupling functions is a class of effective connectivity methods for coupling assessment. The main pillar of the procedure is a method for dynamical inference, often referred to as dynamical modelling or dynamical filtering von Toussaint (2011); Kalman (1960); Sanjeev Arulampalam et al. (2002); Voss et al. (2004); Stankovski et al. (2014b). The dynamical inference procedure starts with the data from two (or more) interacting dynamical systems and uses a method that infers a dynamical model in terms of (ordinary or stochastic) differential equations. The coupling functions are an integral part of, and can be extracted from, the inferred model.

There are a number of other methods for dynamical inference and coupling functions assessment, including those based on least squares and kernel smoothing fits Rosenblum and Pikovsky (2001); Kralemann et al. (2013), dynamical Bayesian inference Stankovski et al. (2012), maximum likelihood (multiple-shooting) methods Tokuda et al. (2007), stochastic modeling Schwabedal and Pikovsky (2010) and the phase resetting curve Levnajić and Pikovsky (2011). Below, the dynamical Bayesian inference Stankovski et al. (2012) will be presented and applied.

The signals under consideration are oscillatory and their interactions can be studied effectively through their phase dynamics. Therefore, a model of two coupled phase oscillators Kuramoto (1984) described by the stochastic differential equation is considered:

[TABLE]

with $i\neq j$ for $i,j=\{1,2\}$ and where $\omega_{i}$ is the parameter for the natural frequency. The deterministic part given by the base functions $q_{i}(\phi_{i},\phi_{j})$ describes the self and the interacting dynamics. The external stochastic dynamics $\xi_{i}(t)$ is considered to be Gaussian white noise $\langle\xi_{i}(t)\xi_{j}(\tau)\rangle=\delta(t-\tau)D_{ij}$. Due to the periodic nature of the deterministic dynamics, the base functions can be decomposed into an infinite Fourier series $q_{i}(\phi_{i},\phi_{j})=\sum_{s=-\infty}^{\infty}\sum_{r=-\infty}^{\infty}\tilde{c}_{i;r,s}\,e^{i2\pi r\phi_{i}}e^{i2\pi s\phi_{j}}$. In practice, however, the dynamics are well-described by a finite number of Fourier terms, so that one can rewrite the phase dynamics as:

[TABLE]

where $\tilde{c}^{(i)}_{0}=\omega_{i}$, and the rest of $\Phi_{i,k}$ and $\tilde{c}^{(i)}_{k}$ are the $K$ most important Fourier components; here we used $K=2$. It is important to note that the net coupling function $q_{i}(\phi_{i},\phi_{j})$ is decomposed into a series of sub-coupling components $\Phi_{i,k}(\phi_{i},\phi_{j})$, which in turn allows for the particular complexity of the interactions to be revealed. In this way, a coupling influence is separated on a large but finite set of sub-functional elements, making it possible to study the coupling contributions of each of them separately or in sub-groups. The assessment of the separate components defines the form of the net coupling function. The Fourier components $\Phi_{i,k}$ act as base functions for the dynamical Bayesian inference, through which the parameters $\tilde{c}_{k}^{(i)}$ are evaluated.

Dynamical Bayesian inference Stankovski et al. (2012); Smelyanskiy et al. (2005); Stankovski et al. (2015) enables one to evaluate the model parameters $\tilde{}{\bf c}$, which give the time-evolving coupling functions and coupling strength in the presence of noise. From Bayes’ theorem one can derive the minus log-likelihood function, which is of quadratic form. Assuming that the parameters are represented as a multivariate normal distribution (with mean $\bar{{\bf c}}$, and covariance matrix ${{\bf\Sigma}\equiv\Xi^{-1}}$), and given such a distribution for the prior knowledge using the likelihood function, one can calculate recursively Stankovski et al. (2012); Smelyanskiy et al. (2005) the posterior distribution of the parameters $\tilde{}{\bf c}_{k}$ using only the following four equations:

[TABLE]

where summation over $n=1,\ldots,N$ is assumed, and summation over repeated indices $k$ and $w$ is implicit. We used informative priors and a special procedure for the propagation of information between consecutive data windows Stankovski et al. (2012), which permitted the inference parameters that varied with time (for implementation and usage see Stankovski et al. (2014c)). Once we have the inferred parameters $\tilde{}{\bf c}$, we can calculate the coupling quantities and characteristics. The coupling functions are evaluated on a $2\pi\times 2\pi$ grid using the relevant base functions i.e. Fourier components scaled by their inferred coupling parameters. The method has been also generalized and applied to networks of oscillators Duggento et al. (2012); Stankovski et al. (2015).

The method is applied on data generated from the model of two interacting phase oscillators (Equ. 1), with the same time-variability of the coupling sub-components $c_{1}(t)$,$c_{2}(t)$ and constant total net coupling $\varepsilon(t)$ (Fig. 3 (a)) and the frequency parameters set to $\omega_{1}=9$ and $\omega_{2}=1.2$ to ensure that there is no synchronization between the systems. The latter is important, if there is strong phase synchronization the coupling method will not be able to correctly infer the couplings Rosenblum and Pikovsky (2001). The inference of the model parameters is very precise and the inferred time-varying coupling sub-components and the net coupling strength are practically indistinguishable from the original simulated parameters – Fig. 3 (a). Thus, the dynamical inference provides the information of the whole model and the variations of the sub-function components. The latter define the form of the net coupling function, hence the method is able to follow its time-evolution as well – compare the coupling functions from left to right in Fig. 3 (b). As discussed, even though the coupling strength is constant these varying coupling functions can cause qualitative transitions in the interactions.

V Discussion and Conclusion

Much effort has been concentrated in the study of coupling functions with significant progress on theory, methods and application. However, many open questions relevant to the coupling functions still exist. One particular aspect is their assessment, or what can be done once the coupling functions have been determined. The effect of time-varying coupling functions relates and further extends this aspect.

The time-variability of the sub-coupling components and the coupling functions are quite pronounced in biological systems. This was first observed in cardiorespiratory interactions and demonstrated here with the neural coupling functions of the human resting state (Fig. 1). It is important to note that the variability of the biological coupling functions can depend on the physiological state or disease. For example, it was observed that the cardiorespiratory coupling functions are less varying in young compared to old human subjects Iatsenko et al. (2013); Kralemann et al. (2013). Similarly, it was found that general anaesthesia reduces the variability of the delta-alpha neural coupling functions Stankovski et al. (2016a). The complexity arising from the varying coupling functions and their analysis is thus of great importance for the assessment of different medical states and diseases.

The theoretical case presented shows the relevance of the time-varying coupling functions in a special situation when the net coupling is invariant and the transition to synchronization is only due to the variations of the form of the coupling functions (Fig. 2). The treatment of the same phenomenon with the non-autonomous theory is an interesting open question for future developments. The particular example (Eq. 1) was chosen to be simple and elementary, as it elaborates rather complicated concepts, however the same phenomenon can be investigated also with more complex systems, like limit-cycle oscillators or higher-dimension chaotic systems. Importantly, the phenomenon could have even more complex implications for networks Albert and Barabási (2002) and various methods for networks inference Timme (2007); Levnajić and Pikovsky (2010); Sysoev et al. (2014); Levnaji and Pikovsky (2014); Timme and Casadiego (2014). Also, the variability is quite simple – only in two sub-coupling components. In real systems, the variations may be less pronounced but spread across more sub-coupling components. Needless to say, the situation of having exactly constant coupling strength while having large coupling function variations is quite a special case. In reality it is expected that the form variations are only partially contributing to the change of the interactions. The coupling strength and the coupling function are strongly related and to some extend dependant, but as shown here they can also affect the outcome of the interactions independently. This is because the form of the coupling function acts as an additional dimension in the assessment of the interactions.

The information theoretic based methods assess the statistical dependencies between the signals from the interacting systems. The most prominent methods of this kind are based on Granger causality, transfer entropy, mutual information and symbolic transfer entropy Granger (1969); Barrett and Barnett (2013); Schreiber (2000); Paluš and Stefanovska (2003); Staniek and Lehnertz (2008). They are statistical measures that can determine the causal relation and the predominant direction of influence, thus measuring a directed functional connectivity. In this way, they usually reveal only the net coupling and direction, thus they are not able to detect the varying sub-coupling components. This is because these methods are designed to perform in such a way i.e.  designed to infer only the net statistical effects Barrett and Barnett (2013). In fact, they are very useful methods for identification of directed connectivity and interaction existence, especially in the case where a model of differential equations for the dynamical inference can not be reconstructed. The presentation in this paper points to the importance of the phenomenon of time-varying sub-coupling components, and hopefully it will stimulate further developments of some aspects of these methods. Here, it is worth noting out that there have recently been efforts to use information based methods to perform coupling decomposition Faes et al. (2015); Porta et al. (2015) by exploiting the conditional functional dependencies. These methods could advance the detection of certain sub-coupling relations, without exploiting a dynamical model.

The dynamical inference methods are by design able to infer a dynamical model including the coupling functions. This was demonstrated by the use of dynamical Bayesian inference which was able to reveal the sub-coupling functions time-variability (Fig. 3). The design and application of dynamical inference methods is rapidly evolving, promising to reveal even more interacting complexities, unique and dependant on the coupling functions. The latter is even more pressing, because as shown the varying coupling functions are a reality for biological systems, and they could have important implications for the interactions in general.

Acknowledgements.

The author thanks Tiago Pereira, Martin Rasmussen, Michael Thompson and Rachel Sparks for valuable discussions on the topic, and the Institute of Pathophysiology and Nuclear Medicine, Faculty of Medicine in Skopje, and the Nonlinear and Biomedical Physics group at Lancaster University for the research and financial support.