Coupling functions: Universal insights into dynamical interaction mechanisms

TL;DR

This review comprehensively discusses the theory, methods, and applications of coupling functions, which describe interactions in complex dynamical systems across various scientific disciplines.

Contribution

It provides a unified overview of recent advances in analyzing, understanding, and applying coupling functions in diverse fields, highlighting their universality and methodological developments.

Findings

Coupling functions reveal detailed interaction mechanisms.

Methods for detecting and reconstructing coupling functions are diverse.

Applications span physics, biology, neuroscience, and social sciences.

Abstract

The dynamical systems found in Nature are rarely isolated. Instead they interact and influence each other. The coupling functions that connect them contain detailed information about the functional mechanisms underlying the interactions and prescribe the physical rule specifying how an interaction occurs. Here, we aim to present a coherent and comprehensive review encompassing the rapid progress made recently in the analysis, understanding and applications of coupling functions. The basic concepts and characteristics of coupling functions are presented through demonstrative examples of different domains, revealing the mechanisms and emphasizing their multivariate nature. The theory of coupling functions is discussed through gradually increasing complexity from strong and weak interactions to globally-coupled systems and networks. A variety of methods that have been developed for the detection and reconstruction of coupling functions from measured data is described. These methods are based on different statistical techniques for dynamical inference. Stemming from physics, such methods are being applied in diverse areas of science and technology, including chemistry, biology, physiology, neuroscience, social sciences, mechanics and secure communications. This breadth of application illustrates the universality of coupling functions for studying the interaction mechanisms of coupled dynamical systems.

Figures

Tables

Table 1. Table 1: Different examples of coupling functions q 𝑞 q . These pairwise coupling functions (CFs) are considered in relation to the system: x ˙ = f ​ ( x ) + q ​ ( x , y ) ˙ 𝑥 𝑓 𝑥 𝑞 𝑥 𝑦 \dot{x}=f(x)+q(x,y) .

Type of CF	Model	Meaning	Reference
Direct	$q(x,y)=q\left(y\right)$	unidirectional influence	Aronson et al. (1990)
Diffusive	$q(x,y)=q\left(y-x\right)$	dependence on state difference	Kuramoto (1984)
Reactive	$q(x,y)=\left(\epsilon +i\beta \right)q\left(x-y\right)$	complex coupling strength	Cross et al. (2006)
Conjugate	$q(x,y)=q\left(x-Py\right)$	$P$ permutes the variables	Karnatak et al. (2007)
Chemical synapse	$q(x,y)=g\left(x\right)S\left(y\right)$	$S$ is a sigmoidal	Cosenza and Parravano (2001)
Environmental	$q(x,y)\approx \epsilon {\int }_0^t{e}^{-\kappa \left(t-s\right)}\left(x\left(s\right)+y\left(s\right)\right)𝑑s$	given by a differential equation	Resmi et al. (2011)

Table 2. Table 2: Comparison between classes of coupling function and the techniques to obtain synchronization. Dashes indicate that typically we are unable to construct the Lyapunov function in such cases.

Coupling Function	Class	Technique	Global	Persistence
$H$ +ve definite	Lyapunov	LF	Yes	Yes
$DH$ +ve definite	Lyapunov	LF	No	Yes
$H$ differentiable	generic	LF	$-$	$-$
$H$ differentiable	generic	MSF	No	No

Table 3. Table 3: The phase sensitiveness for various models and their bifurcation after Brown et al. ( 2004 ) . Izhikevich ( 2000 , 2007 ) obtained the phase sensitivity Z 𝑍 Z for relaxation oscillators. They are discontinuous and are not shown here.

Bifurcation	$Z\left(\varphi \right)$
SNIPER	$1-\cos \varphi$
Hopf	$\sin \left(\varphi -\beta \right)$
Homoclinic	$\exp \left(-\lambda \varphi \right)$
Integrate and Fire	$2\pi$
Leaky Integrate and Fire	$\exp \left(g\varphi \right)$

Table 4. Table 4: Comparison of methods for the inference of coupling functions in terms of four characteristics (columns) including, respectively: stochastic treatment; absence of parameter overfitting; calculation speed; and the size of data windows. The methods (rows) are as described in the previous sections: least square fitting (LSF); dynamical Bayesian inference (DBI); maximum likelihood estimation with multiple shooting (MLE-MS); random phase resetting (RPR); and the stochastic modelling of effective coupling functions (SMECF). The symbol ✓ ✓ \checkmark indicates if a method possesses that characteristic, and a × \times if not.

	Stochastic	No-overfitting	Calc. speed	Data Size
LSF	$\times$	$\times$	$✓$	$\times$
DBI	$✓$	$✓$	$\times$	$✓$
MLE-MS	$\times$	$\times$	$\times$	$✓$
RPR	$\times$	$\times$	$\times$	$✓$
SMECF	$✓$	$\times$	$✓$	$\times$

Equations

$$\dot{x}=f(x;\mu ),$$

$$\dot{x}$$

$$\dot{y}$$

$$\dot{x}$$

$$\dot{y}$$

$$\begin{array}{rl}{\dot{\varphi }}_1& ={\omega }_1+{\epsilon }_1\sin ({\varphi }_2-{\varphi }_1)\\ {\dot{\varphi }}_2& ={\omega }_2+{\epsilon }_2\sin ({\varphi }_1-{\varphi }_2),\end{array}$$

$$\begin{array}{rl}\dot{\psi }& ={\Delta }_{\omega }+\epsilon q(\psi )=({\omega }_2-{\omega }_1)-({\epsilon }_1+{\epsilon }_2)\sin (\psi ).\end{array}$$

$$\begin{array}{r}{q}_1({\varphi }_1,{\varphi }_2)=Z({\varphi }_1)I({\varphi }_2).\end{array}$$

$$\begin{array}{r}{q}_1({\varphi }_1,{\varphi }_2)=\sin ({\varphi }_2-{\varphi }_1).\end{array}$$

$$\begin{array}{rl}\dot{{\varphi }_1}=& {\omega }_1\\ \dot{{\varphi }_2}=& {\omega }_2+q_2({\varphi }_1,{\varphi }_2)={\omega }_2+\cos ({\varphi }_1+\pi /2.5).\end{array}$$

$$\begin{split}\dot{x}_{1}&=\big{(}1-\sqrt{x_{1}^{2}+y_{1}^{2}}\big{)}x_{1}-\omega_{1}\,y_{1},\\ \dot{y}_{1}&=\big{(}1-\sqrt{x_{1}^{2}+y_{1}^{2}}\big{)}y_{1}+\omega_{1}\,x_{1},\end{split}$$

$$\begin{split}\dot{x}_{2}&=\big{(}1-\sqrt{x_{2}^{2}+y_{2}^{2}}\big{)}x_{2}-\omega_{2}\,y_{2}+\varepsilon x_{1}^{2},\\ \dot{y}_{2}&=\big{(}1-\sqrt{x_{2}^{2}+y_{2}^{2}}\big{)}y_{2}+\omega_{2}\,x_{2}.\end{split}$$

$$\dot{x}=f(x),$$

$${\dot{x}}_1$$

$${\dot{x}}_2$$

$$f(x)=Ax(1+|x|^{2})\mbox{ with }A=\left(\begin{array}[]{ccc}1&-1&1\\ 1&0&0\\ -4&2&-3\end{array}\right),$$

$${\dot{x}}_1=f(x_1)+\epsilon (x_2-x_1)$$

$${\dot{x}}_2=f(x_2)+\epsilon (x_1-x_2).$$

$${\dot{x}}_j=f_j(x_j)$$

$$f_{j}(x)=\omega_{j}Ax+(1-|x|^{2})x,\,\,\,\mbox{ with }A=\left(\begin{array}[]{cc}0&-1\\ 1&0\end{array}\right).$$

$$\dot{r}$$

$$\dot{\psi }$$

$$1-\epsilon +\epsilon \cos {\psi }^{\ast }<0.$$

$$1<\epsilon <\frac{4+{\Delta }_{\omega }^2}{4}.$$

$$x_1(t)=x_2(t)$$

$$y=\frac{1}{2}(x_1+x_2)\text{mbox}andz=\frac{1}{2}(x_1-x_2).$$

$$\dot{z}=J(t)z-2\epsilon z,$$

$$\dot{u}=J(t)u$$

$$\Vert z(t)\Vert \propto e^{(\lambda -2\epsilon )t},$$

$$\epsilon >\frac{\lambda }{2},$$

$${\epsilon }_c^1<\epsilon <{\epsilon }_c^2.$$

$$\dot{z}=-2\epsilon Hz.$$

$$\dot{V}(z)=z^T\dot{z}\le -2c\epsilon |z{|}^2.$$

Full text

1. I Introduction

2. I.1 Coupling functions, their nature and uses

3. I.2 Significance for interacting systems more generally

4. I.2.1 Physical effects of interactions: Synchronization, amplitude and oscillation death

5. I.2.2 Coupling strength and directionality

6. I.2.3 Coupling functions in general interactions

7. II Basic Concept of Coupling Functions

8. II.1 Principle meaning

9. II.1.1 Generic form of coupled systems

10. II.1.2 Coupling function definition

11. II.1.3 Example of coupling function and synchronization

12. II.2 History

13. II.3 Different domains and usage

14. II.3.1 Phase coupling functions

15. II.3.2 Amplitude coupling functions

16. II.3.3 Multivariate coupling functions

17. II.3.4 Generality of coupling functions

18. II.4 Coupling functions revealing mechanisms

19. II.5 Synchronization prediction with coupling functions

20. II.6 Unifying nomenclature

21. III Theory

22. III.1 Strong interaction

23. III.1.1 Two coupled oscillators

24. III.1.2 Comparison between approaches

25. III.2 Weak regime

26. III.2.1 Stable Periodic Orbit and its phase

27. III.2.2 Coupling function and phase reduction

28. III.2.3 Synchronization with external forcing

29. III.2.4 Phase response curve

30. III.2.5 Examples of the phase sensitivity function

31. III.3 Globally coupled oscillators

32. III.3.1 Coupling functions leading to multistability

33. III.3.2 Designing coupling functions for cluster states and chimeras

34. III.3.3 Coupling functions with delay

35. III.3.4 Low dimensional dynamics

36. III.3.5 Noise and nonautonomous effects

37. III.4 Networks of oscillators

38. III.4.1 Reduction to phase oscillators

39. III.4.2 Networks of chaotic oscillators

40. IV Methods

41. IV.1 Inferring coupling functions

42. IV.2 Methods for coupling function reconstruction

43. IV.2.1 Modeling by least-squares fitting

44. IV.2.2 Dynamical Bayesian inference

45. IV.2.3 Maximum likelihood estimation: multiple-shooting

46. IV.2.4 Random phase resetting method

47. IV.2.5 Stochastic modeling of effective coupling functions

48. IV.2.6 Comparison and overview of the methods

49. IV.3 Towards coupling function analysis

50. IV.4 Connections to other methodological concepts

51. IV.4.1 Phase reconstruction procedures

52. IV.4.2 Relation to phase response curve in experiments

53. IV.4.3 General effective connectivity modeling

54. V Applications and Experiments

55. V.1 Chemistry

56. V.2 Cardiorespiratory interactions

57. V.3 Neural coupling functions

58. V.4 Social sciences

59. V.5 Mechanical coupling functions

60. V.6 Secure communications

61. VI Outlook and Conclusion

62. VI.1 Future directions and open questions

63. VI.1.1 Theory

64. VI.1.2 Methods

65. VI.1.3 Analysis

66. VI.1.4 Integration theory-applications

67. VI.1.5 Applications

68. VI.2 Conclusion

I Introduction

I.1 Coupling functions, their nature and uses

Interacting dynamical systems abound in science and technology, with examples ranging from physics and chemistry, through biology and population dynamics, to communications and climate Kuramoto (1984); Winfree (1980); Pikovsky et al. (2001); Strogatz (2003); Haken (1983).

The interactions are defined by two main aspects: structure and function. The structural links determine the connections and communications between the systems, or the topology of a network. The functions are quite special from the dynamical systems viewpoint, as they define the laws by which the action and co-evolution of the systems are governed. The functional mechanisms can lead to a variety of qualitative changes in the systems. Depending on the coupling functions, the resultant dynamics can be quite intricate, manifesting a whole range of qualitatively different states, physical effects, phenomena and characteristics, including synchronization Pikovsky et al. (2001); Acebrón et al. (2005); Lehnertz and Elger (1998); Kapitaniak et al. (2012), oscillation and amplitude death Saxena et al. (2012); Koseska et al. (2013a), birth of oscillations Pogromsky et al. (1999); Smale (1976), breathers MacKay and Aubry (1994), coexisting phases Keller et al. (1992), fractal dimensions Aguirre et al. (2009), network dynamics Boccaletti et al. (2006); Arenas et al. (2008), and coupling strength and directionality Rosenblum and Pikovsky (2001); Hlaváčkováá-Schindler et al. (2007); Marwan et al. (2007); Stefanovska and Bračič (1999). Knowledge of such coupling function mechanisms can be used to detect, engineer or predict certain physical effects, to solve some man-made problems and, in living systems, to reveal their state of health and to investigate changes due to disease.

Coupling functions possess unique characteristics carrying implications that go beyond the collective dynamics (e.g. synchronization or oscillation death). In particular, the form of the coupling function can be used, not only to understand, but also to control and predict the interactions. Individual units can be relatively simple, but the nature of the coupling function can make their collective dynamics particular, enabling special behaviour. Additionally, there exist applications which depend just and only on the coupling functions, including examples of applications in social sciences and secure communication.

Given these properties, it is hardly surprising that coupling functions have recently attracted considerable attention within the scientific community. They have mediated applications, not only in different subfields of physics, but also beyond physics, predicated by the development of powerful methods enabling the reconstruction of coupling functions from measured data. The reconstruction within these methods is based on a variety of inference techniques, e.g. least squares and kernel smoothing fits Rosenblum and Pikovsky (2001); Kralemann et al. (2013a), 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 Galán et al. (2005); Levnajić and Pikovsky (2011); Timme (2007).

Both the connectivity between systems, and the associated methods employed for revealing it, are often differentiated into structural, functional and effective connectivity Park and Friston (2013); Friston (2011). Structural connectivity is defined by the existence of a physical link, like anatomical synaptic links in the brain or a conducting wire between electronic systems. Functional connectivity refers to the statistical dependences between systems, like for example correlation or coherence measures. Effective connectivity is defined as the influence one system exerts over another, under a particular model of causal dynamics. Importantly in this context, the methods used for the reconstruction of coupling functions belong to the group of effective connectivity techniques i.e. they exploit a model of differential equations and allow for dynamical mechanisms – like the coupling functions themselves – to be inferred from data.

Coupling function methods have been applied widely (Fig. 1), and to good effect: in chemistry, for understanding, effecting, or predicting interactions between oscillatory electrochemical reactions Kiss et al. (2007); Miyazaki and Kinoshita (2006); Tokuda et al. (2007); Blaha et al. (2011); Kori et al. (2014); in cardiorespiratory physiology Kralemann et al. (2013a); Stankovski et al. (2012); Iatsenko et al. (2013) for reconstruction of the human cardiorespiratory coupling function and phase resetting curve, for assessing cardiorespiratory time-variability and for studying the evolution of the cardiorespiratory coupling functions with age; in neuroscience for revealing the cross-frequency coupling functions between neural oscillations Stankovski et al. (2015); in social sciences for determining the function underlying the interactions between democracy and economic growth Ranganathan et al. (2014); for mechanical interactions between coupled metronomes Kralemann et al. (2008); and in secure communications where a new protocol was developed explicitly based on amplitude coupling functions Stankovski et al. (2014).

In parallel with their use to support experimental work, coupling functions are also at the centre of intense theoretical research Strogatz (2000); Daido (1996a); Crawford (1995); Acebrón et al. (2005). Particular choices of coupling functions can allow for a multiplicity of singular synchronized states Komarov and Pikovsky (2013). Coupling functions are responsible for the overall coherence in complex networks of non-identical oscillators Pereira et al. (2013); Ullner and Politi (2016); Luccioli and Politi (2010) and for the formation of waves and antiwaves in coupled neurons Urban and Ermentrout (2012). Coupling functions play important roles in the phenomena resulting from interaction such as synchronization Kuramoto (1984); Daido (1996a); Maia et al. (2015), amplitude and oscillation death Aronson et al. (1990); Koseska et al. (2013a); Zakharova et al. (2014); Schneider et al. (2015), the low-dimensional dynamics of ensembles Ott and Antonsen (2008); Watanabe and Strogatz (1993), and clustering in networks Ashwin and Timme (2005); Kori et al. (2014). The findings of these theoretical works are fostering further the development of methods for coupling function reconstruction, paving the way to additional applications.

I.2 Significance for interacting systems more generally

An interaction can result from a structural link through which causal information is exchanged between the system and one or more other systems Kuramoto (1984); Winfree (1980); Pikovsky et al. (2001); Strogatz (2003); Haken (1983). Often it is not so much the nature of the individual parts and systems, but how they interact, that determines their collective behaviour. One example is circadian rhythms, which occur across different scales and organisms DeWoskin et al. (2014). The systems themselves can be diverse in nature – for example, they can be either static or dynamical, including oscillatory, nonautonomous, chaotic, or stochastic characteristics Katok and Hasselblatt (1997); Kloeden and Rasmussen (2011); Landa (2013); Strogatz (2001); Suprunenko et al. (2013); Gardiner (2004). From the extensive set of possibilities, we focus in this review on dynamical systems, concentrating especially on nonlinear oscillators because of their particular interest and importance.

I.2.1 Physical effects of interactions: Synchronization, amplitude and oscillation death

An intriguing feature is that their mutual interactions can change the qualitative state of the systems. Thus they can cause transitions into or out of physical states such as synchronization, amplitude or oscillation death, or quasi-synchronized states in networks of oscillators.

The existence of a physical effect is, in essence, defined by the presence of a stable state for the coupled systems. Their stability is often probed through a dimensionally-reduced dynamics, for example the dynamics of their phase difference or of the driven system only. By determining the stability of the reduced dynamics, one can derive useful conclusions about the collective behaviour. In such cases, the coupling functions describe how the stable state is reached and the detailed conditions for the coupled systems to gain or lose stability. In data analysis, the existence of the physical effects is often assessed through measures that quantify – either directly or indirectly – the resultant statistical properties of the state that remains stable under interaction.

The physical effects often converge to a manifold, such as a limit cycle. Even after that, however, coupled dynamical systems can still exhibit their own individual dynamics, making them especially interesting objects for study.

Arguably, synchronization is the most studied of all such physical effects. It is defined as an adjustment of the rhythms of the oscillators, caused by their weak interaction Pikovsky et al. (2001). Synchronization is the underlying qualitative state that results from many cooperative interactions in nature. Examples include cardiorespiratory synchronization Schäfer et al. (1998); Stefanovska et al. (2000); Kenner et al. (1976), brain seizures Lehnertz and Elger (1998), neuromuscular activity Tass et al. (1998), chemistry Kiss et al. (2007); Miyazaki and Kinoshita (2006), the flashing of fireflies Mirollo and Strogatz (1990); Buck and Buck (1968) and ecological synchronization Blasius et al. (1999). Depending on the domain, the observable properties and the underlying phenomena, several different definitions and types of synchronization have been studied. These include phase synchronization, generalized synchronization, frequency synchronization, complete (identical) synchronization, lag synchronization and anomalous synchronization Kuramoto (1984); Brown and Kocarev (2000); Rulkov et al. (1995); Kocarev and Parlitz (1996); Pecora and Carroll (1990); Blasius et al. (2003); Pikovsky et al. (2001); Rosenblum et al. (1996); Ermentrout (1981); Arnhold et al. (1999); Eroglu et al. (2017).

Another important group of physical phenomena attributable to interactions are those associated with oscillation and amplitude deaths Bar-Eli (1985); Mirollo and Strogatz (1990); Prasad (2005); Suárez-Vargas et al. (2009); Zakharova et al. (2013); Schneider et al. (2015); Koseska et al. (2013a). Oscillation death is defined as a complete cessation of oscillation caused by the interactions, when an inhomogeneous steady state is reached. Similarly, in amplitude death, due to the interactions a homogeneous steady state is reached and the oscillations disappear. The mechanisms leading to these two oscillation quenching phenomena are mediated by different coupling functions and conditions of interaction, including strong coupling Mirollo and Strogatz (1990); Zhai et al. (2004), conjugate coupling Karnatak et al. (2007), nonlinear coupling Prasad et al. (2010), repulsive links Hens et al. (2013) and environmental coupling Resmi et al. (2011). These phenomena are mediated, not only by the phase dynamics of the interacting oscillators, but also by their amplitude dynamics, where the shear amplitude terms and the nonisochronicity play significant roles. Coupling functions define the mechanism through which the interaction causes the disappearance of the oscillations.

There is a large body of earlier work in which physical effects, qualitative states, or quantitative characteristics of the interactions were studied, where coupling functions constituted an integral part of the underlying interaction model, regardless of whether or not the term was used explicitly. Physical effects are very important and they are closely connected with the coupling functions. In such investigations, however, the coupling functions themselves were often not assessed, or considered as entities in their own right. In simple words, such investigations posed the question of whether physical effects occur; while for the coupling function investigations the question is rather how they occur. Our emphasis will therefore be on coupling functions as entities, on the exploration and assessment of different coupling functions, and on the consequences of the interactions.

I.2.2 Coupling strength and directionality

The coupling strength gives a quantitative measure of the information flow between the coupled systems. In an information-theoretic context, this is defined as the transfer of information between variables in a given process. In a theoretical treatment the coupling strength is clearly the scaling parameter of the coupling functions. There is great interest in being able to evaluate the coupling strength, for which many effective methods have been designed Smirnov and Bezruchko (2009); Chicharro and Andrzejak (2009); Staniek and Lehnertz (2008); Paluš and Stefanovska (2003); Bahraminasab et al. (2008); Jamšek et al. (2010); Sun et al. (2015); Faes et al. (2011); Marwan et al. (2007); Mormann et al. (2000); Rosenblum and Pikovsky (2001). The dominant direction of influence, i.e. the direction of the stronger coupling, corresponds to the directionality of the interactions. Earlier, it was impossible to detect the absolute value of the coupling strength, and a number of methods exist for detection only of the directionality through measurements of the relative magnitudes of the interactions – for example, when detecting mutual information Smirnov and Bezruchko (2009); Staniek and Lehnertz (2008); Paluš and Stefanovska (2003), but not the physical coupling strength. The assessment of the strength of the coupling 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 truly connected or not.

When the coupling function results from a number of functional components, its net strength is usually evaluated as the Euclidian norm of the individual components’ coupling strengths. Grouping the separate components, for example the Fourier components of periodic phase dynamics, one can evaluate the coupling strengths of the functional groups of interest. The latter could include the coupling strength from either one system or the other, or from both of them. Thus one can detect the strengths of the self, direct and common coupling components, or of the phase response curve Kralemann et al. (2011); Faes et al. (2015); Iatsenko et al. (2013). In a very similar way, these ideas can be generalized for multivariate coupling in networks of interacting systems.

It is worth noting that, when inferring couplings even from completely uncoupled or very weakly-coupled systems, the methods will usually detect non-zero coupling strengths. This results mainly from the statistical properties of the signals. Therefore, one needs to be able to ascertain whether the detected coupling strengths are genuine, or spurious, just resulting from the inference method. To overcome this difficulty, one can apply surrogate testing Schreiber and Schmitz (2000); Kreuz et al. (2004); Andrzejak et al. (2003); Paluš and Hoyer (1998) which generates independent, uncoupled, signals that have the same statistical properties as the original signals. The apparent coupling strength evaluated for the surrogate signals should then reflect a “zero-level” of apparent coupling for the uncoupled signals. By comparison, one can then assess whether the detected couplings are likely to be genuine. This surrogate testing process is also important for coupling function detection – one first needs to establish whether a coupling relation is genuine and then, if so, to try to infer the form of the coupling function.

I.2.3 Coupling functions in general interactions

The present review is focused mainly on coupling functions between interacting dynamical systems, and especially between oscillatory systems, because most studies to date have been developed in that context. However, interactions have also been studied in a broader sense for non-oscillatory, non-dynamical, systems, spread over many different fields, including for example quantum plasma interactions Marklund and Shukla (2006); Shukla and Eliasson (2011), solid state physics Higuchi and Yasuhara (2003); Farid (1997); Zhang (2013), interactions in semiconductor superlattices Bonilla and Grahn (2005), Josephson junction interactions Golubov et al. (2004), laser diagnostics Stepowski (1992), interactions in nuclear physics Mitchell et al. (2010); Guelfi et al. (2007), geophysics Murayama (1982), space science Feldstein (1992); Lifton et al. (2005), cosmology Faraoni et al. (2006); Baldi (2011), biochemistry Khramov and Bielawski (2007), plant science Doidy et al. (2012), oxygenation and pulmonary circulation Ward (2008), cerebral neuroscience Liao et al. (2013), immunology Robertson and Ghazal (2016), biomolecular systems Christen and Van Gunsteren (2008); Stamenović and Ingber (2009); Dong et al. (2014), gap junctions Wei et al. (2004) and protein interactions Jones and Thornton (1996); Teasdale and Jackson (1996); Okamoto et al. (2009); Gaballo et al. (2002). In many such cases, the interactions are different in nature. They are often structural, and not effective connections in the dynamics; or the corresponding coupling functions may not have been studied in this context before. Even though we do not discuss such systems directly in this review, many of the concepts and ideas that we introduce in connection with dynamical systems can also be useful for the investigation of interactions more generally.

II Basic Concept of Coupling Functions

II.1 Principle meaning

II.1.1 Generic form of coupled systems

The main problem of interest is to understand the dynamics of coupled systems from their building blocks. We start from the isolated dynamics:

[TABLE]

where $f:\mathbb{R}^{m}\times\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ is a differentiable vector field with $\mu$ being the set of parameters. For sake of simplicity, whenever there is no risk of confusion, we will omit the parameters. Over the last fifty years, developments in the theory of dynamical systems have illuminated the dynamics of isolated systems. For instance, we understand their bifurcations, including those that generate periodic orbits as well as those giving rise to chaotic motion. Hence we understand the dynamics of isolated systems in some detail.

In contrast, our main interest here is to understand the dynamics of the coupled equations:

[TABLE]

where $f_{1,2}$ are vector fields describing the isolated dynamics (perhaps with different dimensions) and $g_{1,2}$ are the coupling functions. The latter are our main objects of interest. We will assume that they are at least twice differentiable.

Note that we could also study this problem from an abstract point of view by representing the equations as:

[TABLE]

where the functions $q_{1,2}$ incorporate both the isolated dynamics and the coupling functions. This notation for inclusion of coupling functions, with no additive splitting between the interactions and the isolated dynamics, can sometimes be quite useful Aronson et al. (1990); Pereira et al. (2014). Examples include coupled cell networks Ashwin and Timme (2005), or the provision of full Fourier expansions Rosenblum and Pikovsky (2001); Kiss et al. (2005) when inferring coupling functions from data.

II.1.2 Coupling function definition

Coupling functions describe the physical rule specifying how the interactions occur. Being directly connected with the functional dependences, coupling functions focus not so much on whether there are interactions, but more on how they appear and develop. For instance, the magnitude of the phase coupling function affects directly the oscillatory frequency and describes how the oscillations are being accelerated or decelerated by the influence of the other oscillator. Similarly, if one considers the amplitude dynamics of interacting dynamical systems, the magnitude of the coupling function will prescribe how the amplitude is increased or decreased by the interaction.

A coupling function can be described in terms of its strength and form. While the strength is a relatively well-studied quantity, this is not true of the form. It is the functional form that has provided a new dimension and perspective, probing directly the mechanisms of interaction. In other words, the mechanism is defined by the functional form which, in turn, specifies the rule and process through which the input values are translated into output values i.e. in terms of one system (System A) it prescribes how the input influence from another system (System B) gets translated into consequences in the output of System A. In this way the coupling function can describe the qualitative transitions between distinct states of the systems e.g. routes into and out of synchronization. Decomposition of a coupling function provides a description of the functional contributions from each separate subsystem within the coupling relationship. Hence, the use of coupling functions amounts to much more than just a way of investigating correlations and statistical effects: it reveals the mechanisms underlying the functionality of the interactions.

II.1.3 Example of coupling function and synchronization

To illustrate the fundamental role of coupling functions in synchronization, we consider a simple example of two coupled phase oscillators Kuramoto (1984):

[TABLE]

where $\phi_{1},\phi_{2}$ are the phase variables of the oscillators, $\omega_{1},\omega_{2}$ are their natural frequencies, $\varepsilon_{1},\varepsilon_{2}$ are the coupling strength parameters, and the coupling functions of interest are both taken to be sinusoidal. (For further details including, in particular, the choice of the coupling functions, see also section III). Further, we consider coupling that depends only on the phase difference $\psi=\phi_{2}-\phi_{1}$. In this case, from $\dot{\psi}=\dot{\phi}_{2}-\dot{\phi}_{1}$ and Eqs. (5) we can express the interaction in terms of $\psi$ as:

[TABLE]

Synchronization will then occur if the phase difference $\psi$ is bounded, i.e. if Eq. (6) has at least one stable-unstable pair of solutions Kuramoto (1984). Depending on the form of the coupling function, in this case the sine form $q(\psi)=\sin(\psi)$, and on the specific parameter values, a solution may exist. For the coupling function given by Eq. (6) one can determine that the condition for synchronization to occur is $|\varepsilon_{1}+\varepsilon_{2}|\geq|\omega_{2}-\omega_{1}|$.

Fig. 2 illustrates schematically the connection between the coupling function and synchronization. An example of a synchronized state is sketched in Fig. 2(a). The resultant coupling strength $\varepsilon=(\varepsilon_{1}+\varepsilon_{2})$ has larger values of the frequency difference $\Delta_{\omega}=\omega_{2}-\omega_{1}$ at certain points within the oscillation cycle. As the condition $\dot{\psi}=0$ is fulfilled, there is a pair of stable and unstable equilibria, and synchronization exists between the oscillators. Fig. 2(b) shows the same functional form, but the oscillators are not synchronized because the frequency difference is larger than the resultant coupling strength. By comparing Figs. 2(a) and (b) one can note that while the form of the curve defined by the coupling function is the same in each case, the curve can be shifted up or down by choice of the frequency and coupling strength parameters. For certain critical parameters, the system undergoes a saddle-node bifurcation, leading to a stable synchronization.

The coupling functions of real systems are often more complex than the simple sine function presented in Fig. 2(a) and (b). For example, Fig. 2(c) also shows a synchronized state, but with an arbitrary form of coupling function that has two pairs of stable-unstable points. As a result, there could be two critical coupling strengths ($\varepsilon^{\prime}$ and $\varepsilon^{\prime\prime}$) and either one, or both, of them can be larger than the frequency difference $\omega_{2}-\omega_{1}$, leading to stable equilibria and fulfilling the synchronization condition. This complex situation could cause bistability (as will be presented below in relation to chemical experiments Sec. V.1). Thus comparison of Fig. 2(a) and (c) illustrates the fact that, within the synchronization state, there can be different mechanisms defined by different forms of coupling function.

II.2 History

The concepts of coupling functions, and of interactions more generally, had emerged as early as the first studies of the physical effects of interactions, such as the synchronization and oscillation death phenomena. In the seventeenth century, Christiaan Huygens observed and described the interaction phenomenon exhibited by two mechanical clocks Huygens (1673). He noticed that their pendula, which beat differently when the clocks were attached to a rigid wall, would synchronize themselves when the clocks were attached to a thin beam. He realised that the cause of the synchronization was the very small motion of the beam, and that its oscillations communicated some kind of motion to the clocks. In this way, Huygens described the physical notion of the coupling – the small motion of the beam which mediated the mutual motion (information flow) between the clocks that were fixed to it.

In the nineteenth century, John William Strutt, Lord Rayleigh, documented the first comprehensive theory of sound Rayleigh (1896). He observed and described the interaction of two organ pipes with holes distributed in a row. His peculiar observation was that for some cases the pipes could almost reduce one another to silence. He was thus observing the oscillation death phenomenon, as exemplified by the quenching of sound waves.

Theoretical investigations of oscillatory interactions emerged soon after the discovery of the triode generator in 1920 and the ensuing great interest in periodically alternating electrical currents. Appleton and van der Pol considered coupling in electronic systems and attributed it to the effect of synchronizing a generator with a weak external force Appleton (1922); Van Der Pol (1927). Other theoretical works on coupled nonlinear systems included studies of the synchronization of mechanically unbalanced vibrators and rotors Blekhman (1953), and the theory of general nonlinear oscillatory systems Malkin (1956). Further theoretical studies of coupled dynamical systems, explained phenomena ranging from biology, to laser physics, to chemistry Winfree (1967); Kuramoto (1975); Glass and Mackey (1979); Haken (1975); Wiener (1963). Two of these earlier theoretical works Winfree (1967); Kuramoto (1975) have particular importance and impact for the theory of coupling functions.

In his seminal work Winfree (1967) studied biological oscillations and population dynamics of limit-cycle oscillators theoretically. Notably, he considered the phase dynamics of interacting oscillators, where the coupling function was a product of two periodic functions of the form:

[TABLE]

Here, $I(\phi_{2})$ is the influence function through which the second oscillator affects the first, while the sensitivity function $Z(\phi_{1})$ describes how the first observed oscillator responds to the influence of the second one. (This was subsequently generalized for the whole population in terms of a mean field Winfree (1967, 1980)). Thus, the influence and sensitivity functions $I(\phi_{2})$, $Z(\phi_{1})$, as integral components of the coupling function, described the physical meaning of the separate roles within the interaction between the two oscillators. The special case $I(\phi_{2})=1+\cos(\phi_{2})$ and $Z(\phi_{1})=\sin(\phi_{1})$ has often been used Ariaratnam and Strogatz (2001); Winfree (1980).

Arguably, the most studied framework of coupled oscillators is the Kuramoto model. It was originally introduced in 1975 through a short conference paper Kuramoto (1975), followed by a more comprehensive description in an epoch-making book Kuramoto (1984). Today this model is the cornerstone for many studies and applications Acebrón et al. (2005); Strogatz (2000), including neuroscience Breakspear et al. (2010); Cumin and Unsworth (2007); Cabral et al. (2014), Josephson-junction arrays Wiesenfeld et al. (1996, 1998); Filatrella et al. (2000), power grids Dorfler and Bullo (2012); Filatrella et al. (2008), glassy states Iatsenko et al. (2014) and laser arrays Vladimirov et al. (2003). The model reduces the full oscillatory dynamics of the oscillators to their phase dynamics, i.e. to so-called phase oscillators, and it studies synchronization phenomena in a large population of such oscillators Kuramoto (1984). By setting out a mean-field description for the interactions, the model provides an exact analytic solution.

At a recent conference celebrating “40 years of the Kuramoto Model”, held at the Max Planck Institute for the Physics of Complex Systems, Dresden, Germany, Yoshiki Kuramoto presented his own views of how the model was developed, and described its path from initial ignorance on the part of the scientific community to dawning recognition followed by general acceptance: a video message is available Kuramoto (2015). Kuramoto devoted particular attention to the coupling function of his model, noting that:

In the year of 1974, I first came across Art Winfree’s famous paper [Winfree (1967)] …I was instantly fascinated by the first few paragraphs of the introductory section of the paper, and especially my interest was stimulated when he spoke of the analogy between synchronization transitions and phase transitions of ferroelectrics, […]. [There was a] problem that mutual coupling between two magnets (spins) and mutual coupling of oscillators are quite different. For magnetic spins the interaction energy is given by a scalar product of a two spin vectors, which means that in a particular case of planar spins the coupling function is given by a sinusoidal function of phase difference. In contrast, Winfree’s coupling function for two oscillators is given by a product of two periodic functions, […], and it seemed that this product form coupling was a main obstacle to mathematical analysis. […] I knew that product form coupling is more natural and realistic, but I preferred the sinusoidal form of coupling because my interest was in finding out a solvable model.

Kuramoto studied complex equations describing oscillatory chemical reactions Kuramoto and Tsuzuki (1975). In building his model, he considered phase dynamics and all-to-all diffusive coupling rather than local coupling, took the mean-field limit, introduced a random frequency distribution, and assumed that a limit-cycle orbit is strongly attractive Kuramoto (1975). As already mentioned, Kuramoto’s coupling function was a sinusoidal function of the phase difference:

[TABLE]

The use of the phase difference reduces the dimensionality of the two phases and provides a means whereby the synchronization state can be determined analytically in a more convenient way (see also Fig. 2).

The inference of coupling functions from data appeared much later than the theoretical models. The development of these methods was mostly dictated by the increasing accessibility and power of the available computers. One of the first methods for the extraction of coupling functions from data was effectively associated with detection of the directionality of coupling Rosenblum and Pikovsky (2001). Although directionality was the main focus, the method also included the reconstruction of functions that closely resemble coupling functions. Several other methods for coupling function extraction followed, including those by Kiss et al. (2005), Miyazaki and Kinoshita (2006), Tokuda et al. (2007), Kralemann et al. (2008), and Stankovski et al. (2012), and it remains a highly active field of research.

II.3 Different domains and usage

II.3.1 Phase coupling functions

A widely used approach for the study of the coupling functions between interacting oscillators is through their phase dynamics Kuramoto (1984); Winfree (1967); Pikovsky et al. (2001); Ermentrout (1986). If the system has a stable limit-cycle, one can apply phase reduction procedures (see Sec. III.2 for further theoretical details) which systematically approximate the high-dimensional dynamical equation of a perturbed limit cycle oscillator with a one-dimensional reduced-phase equation, with just a single phase variable $\phi$ representing the oscillator state Nakao (2015). In uncoupled or weakly-coupled contexts, the phases are associated with zero Lyapunov stability, which means that they are susceptible to tiny perturbations. In this case, one loses the amplitude dynamics, but gains simplicity in terms of the single dimension phase dynamics, which is often sufficient to treat certain effects of the interactions, e.g. phase synchronization. Thus phase connectivity is defined by the connection and influence between such phase systems.

To present the basic physics underlying a coupling function in the phase domain, we consider an elementary example of two phase oscillators that are unidirectionally phase-coupled:

[TABLE]

Our aim is to describe the effect of the coupling function $q_{2}(\phi_{1},\phi_{2})$ through which the first oscillator influences the second one. From the expression for $\dot{\phi_{2}}$ in Eq. (9) one can appreciate the fundamental role of the coupling function: $q_{2}(\phi_{1},\phi_{2})$ is added to the frequency $\omega_{2}$. Thus changes in the magnitude of $q_{2}(\phi_{1},\phi_{2})$ will contribute to the overall change of the frequency of the second oscillator. Hence, depending on the value of $q_{2}(\phi_{1},\phi_{2})$, the second oscillator will either accelerate or decelerate relative to its uncoupled motion.

The description of the phase coupling function is illustrated schematically in Fig. 3. Because in real situations one measures the amplitude state of signals, we explain how the amplitude signals (Fig. 3(a) and (d)) are affected depending on the specific phase coupling function (Fig. 3(b) and (c)). In all plots, time is scaled relative to the period $T_{1}$ of the amplitude of the signal originating from the first oscillator $x_{1}(t)$ (e.g. $x_{1}(t)=\sin(\phi_{1})$). For convenient visualisation of the effects we set the second oscillator to be fifteen times slower than the first oscillator: $\omega_{2}/\omega_{1}=15$. The particular coupling function $q_{2}(\phi_{1},\phi_{2})=\cos(\phi_{1}+\pi/2.5)$ presented on a $2\pi\times 2\pi$ grid (Fig. 3(b)) resembles a shifted cosine wave, which changes only along the $\phi_{1}$-axis, like a direct coupling component. Because all the changes occur along the $\phi_{1}$-axis, and for easier comparison, we also present in Fig. 3(c) a $\phi_{2}$-averaged projection of $q_{2}(\phi_{1},\phi_{2})$.

Finally, Fig. 3(d) shows how the second oscillator $x_{2}(t)$ is affected by the first oscillator in time in relation to the phase of the coupling function: when the coupling function $q_{2}(\phi_{1},\phi_{2})$ is increasing, the second oscillator $x_{2}(t)$ accelerates; similarly, when $q_{2}(\phi_{1},\phi_{2})$ decreases, $x_{2}(t)$ decelerates. Thus the form of the coupling function $q_{2}(\phi_{1},\phi_{2})$ shows in detail the mechanism through which the dynamics and the oscillations of the second oscillator are affected: in this case they were alternately accelerated or decelerated by the influence of the first oscillator.

Of course, coupling functions can in general be much more complex than the simple example presented ($\cos(\phi_{1}+\pi/2.5)$). This form of phase coupling function with a direct contribution (predominantly) only from the other oscillator is often found as a coupling component in real applications, as will be discussed below. Other characteristic phase coupling functions of that kind could include the coupling functions from the Kuramoto model (Eq. (8)) and the Winfree model (Eq. (7)), as shown in Fig. 4. The sinusoidal function of the phase difference from the Kuramoto model exhibits a diagonal form in Fig. 4(a), while the influence-sensitivity product function of Winfree model is given by a more complex form spread differently along the two-dimensional space in Fig. 4(b). Although these two functions differ from those in the previous example (Fig. 3), the procedure used for their interpretation is the same.

II.3.2 Amplitude coupling functions

Arguably, it is more natural to study amplitude dynamics than phase dynamics, as the former is directly observable while the phase needs to be derived. Real systems often suffer from the “curse of dimensionality” Keogh and Mueen (2011) in that not all of the features of a possible (hidden) higher-dimensional space are necessarily observable through the low-dimensional space of the measurements. Frequently, a delay embedding theorem Takens (1981) is used to reconstruct the multi-dimensional dynamical system from data. In real application with non-autonomous and non-stationary dynamics, however the theorem often does not give the desired result Clemson and Stefanovska (2014). Nevertheless, amplitude state interactions also have a wide range of applications both in theory and methods, especially in the cases of chaotic systems, strong couplings, delayed systems, and large nonlinearities, including cases where complete synchronization Cuomo and Oppenheim (1993); Stankovski et al. (2014); Kocarev and Parlitz (1995) and generalized synchronization Rulkov et al. (1995); Kocarev and Parlitz (1996); Stam et al. (2002); Abarbanel et al. (1993); Arnhold et al. (1999) has been assessed through observation of amplitude state space variables.

Amplitude coupling functions affect the interacting dynamics by increasing or decreasing the state variables. Thus amplitude connectivity is defined by the connection and influence between the amplitude dynamics of the systems. The form of the amplitude coupling function can often be a polynomial function or diffusive difference between the states.

To present the basics of amplitude coupling functions, we discuss a simple example of two interacting Poincaré limit-cycle oscillators. In the autonomous case, each of them is given by the polar (radial $r$ and angular $\phi$) coordinates as: $\dot{r}=r(1-r)$ and $\dot{\phi}=\omega$. In this way, a Poincaré oscillator is given by a circular limit-cycle and monotonically growing (isochronous) phase defined by the frequency parameter. In our example, we transform the polar variables to Cartesian (state space) coordinates $x=r\cos(\phi)$, $y=r\sin(\phi)$, and we set unidirectional coupling, such that the first (autonomous) oscillator:

[TABLE]

is influencing the $x_{2}$ state of the second oscillator through the quadratic coupling function $q_{2}(x_{1},y_{1},x_{2},y_{2})=x_{1}^{2}$:

[TABLE]

For simpler visual presentation we choose the first oscillator to be twenty times faster than the second one, i.e. their frequencies are in the ratio $\omega_{2}/\omega_{1}=20$, and we set a relatively high coupling strength $\varepsilon=5$.

The description of the amplitude coupling function is illustrated schematically in Fig. 5. In theory, the coupling function $q_{2}(x_{1},y_{1},x_{2},y_{2})$ has four variables, but for better visual illustration, and because the dependence is only on $x_{1}$, we show it only in respect to the two variables $x_{1}$ and $x_{2}$ i.e. $q_{2}(x_{1},x_{2})$. The form of the coupling function is quadratic, and it changes only along the $x_{1}$-axis, as shown in Figs. 5(b) and (c). Finally, Fig. 5(d) shows how the second oscillator $x_{2}(t)$ is affected by the first oscillator in time via the coupling function: when the quadratic coupling function $q_{2}(x_{1},x_{2})$ is increasing, the amplitude of the second oscillator $x_{2}(t)$ increases; similarly, when $q_{2}(x_{1},x_{2})$ decreases, $x_{2}(t)$ decreases as well.

The particular example chosen for presentation used a quadratic function $x_{1}^{2}$; other examples include a direct linear coupling function e.g. $x_{1}$, or a diffusive coupling e.g. $x_{2}-x_{1}$ Aronson et al. (1990); Kocarev and Parlitz (1996); Mirollo and Strogatz (1990). There are a number of methods which have inferred models that include amplitude coupling functions inherently Voss et al. (2004); Smelyanskiy et al. (2005); Friston (2002) or have pre-estimated most probable models Berger and Pericchi (1996), but without including explicit assessment of the coupling functions. Due to the multi-dimensionality and the lack of a general property in a dynamical system (like for example the periodicity in phase dynamics), there are countless possibilities for generalization of the coupling function. In a sense, this lack of general models is a deficiency in relation to the wider treatment of amplitude coupling functions. There are open questions here and much room for further work on generalising such models, in terms both of theory and methods, taking into account the amplitude properties of subgroups of dynamical systems, including for example the chaotic, oscillatory, or reaction-diffusion nature of the systems.

II.3.3 Multivariate coupling functions

Thus far, we have been discussing pairwise coupling functions between two systems. In general, when interactions occur between more than two dynamical systems, in a network (Sec. III.4), there may be multivariate coupling functions with more than two input variables. For example, a multivariate phase coupling function could be $q_{1}(\phi_{1},\phi_{2},\phi_{3})$, which is a triplet function of influence in the dynamics of the first phase oscillator caused by a common dependence on three other phase oscillators. Such joint functional dependences can appear as clusters of subnetworks within a network Albert and Barabási (2002).

Multivariate interactions have been the subject of much attention recently, especially in developing methods for detecting the couplings Baselli et al. (1997); Nawrath et al. (2010); Frenzel and Pompe (2007); Paluš and Vejmelka (2007); Kralemann et al. (2011); Duggento et al. (2012); Faes et al. (2011). This is particularly relevant in networks, where one can miss part of the interactions if only pairwise links are inferred, or a spurious pairwise link can be inferred as being independent when they are actually part of a multivariate joint function. In terms of networks and graph theory, the multivariate coupling functions relate to hypergraph, which is defined as a generalization of a graph where an edge (or connection) can connect any number of nodes (or vertices) Karypis and Kumar (2000); Zass and Shashua (2008); Weighill and Jacobson (2015).

Multivariate coupling functions have been studied by inference of small-scale networks where the structural coupling can differ from the inferred effective coupling Kralemann et al. (2011). The authors considered a network of three van der Pol oscillators where, in addition to pairwise couplings, there was also a joint multivariate cross-coupling function, for example of the form $\varepsilon x_{2}x_{3}$ in the dynamics of the first oscillator $\ddot{x}_{1}$. Due to the latter coupling, the effective phase coupling function is of a multivariate triplet nature. By extracting the phases and applying an inference method, the effective phase coupling was reconstructed, as illustrated by the example in Fig. 6. Comparing the true (Fig. 6 left) and the inferred effective (Fig. 6 right) diagrams, one can see that an additional pairwise link from the third to the first oscillator has been inferred. If the pairwise inference alone was being investigated one might conclude, wrongly, that this direct pairwise coupling was genuine and the only link – whereas in reality it is just an indirect effect from the actual joint multivariate coupling. In this way, the inference of multivariate coupling functions can provide a deeper insight into the connections in the network.

A corollary is the detection of triplet synchronization Kralemann et al. (2013b); Jia et al. (2015). This is a synchronization phenomenon which has an explicit multivariate coupling function of the form $q_{1}(\phi_{1},\phi_{2},\phi_{3})$ and which is tested in respect of the condition $|m\phi_{1}+n\phi_{2}+l\phi_{3}|\leq{\rm const}$, for $n,m,l$ negative or positive. It is shown that the state of triplet synchronization can exist, even though each pair of systems remains asynchronous.

The brain mediates many oscillations and interactions on different levels Park and Friston (2013). Interactions between oscillations in different frequency bands are referred to as cross-frequency coupling in neuroscience Jensen and Colgin (2007). Recently, neural cross-frequency coupling functions were extracted from multivariate networks Stankovski et al. (2015) (see also Sec. V.3). The network interactions between the five brainwave oscillations $\delta$, $\theta$, $\alpha$, $\beta$ and $\gamma$ were analysed by reconstruction of the multivariate phase dynamics, including the inference of triplet and quadruplet coupling functions. Fig. 7 shows a triplet coupling function of how the $\theta$ and $\alpha$ influence $\gamma$ brain oscillations. It was found that the influence from theta oscillations is greater than from alpha, and that there is significant acceleration of gamma oscillations when the theta phase cycle changes from $\pi$ to $2\pi$.

Very recently, Bick et al. (2016) have shown theoretically that symmetrically-coupled phase oscillators with multivariate (or non-pairwise) coupling functions can yield rich dynamics, including the emergence of chaos. This was observed even for as few as $N=4$ oscillators. In contrast to the Kuramoto-Sakaguchi equations, the additional multivariate coupling functions mean that one can find attracting chaos for a range of normal-form parameter values. Similarly, it was found that even the standard Kuramoto model can be chaotic with a finite number of oscillators Popovych et al. (2005).

II.3.4 Generality of coupling functions

The coupling function is well-defined from a theoretical perspective. That is, once we have the model (as in Eqs. 1), the coupling function is unique and fixed. The solutions of the equations also depend continuously on the coupling function. Small changes in the coupling function will cause only small changes in the solutions over finite time intervals. If solutions are attracted to some set exponentially and uniformly fast, then small changes in the coupling do not affect the stability of the system.

When we want to infer the coupling function from data we can face a number of challenges in obtaining a unique result (Sec. IV). Typically, we measure only projections of the coupling function, which might in itself lead to non-uniqueness of the estimate. That is, we project the function (which is infinite-dimensional) onto a finite-dimensional vector space. In doing so, we could lose some information and, generically, it is not possible to estimate the function uniquely (even without taking account of noise and perturbations). Furthermore, the final form of the estimated function will depend on the choice and number of base functions. For example, the choice of Fourier series or general orthogonal polynomials as base functions can affect slightly the final estimate of the coupling function. The choice of which base functions to be used is infinite. Even though many aspects of coupling functions (like the number of arguments, decomposition under an appropriate model, analysis of coupling function components, prediction with coupling functions, etc.), can be applied with great generality, the coupling functions themselves cannot be determined uniquely.

In the literature, authors often speak of the commonly-used coupling functions including, but not limited to, those listed in Table 1. Note that reactive and diffusive coupling have functionally the same form, the difference being that the reactive case includes complex amplitudes. This results in a phase difference between the coupling and the dynamics. Also in the literature, a diffusive coupling function $q(y-x)$ satisfying a local condition $q^{\prime}(0)<0$ is called dissipative coupling Rul’Kov et al. (1992). This condition resembles Fick’s law as the coupling forces the coupled system to converge towards the same state. When $q^{\prime}(0)>0$ the coupling is called repulsive Hens et al. (2013). Chemical synapses are an important form of coupling where the influences of $x$ and $y$ appear together as a product. There are also other interesting forms of coupling such as the geometric mean and further generalizations Petereit and Pikovsky (2017); Prasad et al. (2010). In environmental coupling, the function is given by the solution of a differential equation. In this case one can consider $\dot{y}=-\kappa y+\varepsilon(x(t)+y(t))$ for $\kappa>0$, so that the variables are considered as external fields driving the equation. Its solution $y(t)=y(t;x,y)$ is taken as the coupling function $q(x,y)$ and, for $t\gg 1$, is given in the table. The generality of coupling functions, and the fact that the form can come from an unbounded set of functions, were used to construct the encryption key in a secure communications protocol Stankovski et al. (2014) (see Sec. V.6).

II.4 Coupling functions revealing mechanisms

The functional form is a qualitative property that defines the mechanism and acts as an additional dimension to complement the quantitative characteristics such as the coupling strength, directionality, frequency parameter and limit-cycle shape parameters. By definition, the mechanism involves some kind of function or process leading to a change in the affected system. Its significance is that it may lead to qualitative transitions and induce or reduce physical effects, including synchronization, instability, amplitude death, or oscillation death.

But why is the mechanism important, and how it can be used? The first and foremost use of the coupling function mechanism is to illuminate the nature of the interactions themselves. For example, the coupling function of the Belousov-Zhabotinsky chemical oscillator has been reconstructed Miyazaki and Kinoshita (2006) with the help of a method for the inference of phase dynamics. Fig. 8 shows such a coupling function, demonstrating a form that is very far from a sinusoidal function: a curve that gradually decreases in the region of a small $\psi$ and abruptly increases at a larger $\psi$, with its minimum and maximum at around $5/4\pi$ and $7/4\pi$, respectively.

Another important set of examples is the class of coupling functions and phase response curves used in neuroscience. In neuronal interactions, some variables are very spike-like i.e. they resemble delta functions. Consequently, neuronal coupling functions (which are convolution of phase response curves and perturbation functions) then depend only, or mainly, on the phase response curves. So the interaction mechanism is defined by the phase response curves: quite a lot of work has been done in this direction Schultheiss et al. (2011); Tateno and Robinson (2007); Gouwens et al. (2010); Ermentrout (1996); see also Sec. IV.4.2. For example, Tateno and Robinson (2007) and Gouwens et al. (2010) reconstructed experimentally the phase response curves for different types of interneurons in rat cortex, in order to better understand the mechanisms of neural synchronization.

The mechanism of a coupling function depends on the differing contributions from individual oscillators. Changes in form may depend predominantly on only one of the phases (along one-axis), or they may depend on both phases, often resulting in a complicated and intuitively unclear dependance. The mechanism specified by the form of the coupling function can be used to distinguish the individual functional contributions to a coupling. One can decompose the net coupling function into components describing the self, direct and indirect couplings Iatsenko et al. (2013). The self-coupling describes the inner dynamics of an oscillator which results from the interactions and has little physical meaning. Direct-coupling describes the influence of the direct (unidirectional) driving that one oscillator exerts on the other. The last component, indirect-coupling, often called common-coupling, depends on the shared contributions of the two oscillators e.g. the diffusive coupling given with the phase difference terms. This functional coupling decomposition can be further generalized for multivariate coupling functions, where for example, a direct coupling from two oscillators to a third one can be determined Stankovski et al. (2015).

After learning the details of the reconstructed coupling function, one can use this knowledge to study or detect the physical effects of the interactions. In this way, the synchronous behavior of the two coupled Belousov-Zhabotinsky reactors can be explained in terms of the coupling function as illustrated by the examples given in Fig. 8 Miyazaki and Kinoshita (2006) and in Sec. V.1. Furthermore, the mechanisms and form of the coupling functions can be used to engineer and construct a particular complex dynamical structure, including sequential patterns and desynchronization of electrochemical oscillations Kiss et al. (2007). Even more importantly, one can use knowledge about the mechanism of the reconstructed coupling function to predict transitions of the physical effects – an important property described in detail for synchronization in the following section.

II.5 Synchronization prediction with coupling functions

Synchronization is a widespread phenomenon whose occurrence and disappearance can be of great importance. For example, epileptic seizures in the brain are associated with excessive synchronization between a large number of neurons, so there is a need to control synchronization to provide a means of stopping or preventing seizures Schindler et al. (2007); while in power grids the maintenance of synchronization is of crucial importance Rubido (2015). Therefore, one often needs to be able to control and predict the onset and disappearance of synchronization.

A seminal work on coupling functions by Kiss et al. (2005) uses the inferred knowledge of the coupling function to predict characteristic synchronization phenomena in electrochemical oscillators. In particular, the authors demonstrated the power of phase coupling functions, obtained from direct experiments on a single oscillator, to predict the dependence of synchronization characteristics such as order-disorder transitions on system parameters, both in small sets and in large populations of interacting electrochemical oscillators.

The authors investigated the parametric dependence of mutual entrainment using an electrochemical reaction system, the electrodissolution of nickel in sulfuric acid (see also Sec. V.1 for further applications on chemical coupling functions). A single nickel electrodissolution oscillator can have two main characteristic waveforms of periodic oscillation – the smooth type and the relaxation oscillation type. The phase response curve is of the smooth type and is nearly sinusoidal, while being more asymmetric for the relaxation oscillations.

The coupling functions are calculated using the phase response curve obtained from experimental data for the variable through which the oscillators are coupled. The coupling functions $q(\psi)$ of two coupled oscillators are reconstructed for three characteristic cases, as shown in Fig. 9(a)-(c), left panels. The right panels in Fig. 9 show the corresponding odd (antisymmetric) part of the coupling functions $q_{-}(\psi)=[q(\psi)-q(-\psi)]/2$, which is important for determination of the synchronization. The coupling functions $q(\psi)$ Fig. 9 (a)-(c) have predominantly positive values, so the interactions contribute to the acceleration of the affected oscillators. The first coupling function Fig. 9(a) for smooth oscillations has a sinusoidal $q_{-}(\psi)$ which can lead to in-phase synchronization at the phase difference of $\psi^{*}=0$. The third case of relaxation oscillations Fig. 9(c) has an inverted sinusoidal form $q_{-}(\psi)$, leading to stable anti-phase synchronization at $\psi^{*}=\pi$. The most peculiar case is the second one Fig. 9(b) of relaxation oscillations, where the odd coupling function $q_{-}(\psi)$ takes the form of a second harmonic ($q_{-}(\psi)\approx\sin(2\psi)$) and both the in-phase ($\psi^{*}=0$) and anti-phase ($\psi^{*}=\pi$) entrainments are stable, in which case the actual state attained will depend on the initial conditions.

Next, the knowledge obtained from experiments with a single oscillator was applied to predict the onset of synchronization in experiments with 64 globally coupled oscillators. The experiments confirmed that for smooth oscillators the interactions converge to a single cluster, and for relaxational oscillators they converge to a two-cluster synchronized state. Experiments in a parameter region between these states, in which bistability is predicted, are shown in Fig. 10. A small perturbation of the stable one-cluster state (left panel of Fig. 10) yields a stable two-cluster state (right panel of Fig. 10). Therefore, all the synchronization behavior seen in the experiments was in agreement with prior predictions based on the coupling functions.

In a separate line of work, synchronization was also predicted in neuroscience: interaction mechanisms involving individual neurons, usually in terms of phase-response curves (PRCs) or spike-time response-curves (STRCs), were used to understand and predict the synchronous behavior of networks of neurons Acker et al. (2003); Netoff et al. (2005); Schultheiss et al. (2011). For example, Netoff et al. (2005) studied experimentally the spike-time response-curves of individual neuronal cells. Results from these single-cell experiments were then used to predict the multi-cell network behaviors, which were found to be compatible with previous model-based predictions of how specific membrane mechanisms give rise to the empirically measured synchronization behavior.

II.6 Unifying nomenclature

Over the course of time, physicists have used a range of different terminology for coupling functions. For example, some publications refer to them as interaction functions and some as coupling functions. This inconsistency needs to be overcome by adopting a common nomenclature for the future.

The terms interaction function and coupling function have both been used to describe the physical and mathematical links between interacting dynamical systems. Of these, coupling function has been used about twice as often in the literature, including the most recent. The term coupling is closer to describing a connection between two systems, while the term interaction is more general. Coupling implies causality, whereas interaction does not necessarily do so. Often correlation and coherence are considered as signatures of interactions, while they do not necessarily imply the existence of couplings. We therefore propose that the terminology be unified, and the term coupling function be used henceforth to characterise the link between two dynamical systems whose interaction is also causal.

III Theory

In physics one is likely to examine stable static configurations whereas, in dynamical interaction between oscillators, solutions will converge to a subspace. For example, if two oscillators are in complete synchronization the subspace is called the synchronization manifold and corresponds to the case where the oscillators are in the same state for all time Pecora and Carroll (1990); Fujisaka and Yamada (1983). So, within the subspace, the oscillators have their own dynamics and finer information on the coupling function is needed.

The analytical techniques and methods needed to analyze the dynamics will depend on whether the coupling strength is strong or weak. Roughly speaking, in the strong coupling regime, we will have to tackle the fully-coupled oscillators whereas in the weak coupling we can reduce the analysis to lower-dimensional equations.

III.1 Strong interaction

To illustrate the main ideas and challenges of treating the case of strong interaction, while keeping technicalities to a minimum, we will first discuss the case of two coupled oscillators. These examples contain the main ideas and reveal the role of the coupling function and how it guides the system towards synchronization.

III.1.1 Two coupled oscillators

We start by illustrating the variety of dynamical phenomena that can be encountered and the role played by the coupling function in the strong coupling regime.

Diffusion driven oscillations

When two systems interact they may display oscillations solely because of the interaction. This is the nature of the problem posed by Smale (1976) based on Turing’s idea of morphogenesis Turing (1952). We consider two identical systems which, when isolated, each exhibit a globally asymptotically stable equilibrium, but which oscillate when diffusively coupled. This phenomenon is called diffusion driven oscillation.

Assume that the system

[TABLE]

where $f:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ is a differentiable vector field with a globally stable attraction with point – all trajectories will converge to this point. Now consider two of such systems coupled diffusively

[TABLE]

The problem proposed by Smale was to find (if possible) a coupling function (positive definite matrix) $H$ such that the diffusively coupled system undergoes a Hopf bifurcation. Loosely speaking, one may think of two cells that by themselves are inert but which, when they interact diffusively, become alive in a dynamical sense and start to oscillate.

Interestingly, the dimension of the uncoupled systems comes into play. Smale constructed an example in four dimensions. Pogromsky et al. (1999) constructed examples in three dimensions and also showed that, under suitable conditions, the minimum dimension for diffusive coupling to result in oscillation is $n=3$. The following example illustrates the main ideas. Consider

[TABLE]

where $|x|^{2}=x^{T}x$. Note that all the eigenvalues of $A$ have negative real parts. So the origin of the system Eq. (14) is exponentially attracting.

Consider the coupling function to be the identity

[TABLE]

For $\varepsilon=0$ the origin is globally attracting; the uniform attraction persists when $\varepsilon$ is very small, and so the origin is still globally attracting. However, for large values of the coupling $\varepsilon>0.6512$ the coupled systems exhibit oscillatory solutions (the origin has undergone a Hopf bifurcation).

Generalizations: In this example the coupling function was the identity. Pogromsky et al. (1999) discussed further coupling functions, such as coupling functions of rank two that generate diffusion-driven oscillators. Further oscillations in originally passive systems have been reported in spatially extended systems Gomez-Marin et al. (2007). In diffusively coupled membranes, collective oscillation in a group of nonoscillatory cells can also occur as a result of spatially inhomogeneous activation factor Ma and Yoshikawa (2009). These ideas of diffusion leading to chemical differentiation have also been observed experimentally and generalized by including heterogeneity in the model Tompkins et al. (2014).

Oscillation death

We now consider the opposite problem: Systems which when isolated exhibit oscillatory behaviour but which, when coupled diffusively, cease to oscillate and where the solutions converge to an equilibrium point.

As mentioned above in Sec. IB, this phenomenon is called oscillation death Bar-Eli (1985); Koseska et al. (2013a); Mirollo and Strogatz (1990); Ermentrout and Kopell (1990). To illustrate the essential features we consider a normal form of the Hopf bifurcation

[TABLE]

where

[TABLE]

So, each isolated system has a limit cycle of amplitude $|x|^{2}=1$ and a frequency of $\omega_{j}$. Note that the origin $x=0$ is an unstable equilibrium point. In oscillation death when the systems are coupled, the origin may become stable.

Focusing on diffusive coupling, again, the question concerns the nature of the coupling function. Aronson et al. (1990) remarked that the simplest coupling function to have the desired properties is the identity with strength $\varepsilon$. The equations have the same form as Eq. (13) with $H$ being the identity.

The effect can be better understood in terms of phase and amplitude variables. Let $r_{1}$, $r_{2}$ be the amplitudes and $\phi_{1},\phi_{2}$ the phases of $x_{1}$ and $x_{2}$, respectively. We consider $r_{1}=r_{2}=r$ which captures the main causes of the effect, as well as the phase difference $\psi=\phi_{1}-\phi_{2}$. Then the equations in these variables can be well approximated as

[TABLE]

The conditions for oscillation death are a stable fixed point at $r=0$ along with a stable fixed point for the phase dynamics. The above equations provide the main mechanism for oscillation death.

First, we can determine the stable fixed point for the phase dynamics, as illustrated in Fig. 2. There is a fixed point $\psi^{*}$ if $\varepsilon>\Delta_{\omega}/2$ and $\sin\psi^{*}=\Delta_{\omega}/(2\varepsilon)$. We will assume that $\Delta_{\omega}>2$ which implies that, when the fixed point $\psi^{*}$ exists, $\varepsilon>1$.

Next, we analyse the stability of the fixed point $r^{*}=0$. This is determined by the linear part of Eq. 15. Hence, the condition for stability is

[TABLE]

Using the equation for the fixed point we have $\cos\psi^{*}=\sqrt{1-[\Delta_{\omega}/(2\varepsilon)]^{2}}$. Replacing this in the stability condition we obtain $\varepsilon<(4+\Delta_{\omega}^{2})/8$. The analysis reveals that the system will exhibit oscillation death if the coupling is neither too weak nor too strong. Because we are assuming that the mismatch is large enough, $\Delta_{\omega}>2$, then there are minimum and maximum coupling strengths for oscillation death

[TABLE]

Within this range, there are no stable limit cycles: the only attracting point is the origin, and so the oscillations are dead.

The full equation is tackled in Aronson et al. (1990). The main principle is that the eigenvalues of the coupling function modify the original eigenvalues of the system and change their stability. It is possible to generalize these claims to coupling functions that are far from the identity Koseska et al. (2013a). The system may converge, not only to a single fixed point, but to many Koseska et al. (2013b).

Synchronization

One of the main roles of coupling functions is to facilitate collective dynamics. Consider the diffusively coupled oscillators described by Eq. (13). We say that the diagonal

[TABLE]

is the complete synchronization manifold Brown and Kocarev (2000). Note that the synchronization manifold is an invariant subspace of the equations of motion for all values of the coupling strength. Indeed, when the oscillators synchronize the coupling term vanishes. So, they will be synchronized for all future time. The main question is whether the synchronization manifold is attractive, that is, if the oscillators are not precisely synchronized will they converge towards synchronization? Similarly, if they are synchronized, and one perturbs the synchronization, will they return to synchronization?

Let us first consider the case where the coupling is the identity $H(x)=x$, and discuss the key mechanism for synchronization. Note that there are natural coordinates to analyze synchronization

[TABLE]

These coordinates have a natural meaning. If the system synchronizes, $z\rightarrow 0$ and $y\rightarrow s$ with $\dot{s}=f(s)$. Hence, we refer to $y$ as the coordinate parallel to the synchronization subspace $x_{1}=x_{2}$, and to $z$ as the coordinate transverse to the synchronization subspace, as illustrated in Fig. 11.

The synchronization analysis follows two steps: (i) Obtaining a governing equation for the modes $z$ transverse to the synchronization subspace; and (ii) using the coupling function to damp instabilities in the transverse modes.

**(i) Obtain an equation for **z. Let us assume that the initial disturbance of $z$ is small. Then we can obtain a linear equation for $z$ by neglecting the high order terms proportional to $|z|^{2}$. Noting that $\dot{z}=(\dot{x}_{1}-\dot{x}_{2})/2$, using Eqs. (13) for $x_{1}$ and $x_{2}$, and expanding $f$ in a Taylor series we obtain

[TABLE]

where $J(t)=Df(x_{2}(t))$ is the Jacobian evaluated along a solution of $x_{2}$.

(ii) Coupling function to provide damping. The term $-2\varepsilon z$ coming from the coupling now plays the role of a damping term. So we expect that the coupling will win the competition with $J$ and will force the solutions of $z$ to decay exponentially fast to zero. To see this, we observe that the first term

[TABLE]

depends on the dynamics of $x_{2}$ alone. Typically, $\|u(t)\|\propto e^{\lambda t}$ for $\lambda>0$. Now, to obtain a bound on the solution of Eq. (17), we consider the ansatz $z=ue^{-2\varepsilon t},$ and notice that differentiating $z$ we obtain Eq. (17). Hence

[TABLE]

from the growth behaviour of the disturbance $z$ we can also obtain the critical coupling strength to observe synchronization.

Critical Coupling: From this estimate, we can also obtain the critical coupling such that the solutions $z$ decay to zero. For coupling strengths

[TABLE]

the oscillators synchronize.

Meaning of $\lambda>0$: This corresponds to chaotic behaviour in the synchronization manifold. If $z\rightarrow 0$ then, $J(t)$ will be the Jacobian along a solution of $\dot{s}=f(s).$ So $\lambda$ depends on the dynamics on the synchronization manifold. If $\lambda>0$ and the solutions are bounded, the dynamics of the synchronized system is chaotic. Roughly speaking, $\lambda>0$ means that two nearby trajectories will diverge exponentially fast for small times and, because the solutions are bounded, they will subsequently come close together again. So, the coupled systems can synchronize even if the dynamics of the synchronized system is chaotic, as shown in Fig. 11 for the chaotic Lorenz attractor. The number $\lambda$ is the maximum Lyapunov exponent of the synchronization subspace.

There are intrinsic challenges associated with the analysis, and more when we attempt to generalize these ideas and also because of the nonlinearities that we neglected during the analysis.

General coupling functions: From a mathematical perspective the argument above worked because the identity commutes with all matrices. For other coupling functions, the argument above cannot be applied, and we encounter three possible scenarios:

i) The coupling function does not damp instabilities and the system never synchronizes Pecora and Carroll (1998); Boccaletti et al. (2002).

ii) The coupling function damps out instabilities only for a finite range of coupling strengths.

[TABLE]

For instance, this is the case for the Rössler system with coupling only in the first variable Huang et al. (2009).

iii) The coupling function damps instabilities and there is a single critical coupling $\varepsilon_{c}$. This is the case, when the coupling function eigenvalues have positive real parts. Pereira et al. (2014). 2. 2.

Local versus global results: In the above argument we have expanded the vector field in a Taylor series and obtained a linear equation to describe how the systems synchronize. This means that any claim on synchronization is local. It is still an open question how to obtain global results. 3. 3.

Nonlinear effects: We have neglected the nonlinear terms (the Taylor remainders), which can make synchronization unstable. Many researchers have observed this phenomenon through the bubbling transition Ashwin et al. (1994); Viana et al. (2005); Venkataramani et al. (1996), intermittent loss of synchronization Gauthier and Bienfang (1996); Yanchuk et al. (2001), and the riddling basin Heagy et al. (1994); Ashwin and Timme (2005).

To highlight the role of the coupling function and illustrate the above challenges, we will show how to obtain global results depending on the coupling function and discuss how local and global results are related.

Global argument: Assume that $H$ is a Hermitian positive definite matrix. The main idea is to turn the problem upside down. That is, we see the vector field as perturbing the coupling function. So, consider the system with only the coupling function and use the transverse coordinates

[TABLE]

Since $H$ is positive definite we obtain $-z^{T}Hz\leq-2c\varepsilon|z|^{2},$ where $c=c(H)$ is the smallest eigenvalue of $H$. The global stability of the system can be obtained by constructing a Lyapunov function $V$. The system will be stable if $V$ is positive and its derivative $\dot{V}$ is negative. This system admits a quadratic Lyapunov function $V(z)=\frac{1}{2}z^{T}z.$ Indeed, taking the derivative

[TABLE]

Hence all solutions of Eq. (19) will converge to zero exponentially fast. Next, consider the coupled system

[TABLE]

where by the mean value theorem we obtain

[TABLE]

Because we did not Taylor-expand the vector fields, the equation is globally valid. Assuming that the Jacobian is bounded by a constant $M_{f}>0$, we obtain $|J(t,z)|\leq M_{f}|z|$.

Computing again the Lyapunov function for the coupled system (including the vector fields) we obtain

[TABLE]

The system will synchronize if $\dot{V}$ is negative. So synchronization is attained if

[TABLE]

Again the critical coupling has the same form as before. The coupling function came into play via the constant $c$, and instead of $\lambda$ we have $M_{f}$. Typically, $M_{f}$ is much larger than $\lambda$. So, global bounds are not sharp. This conservative bound guarantees that the coupling function can damp all possible instabilities transverse to the synchronization manifold. Moreover, they are persistent under perturbation.

Local results: First, we Taylor-expand the system to obtain

[TABLE]

in just the same form as before. Note however that the trick we used previously, by defining $\dot{u}=Ju$, is no longer applicable. Indeed, we use the ansatz $z=ue^{-2\varepsilon Ht},$ to obtain $\dot{z}=-2\varepsilon Hz+e^{-2\varepsilon Ht}J(t)u,$ and, since $H$ and $J(t)$ do not commute,

[TABLE]

the ansatz cannot be used. Thus we need a better way forward. So in the same way as we calculated the expansion rate $\lambda$ for $J(t)$, we calculate the expansion rate for Eq. (22). Such Lyapunov exponents are very important in a variety of contexts. For us, it suffices to know that there are various ways to compute them Dieci and Van Vleck (2002); Pikovsky and Politi (2016). We calculate the Lyapunov exponent for each value of the coupling strength $\varepsilon$ to obtain a function

[TABLE]

This function is called the Master Stability Function (MSF). We will extract the synchronization properties from $\Lambda(\varepsilon)$. As we already discussed, the solutions of Eq. (19) will behave as

[TABLE]

Now note that $\Lambda(0)=\lambda>0$ (the expansion rate of the uncoupled equation), since we considered the case of chaotic oscillators. Because of our assumptions, we know that there is a $\varepsilon_{c}$ such that

[TABLE]

and for which $\Lambda(\varepsilon)$ will become negative, and $z$ will converge to zero.

Global versus local results. In the global analysis, the critical coupling depends on $M_{f}$, which is an upper bound for the Jacobian. This approach is rigorous and guarantees that all solution will synchronize. In the local analysis, we linearized the dynamics about the synchronization manifold and computed the Lyapunov exponent associated with the transverse coordinate $z$. The critical coupling was then obtained by analysing the sign of the Lyapunov exponent. Generically, $\varepsilon_{c}^{\tt msf}\ll M_{f}/2c$. The main reasoning is as follows. The Lyapunov exponents measure the mean instability whereas, in the global argument, we consider the worst possible instability. So the local method allows us to obtain a sharp estimate for the onset of synchronization.

The pitfalls of the local results. The main challenge of the local method lies in the intricacies of the theory of Lyapunov exponents Barreira and Pesin (2002); Pikovsky and Politi (2016). These can be discontinuous functions of the vector field. In other words, the nonlinear terms we threw away as Taylor remainders can make the Lyapunov exponent jump from negative to positive. Moreover, in the local case we cannot guarantee that all trajectories will be uniformly attracted to the synchronization manifold. In fact, for some initial conditions trajectories are attracted to the synchronization manifold, whereas nearby initial conditions are not. This phenomenon is called riddling Heagy et al. (1994).

III.1.2 Comparison between approaches

As discussed above, there is a dichotomy between global versus local results, and sharp bounds for critical coupling. These issues depend on the coupling function. Some coupling functions allow one to employ a given technique and thereby obtain global or local results.

First, we compare the two main techniques used in the literature, that is, Lyapunov functions (LFs) and the master stability function (MSF). For a generic coupling function, the LFs are unknown; but Lyapunov exponents can be estimated efficiently by numerical methods Dieci and Van Vleck (2002); Froyland et al. (2013); Ginelli et al. (2007).

Given additional information on the coupling, we can further compare the techniques. Note that the coupling function $H$ can be nonlinear. In this case, we consider the Jacobian

[TABLE]

Moreover, we say that $\Gamma$ belongs to the Lyapunov class if there are positive matrices $Q$ and $P$ such that

[TABLE]

Whenever the matrix $\Gamma$ is in the Lyapunov class we can construct the Lyapunov function algorithmically.

Table 2 reveals that the MSF method is very versatile. Although it may not encompass nonlinear perturbations it provides a framework to tackle a generic class of coupling functions Huang et al. (2009). In the theory of chaotic synchronization, therefore, this has been the preferred approach. However, it should be used with caution.

III.2 Weak regime

In the weak-coupling regime, the coupling strength is by definition insufficient to affect the amplitudes significantly; however the coupling can still cause the phases to adapt and adjust their dynamics Kuramoto (1984). Many of the phenomena observed in nature relate to the weak coupling regime.

Mathematical descriptions of coupled oscillators in terms of their phases offer two advantages: first, it reduces the dimension of the problem, and secondly, it can reveal principles of collective dynamics and other phenomena.

The theory for the weak coupling regime is well-developed. In the seventies and early eighties Winfree (1967, 1980) and Kuramoto (1975, 1984) developed the idea of asymptotic phase and phase reduction. Also starting from the seventies, the mathematical theory for such phase reduction was brought to completion in terms of normally hyperbolic invariant manifolds Wiggins (2013); Eldering (2013); Hirsch et al. (1977). Since then, the phase reduction theory Nakao (2015) has been significantly extended and generalized, for inclusion of phase reduction in the case of strongly perturbed oscillations, for stochastic treatment of interacting oscillators subject to noise of different kinds, for oscillating neuronal populations, and for spatiotemporal oscillations in reaction-diffusion systems Nakao et al. (2014); Goldobin et al. (2010); Teramae et al. (2009); Yoshimura and Arai (2008); Kurebayashi et al. (2013); Brown et al. (2004). The main ingredient in this approach is an attracting periodic orbit.

III.2.1 Stable Periodic Orbit and its phase

If the system in question has an exponentially stable periodic orbit, the theory guarantees the existence of the reduction and provides a method to obtain it. Thanks to the works of Izhikevich (2007); Hoppensteadt and Izhikevich (2012); Ermentrout (1996); Ermentrout et al. (2008); Rinzel and Ermentrout (1998); Ermentrout and Terman (2010) we now have a phase description for certain classes of neurons and we understand its limitations Smeal et al. (2010). The strategy is as follows. We assume that the system

[TABLE]

where $f:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$, has a uniformly exponentially attracting periodic orbit $\gamma$ with period $T$, that is, $\gamma(t+T)=\gamma(t)$.

The orbit is exponentially stable if the trajectories of the system approach it exponentially fast and the rate of convergence does not depend on the initial time or on initial conditions (for points sufficiently close to the orbit).

We can parameterize the orbit by its phase $\phi$, $\gamma(\phi+2\pi)=\gamma(\phi)$. We can also re-parameterize time such that that phase $\phi$ increases uniformly along the orbit $\gamma$. That is, the phase is uniform frequency equal to unity. By the chain rule we then have

[TABLE]

The key idea here is that weak coupling can adjust the rhythm of the phase dynamics. The goal is to obtain the phase reduction solely on the basis of information about the isolated system (the orbit $\gamma$). To this end we need to extend the phase $\varphi$ to a neighborhood of the orbit. The main ingredient necessary for the reduction of the problem to its phase dynamics is the concept of asymptotic phase Winfree (1967, 1980), which will provide us with the coupling function.

Asymptotic phase: Right now, the phase $\phi$ is defined only along the orbit $\gamma$. Our first step is to extend $\phi$ to a neighbourhood of $\gamma$. Since the periodic orbit is exponentially and uniformly attracting, it will attract an open neighbourhood of $\gamma$. We call this set the basin of attraction of the periodic orbit. Note that every initial point ${x}_{0}$ in the basin of attraction of the orbit will converge to the orbit. Hence, we have a $\phi(x_{0})$ such that

[TABLE]

where ${x}(t,{x}_{0})$ is the solution of the system with initial condition ${x}_{0}$. For each initial point in the basin of attraction of $\gamma$ we can assign a unique point in the orbit $\theta$. This $\phi\in[0,2\pi]$ is called the asymptotic phase.

Isochron. For each value of phase $\phi$ in the orbit $\gamma$ we have a curve passing through this phase values. And along this curve every initial will have the same asymptotic phase. This set is called isochron. That is, the isochron is a level set of $\phi({x})$. So points in the isochron have the same value of phase and will move at the same speed. See Fig. 12 where points in the same isochron approach the orbit along the same isochron. The set of points where the isochron cannot be defined is called phaseless. Once we find the isochron we can perform the phase reduction.

III.2.2 Coupling function and phase reduction

Consider Eq. (23) with a stable periodic orbit $\gamma$ being perturbed

[TABLE]

where $\vartheta=\omega t$ is the phase of the external influence, and $I$ the influence is a periodic on $\vartheta$. The weak coupling implies $\varepsilon\ll 1$.

One of the cornerstones of the theory of invariant manifolds is to guarantee that, when the system is perturbed and the coupling strength is weak $\varepsilon\ll 1$, there is a new attracting periodic orbit $\tilde{\gamma}$ close to the orbit $\gamma$, the difference between the orbits being of order $\varepsilon$. Moreover, $\tilde{\gamma}$ is exponentially attractive and the isochrons also persist. So, while the amplitudes are basically unaffected, the dynamics of the phases change greatly.

With the help of the asymptotic phase, we define the phase in a neighborhood of the orbit $\gamma$. This neighborhood contains the new orbit $\tilde{\gamma}$ as we consider small $\varepsilon>0$. So, calculating the phase along $\tilde{\gamma}$, by the chain rule we obtain

[TABLE]

But by construction $\nabla\phi\cdot f=1$ in a neighborhood of $\gamma$. Because $\tilde{\gamma}$ is $O(\varepsilon)$ distant from $\gamma$ we can expand both $f$ and $I$ and evaluate them along $\gamma$ at the expense of a perturbation of order $\varepsilon$. It is standard to denote $Z=\nabla_{\gamma}\phi$. In this setting we have

[TABLE]

and we have successfully reduced the problem to the phase of the unperturbed orbit $\gamma$. In general terms, we study problems of the following type

[TABLE]

The main insight was to obtain the coupling function in terms of how the phase of the unperturbed orbit behaves near the orbit $\gamma$. We performed the following steps:

(1) Phase sensitivity $Z$ of the unperturbed system.

Once we have the asymptotic phase, we can use it as the new phase variable $\phi$, extending the definition of phase along the orbit to a neighborhood of the orbit. From the phase, we can in turn compute the phase sensitivity function

[TABLE]

where the gradient is evaluated along the orbit $\gamma$.

(2) Obtain coupling function by $q=Z\cdot I$.

For this step, we need to take the inner product of $Z$ with the perturbation $p$. When studying collective phenomena $q$ will contain fast variables and slow variables. Typically, only the slow variable are of interest, so we will average $q$ over the fast variables.

Meaning of $q$. In this approach we have a strong underlying assumption: that the phase responds linearly to perturbations. That is, the coupling function is linear in the perturbations. If the phase is perturbed by $I_{1}$ and $I_{2}$ the net effect will be the sum of $I_{1}$ and $I_{2}$. Notice, that the linearity is only in terms of the perturbations. The equation itself is nonlinear in the phase variable $\varphi$. The linearity with respect to perturbations is because we have discarded all nonlinear terms and terms of order $\varepsilon^{2}$ (by computing $Z$ along the unperturbed orbit). We will discuss these issues in an example below. This linearity will facilitate the study of networks and large ensembles of oscillators.

These two steps will provide the phase description for weakly coupled oscillators. Using these steps, it is possible to explain the collective behaviour of neurons Ermentrout (1996), and circadian dynamics Winfree (1980), among other processes Kuramoto (1984); Ermentrout and Terman (2010).

III.2.3 Synchronization with external forcing

We will illustrate and discuss how the above ideas can be applied to study the problem of synchronization with external forcing. Consider the system

[TABLE]

By inspection, it is clear that $x=0$ is an unstable point and the system has an attracting periodic orbit $\gamma$ of radius 1. This can be better seen by changing to polar coordinates $x_{1}(r,\phi)=r\sin(\varphi)$ and $x_{2}(r,\varphi)=r\cos(\varphi)$ using which, we obtain

[TABLE]

The orbit $\gamma$ corresponds to $r=1$ and is shown in Fig. 12.

Asymptotic phase: The phase $\varphi$ as defined in the orbit $\gamma$ has a constant frequency equal to unity. Along the orbit ${\gamma}$, we therefore have $\dot{\varphi}=1$ (by inspection of the equations). For points outside the orbit, however, this is no longer true. The asymptotic phase $\phi$ will fix this issue because the points then move at the same speed as the corresponding points in the orbit, so that $\dot{\phi}=1$ for points outside the orbit.

Because of the symmetry ($r$ does not depend on $\varphi$) we can use the ansatz

[TABLE]

where we aim to find the function $\zeta$. Differentiating we obtain $\dot{\phi}=\dot{\varphi}+\frac{d\zeta}{dr}\frac{dr}{dt}$ and, using the isochron’s properties together with the equations for $r$ and $\varphi$, we obtain $\dot{\zeta}=1/r$ so $\zeta=\log r+C.$ Since we want to extend the phase continuously from the orbit, if ${x}\in{\gamma}$ then $\phi(x)=\varphi(x)$. We choose the constant $C=0.$ Therefore,

[TABLE]

and we can define the isochron with asymptotic phase $\phi(r,\varphi)=c.$ In Fig. 12(top) we show four level sets of the asymptotic phase corresponding to $\phi=0,\pi/2,\pi,$ and $3\pi/2$.

We can use the asymptotic phase to obtain a coordinate that decouples the phase dynamics from the other coordinates. Note that, by defining a new coordinate $\phi=\varphi-\eta(r),$ we obtain

[TABLE]

which is valid, not only along the orbit $\gamma$ via Eq. 26, but also in a neighborhood of the orbit. In the first equality we just stressed the identity between the frequency (applying the chain rule) and the gradient of $\phi$.

We can now readily take the gradient (in polar coordinates), yielding $\nabla\phi=e_{1}(\sin(\phi-\log(r))+\cos(\phi-\log(r))/r.$ Along the unperturbed orbit ${\gamma}$ we have $Z(\phi)=\nabla_{\gamma}\phi$ so that the phase sensitivity functions are

[TABLE]

Next we obtain the coupling function.

Obtaining the coupling function of external forcing. Now we consider the system being forced at frequency:

[TABLE]

where $\vartheta=\omega t$. We obtain the coupling function through the isochron. We now justify in detail why we have discarded the corrections in $\varepsilon^{2}$.

We compute the equation for the phase dynamics. Note that, by chain rule $\dot{\phi}=\nabla\phi\cdot{\dot{x}}=\nabla\phi\cdot(f(x)+\varepsilon{I})$. Using Eq. (27) and evaluating the gradient along the orbit $\tilde{\gamma}$, we obtain

[TABLE]

For small $\varepsilon$, we know that the difference between $\nabla_{\gamma}\phi$ and $\nabla_{\tilde{\gamma}}\phi$ is of order $\varepsilon$, so we can replace the gradient along the perturbed orbit and unperturbed orbit with corrections of order $\varepsilon^{2}$ (because $\varepsilon$ is already multiplying the function). Hence,

[TABLE]

Synchronization and coupling function. The main idea is that the coupling function $q$ can help in adjusting the frequency of the system to the frequency $\omega$ of the forcing. As we discussed above, we will neglect the terms ${O}(\varepsilon^{2})$. Introducing the phase difference

[TABLE]

and considering $1-\omega=\Delta_{\omega}$ we obtain

[TABLE]

If $\Delta_{\omega}$ is of order $\varepsilon$ then the dynamics of $\psi$ will be slow in comparison with the dynamics of $\theta$. Roughly speaking, for each cycle of $\psi$ we have $1/\varepsilon$ cycles of $\theta$. Because, the dynamics of $\vartheta$ is faster than that of $\psi$, we use the averaging method to obtain the coupling function

[TABLE]

where $T=2\pi$ is the period of $p$ as a function of $\vartheta$. Note that, for our result, $Z$ is sinusoidal so that by integrating over $\vartheta$ while keeping $\psi$ fixed, we obtain $q(\psi)=A\sin(\psi+\beta)$. Hence we obtain the dynamics in terms of the phase difference

[TABLE]

which is exactly the equation shown in Fig. 2.

We are now ready to study collective phenomena between the driving and the system. For instance, the system will phase-lock with the driving when $\Delta+\varepsilon q(\psi^{*})=0$. In this case, the oscillators will have the same frequency. Moreover, because $q$ is a periodic function, the fixed point $\psi^{*}$ will exist only when $|\Delta/\varepsilon|\leq\max q$.

Higher order n:m phase locking. Our assumption is that $\Delta=O(\varepsilon)$, so that Eq. (29) for the phase difference $\psi=\phi-\vartheta$ is a slow variable. It may happen that $\psi_{mn}=m\varphi-n\vartheta$ gives rise to a slow variable. In such cases, we perform the same analysis for $\psi_{mn}$ and further information on the higher order phase locking can be obtained Ermentrout (1981).

III.2.4 Phase response curve

The phase sensitivity function $Z$ plays a major role in this analysis. It also has many names: infinitesimal phase response curve (iPRC), linear response function, infinitesimal phase resetting curve. It is deeply related to the so-called phase response curve or phase resetting curve (PRC). For an oscillator to be able to adjust its rhythm and synchronize, it must respond differently to the perturbations at different phases $\phi$. So the phase can advance or retard to adjust its rhythm to the external forcing. The PRC is a natural way of displaying the response of oscillators to perturbations and thereby to gain insight into the collective dynamics.

The main idea of the PRC is as follows: If we perform a small and short perturbation of the orbit, the phase may complete its cycle before expectation (in the absence of perturbations), or it may be delayed. The unperturbed period of the orbit $\gamma$ is $T_{0}$. Every point on the orbit can be uniquely described by the phase $\phi$. A small perturbation applied at a phase $\phi_{s}$ can cause the phase to complete its full cycle at a time $T_{1}$. The normalized phase difference between the cycles is

[TABLE]

Note that the PRC depends on the phase $\phi_{s}$ at which the small perturbation was applied, that is, $PRC=PRC(\phi_{s})$. This is the so-called phase response curve.

In the theory of weakly coupled oscillators, we use the concept of an infinitesimal PRC (iPRC). It is equivalent to the gradient of the phase $Z$, and it is defined as the PRC normalized by the amplitude of the perturbation $A$

[TABLE]

Indeed, the isochrons and the PRC are closely related. If a point moving along the orbit $\gamma$ is instantaneous and the perturbation is small, the point will land on an isochron, which tells us the new phase $\phi$ of the point once it comes back to the orbit. Further considerations and the relationship of the PRC to experiments are given in Sec. IV.4.2.

In neuroscience, pulse-coupled oscillators are an important class of models. Here the interactions happen in instantaneous pulses of communication. The collective dynamics of such models are of great interest Mirollo and Strogatz (1990). The relationship between pulse-coupled oscillators and the phase reduction presented above has recently been elucidated by Politi and Rosenblum (2015) who showed that the models are equivalent.

III.2.5 Examples of the phase sensitivity function

Because of the works of Winfree (1967, 1980), Kuramoto (1984), Ermentrout (1996); Stiefel et al. (2008), we now have a good understanding of the phase sensitivity $Z$ for many classes of systems such as heartbeats, circadian rhythms, and in some neurons (with stable repetitive firing).

The iPRC and PRC are closely related to the bifurcation that led to the oscillatory behaviour Ermentrout (1996); Brown et al. (2004). $Z(\phi)$ is a vector and, in our example in Sec. III.2.3, the norm of $Z$ was proportional to $\sin(\phi+\beta)$. This is typical of Hopf bifurcations. In Fig. 13 we present typical bifurcations in neuronal models for which the iPRC and PRC are relevant.

In neuron models, the coupling is in one single variable: the membrane potential $V$. So we only need to compute the derivative with respect to $V$. Thus, $Z(\phi)=\partial\phi/\partial V$ 111We are abusing the notation by using $Z$ to represent both the full gradient and the derivative with respect to a single variable. In theoretical neuroscience this convention is standard.. Izhikevich (2000) derived a phase model for weakly coupled relaxation oscillators and burster neurons Izhikevich (2007). Brown et al. (2004) obtained the phase sensitivity $Z$ for other interesting cases, including homoclinic oscillators. Neurons with stable repetitive firing (corresponding to a stable orbit) can be classified as having PRC type I dynamics corresponding to a SNIPER bifurcation, or PRC type II dynamics corresponding to a Hopf bifurcation. The phase portraits for these two bifurcations are illustrated in Fig. 13. PRCs of type I are always positive whereas PRCs of type II have both negative and positive parts, as shown in Table 3. The PRC type is indicative of the neuron’s ability to synchronize: networks of neurons with PRC type II can synchronize via mutual excitatory coupling, but those of PRC type I cannot Ermentrout (1996).

III.3 Globally coupled oscillators

Now suppose that we have $N$ coupled oscillators

[TABLE]

We assume that, when they are uncoupled $\varepsilon=0$, each system has an exponentially attracting periodic orbit. So the dynamics of the uncoupled system occurs on a torus $\mathbb{T}^{N}$ that is exponentially attracting. Moreover, we also assume that $f_{i}$ is close to $f$. As we turn the coupling on, the dynamics changes. The theory of normally hyperbolic invariant manifolds guarantees that the dynamics of system with small coupling will also take place on a torus Eldering (2013). So the amplitudes remain roughly the same. But the dynamics on the torus, that is, the phases can change a lot Turaev (2015).

Again if we know the isochrons for the phases we can use the same arguments to describe the system in terms of the phases. The corresponding phase model

[TABLE]

where each oscillator has its own period $T_{i}$, and $q_{ij}$ is the coupling function describing the influence of the $j$-th oscillator on the $i$-th oscillator. Here

[TABLE]

Note that $Z$ is independent of the index $i$ because we assumed that $f_{i}$’s are all close to $f$. Again, we can average over the fast variables to obtain equations in terms of the phase difference Daido (1996b).

In many cases, the isolated oscillators are close to a Hopf bifurcation. As discussed above, the coupling function after averaging takes the form

[TABLE]

This is by far the best-studied coupling function, and it has offered deep insights into collective properties for both globally coupled oscillators Acebrón et al. (2005) and complex networks Arenas et al. (2008); Rodrigues et al. (2016).

First, we consider $\beta=0$. The model is then written as

[TABLE]

If the oscillators are identical $\omega_{i}=\omega$ then any small coupling $\varepsilon>0$ leads to synchronization (the phases will converge to the same value). If the distribution $g$ of natural frequencies $\omega_{i}$ is broad, then at a critical coupling $\varepsilon_{c}$ a large cluster of synchronized oscillators appears and, with further increase of coupling, additional oscillators join the cluster.

The main idea of the analysis is to introduce an order parameter

[TABLE]

and to rewrite the equations in terms of the parameters $r$ and $\psi$ (which are now mean-field parameters)

[TABLE]

Taking the limit $N\rightarrow\infty$, we can write the model in terms of self-consistent equations Kuramoto (1984).

When the distribution $g$ of the natural frequencies is an even, unimodal, and non-increasing function, and the coupling is weak, the incoherent state is neutral Strogatz and Mirollo (1991), but the order parameter $r$ vanishes (at a polynomial rate) if $g$ is smooth Fernandez et al. (2014). Moreover, on increasing the coupling, the incoherent solution $r=0$ bifurcates for

[TABLE]

A systematic review of the critical coupling, including bi-modal distributions (keeping the coupling function purely harmonic) has been given by Acebrón et al. (2005).

III.3.1 Coupling functions leading to multistability

As we have seen, if the oscillators are close to a Hopf bifurcation as is typical, then the corresponding phase sensitivity $Z(\phi)\propto\sin(\phi+\beta)$. This means that the coupling function will have a phase shift

[TABLE]

This coupling function is called the Sakaguchi-Kuramoto coupling Sakaguchi and Kuramoto (1986). The slight modification can lead to non-monotonic behaviour of synchronization Omel’chenko and Wolfrum (2012). For certain unimodal frequency distributions $g$, the order parameter can decay as the coupling increases above the critical coupling, and the incoherent state can regain stability. Likewise multistability between partially synchronized states and/or the incoherent state can also appear.

Although the dynamics of the model with this slight modification can be intricate, it is still possible to treat a more general version of the Sakaguchi-Kuramoto coupling

[TABLE]

This coupling function generalizes the standard Sakaguchi-Kuramoto model as it allows for different contributions of oscillators to the mean field, on account of the phase shifts $\alpha_{i}$ and $\beta_{j}$ and coupling factors $B_{j}$. In turn, the mean field acts on each oscillator differently. This scenario is tractable in terms of the self-consistency equations for the amplitude and frequency of the mean field Vlasov et al. (2014). Also in this setting, solutions of the coupled phase oscillators approximate solutions of phase oscillators with an inertial term Dorfler and Bullo (2012) which plays a major role in power-grids.

Higher harmonics. In the previous discussion, the coupling function $q$ contained one harmonic. $q(\phi)=\sin(\phi+\beta)$. Depending on the underlying bifurcation we must now include further, higher-order, Fourier components,

[TABLE]

where $a_{i}$ and $b_{i}$ are parameters. For example, synchronization of weakly-coupled Hodgkin-Huxley neurons can be replicated using coupling functions consisting of the first four Fourier components Hansel et al. (1993b).

Moreover, considering the coupling function $q$ to be a biharmonic coupling function Hansel et al. (1993a), there is a multiplicity of such states, which differ microscopically in the distributions of locked phases Komarov and Pikovsky (2013). Higher harmonics in the coupling function can also lead to the onset of chaotic fluctuations in the order parameter Bick et al. (2011). Indeed, the coupling function alone can generate chaos, that is, even keeping the frequencies identical and having no amplitude variations.

III.3.2 Designing coupling functions for cluster states and chimeras

Ashwin and co-workers tailored the coupling function to obtain clusters states in identically and globally coupled phase oscillators Ismail and Ashwin (2015); Orosz et al. (2009). In this situation, because the oscillators are identical, the behaviour of the system is determined by the number of oscillators $N$ and the coupling function $q$. By carefully choosing the Fourier coefficients of the coupling function we obtain two major results: (a) any clustering can appear and be stable, and (b) open sets of coupling functions can generate heteroclinic network attractors. Heteroclinic networks are not confined to globally-coupled oscillators and they can appear robustly in complex networks Aguiar et al. (2011); Field (2015).

In networks of identical oscillators a chimera state is defined as a spatio-temporal pattern in which synchronous and asynchronous oscillations coexist Abrams and Strogatz (2004); Hagerstrom et al. (2012). Chimera states among phase oscillators appear only when a spatial (long-range) coupling is included. The approach described above, in which the coupling is tailored to obtain clusters and complicated attractors in identical and globally-coupled oscillators, can be extrapolated to construct chimeras. They can be obtained either by consideration of higher harmonics Ashwin and Burylko (2015), or by perturbation of the coupling function in specific ways. A partially coherent inhomogeneous pattern called chimera death, which combines the features of chimera states and oscillation death, has been also established Zakharova et al. (2014).

III.3.3 Coupling functions with delay

A natural generalization of the coupling function is to introduce delay. A common case is the inclusion of transmission delays

[TABLE]

The addition of delays makes the model infinite-dimensional and leads to a series of new phenomena such as bistability between synchronized and incoherent states, and unsteady solutions with time-dependent order parameters Yeung and Strogatz (1999). Multistability is very common in the presence of delays. In particular, it can also be observed in small variations of the previous coupling function

[TABLE]

and multistablity can also be observed Kim et al. (1997).

Carefully chosen communication delays can also be used to encode patterns in the temporal coding by spikes. These patterns can be obtained by a modulation of the multiple, coexisting, stable, in-phase synchronized states, or traveling waves propagating along or against the direction of coupling Popovych et al. (2011). Coupling functions with delay can also be used for controlling the state of oscillation death Zakharova et al. (2013). Two limiting cases of delay can be treated. First, for very small delays, the theory of an invariant manifold can be applied. Secondly, in the case of large delays, developments due to Lichtner et al. (2011) and Flunkert et al. (2010) can be used to determine the collective properties of ensembles of oscillators.

III.3.4 Low dimensional dynamics

A particularly striking observation is the low-dimensional dynamics of identical globally-coupled phase oscillators under the Sakaguchi-Kuramoto coupling function. Note that, in this case, writing the sinusoidal coupling in exponential form, we can express the coupled equations as

[TABLE]

where $f$ and $g$ are smooth function of the phases which can also depend on time.

These identically coupled oscillators evolve under the action of the Moebius symmetry group $M$ (actually a Moebius subgroup). So, ensembles of identical, globally-coupled oscillators have $N-3$ constants of motion and their dynamics is three-dimensional Watanabe and Strogatz (1993); Marvel et al. (2009). That is, all phases evolve according to the action of the same Moebius transformation

[TABLE]

where the $\theta_{j}$ are constants and $\alpha\in\mathbb{C}$ and $\psi\in\mathbb{S}^{1}$ are the parameters of the Moebius group.

This approach can do more. In the limit of large $N$ it is possible to obtain nonlinear equations for the order parameter. Choosing uniformly-distributed constants of motion $\theta_{j}$, the complex order parameter $z$ follows a Riccati equation

[TABLE]

This reduction was applied to study a number of nonlinear problems in arrays of Josephson junctions Vlasov and Pikovsky (2013); Marvel and Strogatz (2009), discontinuous transitions in explosive synchronization Vlasov et al. (2015) and to classify the attractors in the ensemble of oscillators. Indeed, the only attractors are fixed points or limit cycles where all but one oscillator are synchronized Engelbrecht and Mirollo (2014).

The above reduction follows from the group symmetry of the equations and it is valid only for identical frequencies. Ott and Antonsen (2008) put forward a scheme allowing dimensional reduction for nonidentical frequencies. In the limit $N\rightarrow\infty$, the state of the oscillator system is described by a distribution

[TABLE]

where c.c. stands for complex conjugate. Next we assume that

[TABLE]

that is, the whole distribution $f$ is determined by only one function with $|\alpha|<1$. It is possible to show that the evolution of the system preserves this form of $f$. For various classes of distribution $g$ it is possible to obtain equations for $\alpha$ and thereby for the order parameter.

So the scheme will give low dimensional equations for the order parameter. This ansatz of Eq. (31) has been successfully applied to understand the dynamics of globally coupled oscillators and the second order Kuramoto model Rodrigues et al. (2016), and to understand the formation of clusters when higher order harmonics are included in the coupling function Skardal et al. (2011). The approach can also be used to study nonautonomous globally coupled ensembles of phase oscillators Petkoski and Stefanovska (2012).

Pikovsky and Rosenblum (2011) made a generalization to heterogeneous ensembles of phase oscillators, and connected the Watanabe and Strogatz reduction to the Ott and Antonsen ansatz. In the limit of infinitely many oscillators, the Kuramoto order parameter $z$ can be written as an integral over the stationary distribution of phases $\nu$. Clearly $z$ does not characterize the distribution $\nu$, so one may consider the generalized order parameters

[TABLE]

which are the Fourier coefficients of the distribution $\rho$. Clearly, $z_{1}=z$ is the standard order parameter. If the distribution of the constants of motion is uniform, then

[TABLE]

and, for this particular case, the order parameter $z$ completely determines the distribution. The Ott and Antonsen ansatz corresponds to the special case where the generalized order parameters are expressed via the powers of order parameter.

Sensitivity to the coupling function. This approach to finding low-dimensional dynamics is dependent on the coupling function being sinusoidal in shape. The Watanabe and Strogatz reduction for $N$ globally coupled oscillators gives $N-3$ constants of motion. Therefore, the dynamics of the ensemble is neutral. The dynamics on these subspaces evolves under the identity map. Recent results show that, by perturbing the identity map, we can generate any dynamics Turaev (2015). So, small perturbations in the coupling function can lead to abrupt changes in the dynamics of the ensemble.

III.3.5 Noise and nonautonomous effects

If the oscillators are subject to noise, the phase reduction scheme can still be applied but with some minor modifications Balanov et al. (2008); Ermentrout and Saunders (2006). Even in the absence of coupling, the oscillators can synchronize if driven by a common noise. This is a general result by Le Jan (1987) who showed that, when phase oscillators are driven by noise, the trajectories converge to a random fixed point (corresponding to the two oscillators going to the same trajectory). This result is well appreciated in the physics community as are also the differing effects of common and independent noises Lindner et al. (2004).

In the context of interacting oscillators we can analyze the contributions of the coupling function and common noise in driving the oscillators towards synchronization Garcia-Alvarez et al. (2009). Consider the following model of two coupled (or uncoupled) phase oscillators with common and independent noises:

[TABLE]

where $\xi,\xi_{1}$ and $\xi_{2}$ are Gaussian noises of unit variance and, again, $\psi=\phi_{1}-\phi_{2}$. It is possible to obtain a stochastic differential equation for the generalised phase difference $\psi$. This equation is nonautonomous. An analytical approach is to write a Fokker-Planck equation for the probability density of $\psi$.

The probability density is almost independent of the fast variables $\phi_{1}$ and $\phi_{2}$, so a good approximation is to integrate over these variables to obtain a proxy for a stationary probability distribution. This approach reveals three important effects: $(i)$ independent noises $\xi_{1},\xi_{2}$ hinder synchronisation; $(ii)$ coupling-induced synchronisation takes place for low noise intensity and large coupling strengths; and $(iii)$ common-noise-induced synchronisation occurs for large common-noise intensities and small coupling strengths.

Nonautonomous effects. If the frequencies of the oscillators are nonautonomous but the oscillators are identical the reduction techniques can still be applied Watanabe and Strogatz (1993); Marvel et al. (2009); Petkoski and Stefanovska (2012), so that phenomena such as synchronization can be studied. A new class of systems described by nonautonomous differential equations are chronotaxic systems. These are defined as dissipative dynamical systems with internal sources of energy. In such cases the coupling function is nonautonomous and the systems retain stable (time-dependent) amplitude and phase under external perturbation Suprunenko et al. (2013).

III.4 Networks of oscillators

In this section, we generalize the discussion to networks of interacting systems with pairwise interaction. That is, we consider

[TABLE]

where $W_{ij}$ is the matrix encoding the strength of interaction between $j$ and $i$, and $H$ is the coupling function. Note that we allow each isolated vector field to be distinct. To be able to draw conclusions about the overall dynamics from the microscopic data for $f$, $W$ and $H$, we will consider a subclass of vector fields and coupling functions.

III.4.1 Reduction to phase oscillators

Assume that for $\varepsilon=0$, each isolated system has an exponentially-attracting periodic orbit. A typical assumption is that $H_{ij}=H$ (i.e. all coupling functions are identical). Then proceeding in the same way as in Sec. III.3 for globally coupled oscillators in the limit of small coupling strengths, we can reduce the dynamics to the phases

[TABLE]

This model describes the dynamics of the phase oscillators in terms of complex networks of interactions. Most results relate to the sinusoidal coupling function $q(\varphi_{i},\varphi_{j})=\sin(\varphi_{j}-\varphi_{i})$. The main questions lie in the realm of collective dynamics and transitions from incoherent to coherent states Rodrigues et al. (2016). The situation here is less-well-understood. For instance, it is unclear how to generalize the low-dimensional reduction approach.

III.4.2 Networks of chaotic oscillators

A subclass of this model offers insight. Consider $W_{ij}=A_{ij}$, where $A_{ij}=1$ if $i$ receives a connection from $j$, and $A_{ij}=0$ otherwise. Moreover, consider the diffusive coupling functions $H_{ij}(x,y)=H(x-y)$. Suppose also that $f_{i}=f$, that is, all isolated nodes are identical. We also assume that the isolated systems are chaotic. This assumption is not necessary but æsthetically pleasant, because in this situation the only possible source for collective dynamics is through the coupling. This model corresponds to identical oscillators interacting diffusively, and it can then be rewritten as

[TABLE]

The role of the coupling function is to attempt to bring the system towards synchronization $x_{1}=x_{2}=\cdots=x_{N}$. The main questions now are: (a) when will the coupling function $H$ bring the system towards synchronization; and (b) how will the interaction structure $A_{ij}$ influence the system? We should analyse the growth of small perturbations $x_{i}=s+\xi_{i}$, where $\dot{s}=f(s)$. But we face the challenges of having too many equations and, moreover, of all the $\xi_{i}$ being coupled.

Global Results: Here we want to find the conditions on the coupling function guaranteeing that the network dynamics will converge to synchronization, regardless the initial conditions. The challenge is to construct a Lyapunov function, whose existence is a sufficient condition for a globally stable synchronous state.

Pogromsky and Nijmeijer (2001) used control techniques and concepts of passive systems to obtain global synchronization results for arrays of interacting systems. Assuming that the coupling function is positive definite, they were able to construct a Lyapunov function for the array and to express its construction in terms of the spectrum of the network. They showed thereby that all solutions of the coupled equation are bounded and that, if the coupling is large enough, the network synchronizes. They also showed how diffusion-driven instabilities can appear in such arrays. This approach to passive systems was subsequently applied to neuron models to study their synchronization properties Steur et al. (2009).

In a similar spirit to constructing Lyapunov functions, Belykh et al. (2004) developed the connection graph stability method. At its very heart, the method requires the existence of a Lyapunov function for the nonlinear equations of the perturbations $\xi_{i}$. The existence of this function is unclear from the beginning, however, and it may depend on the vector field $f$. The method relates the critical coupling necessary to attain synchronization to the total length of all paths passing through an edge on the network connection graph.

Another approach to studying synchronization is to tackle the equations for the perturbations $\xi_{i}$ using the theories of contraction Russo and Di Bernardo (2009) and exponential dichotomy Pereira et al. (2014). Here we use the coupling function to construct differential inequalities. At their cores, these approaches are equivalent to the construction of local Lyapunov functions. However, if the coupling function is such that the contraction theory can be applied (for example the coupling function is positive definite) then much information on the synchronization can be extracted. In particular, even if the network structure is time-varying the network may synchronize Lü and Chen (2005).

Local Results: In the section above, we took account of the nonlinear behaviour of the perturbation $\xi_{i}$. In the local approach we consider only the linear terms in $\xi$. Pecora and Carroll (1998) had the idea of block-diagonalizing the perturbations $\xi_{i}$ via a change of coordinates where $\xi_{i}$ goes to $\zeta_{i}$. In the new variables $\zeta_{i}$ the perturbations decouple and they all have the same form

[TABLE]

where $\Gamma=DH(0)$. To recover the equation for $\zeta_{i}$ we only need to set $\alpha=\varepsilon\lambda_{i}$ where $\lambda_{i}$ is the $i$th eigenvalue of the Laplacian of the network. So the problem reduces to the case of two coupled oscillators. Obviously, there are additional challenges in understanding the graph structure via $\lambda_{i}$, but the main idea now boils down to the case of two oscillators. We classify the stability of the variational equation for $z$. The commonest criterion used for stability is the Lyapunov exponent, which gives rise to the master stability function $\Lambda$ (just as in the two oscillator case). This approach showed that the topology of the networks can exert systematic influences on the synchronization Barahona and Pecora (2002) and can be used to predict the onset of synchronization clusters Pecora et al. (2014); Williams et al. (2013). In the last two decades this approach has been popular and it has been applied to a variety of network structures Arenas et al. (2008) and to problems of pinning control Sorrentino et al. (2007).

The MSF approach was also extended to the case where the coupling function has time delays Li and Chen (2004). Moreover, in the limit of large delays it is possible to understand the behaviour of the level sets of the master stability function $\Lambda$. Indeed, the level sets tend to be circles whose radii increase monotonically in the complex plane Flunkert et al. (2010). Some networks also have two types of coupling function. Typical examples are neural networks where electrical and chemical synapses coexist. If the underlying matrices defining the chemical and electrical coupling commute, then the MSF can be used to understand the net effects of the coupling function on the synchronization Baptista et al. (2010).

Generalizations: So far, we have discussed networks of identical oscillators. If the network is composed of slightly nonidentical nodes, the MSF approach can still be applied, via a perturbation analysis Sun et al. (2009). In general, to understand the effect of the network, combined with the effects of nonidentical nodes, further information about the coupling function is necessary. If the coupling function $\Gamma$ (linearized about the synchronized manifold) has a spectrum with a positive real part, then we can extract a great deal of information it. Adding random links to a network of nonidentical oscillators can substantially improve the coherence Pereira et al. (2013). For directed networks, depending on the coupling function, improvements in the network topology such as link addition can destabilise synchronization Pade and Pereira (2015). This phenomenon can also be observed in experiments with lasers when the coupling function has a time-delay Hart et al. (2015).

Moreover, if one adds a small perturbation on the nodes of systems interacting in a fully connected network, the collective dynamics will lead to smaller fluctuations than those expected if the oscillators were uncoupled and one applies the central limit theorem Masuda et al. (2010). That is, the central limit theorem would predict fluctuations of order $O(N^{-1/2})$, whereas the collective dynamics gives $O(N^{-1})$. One can classify the regular networks such that the improvement is significant Pereira et al. (2014). For this class of coupling function, one can also gain insight into the speed of synchronization, that is, how fast the system converges towards synchronization. The speed is well known to depend on the network structure Timme et al. (2004); Grabow et al. (2010). For the above class of coupling functions, it is possible to show how the properties of the coupling function and network structure combine to determine the speed of convergence. For example, if the spectrum of the coupling function is real, then the speed towards synchronization is dictated by the real part of the Laplacian spectral gap.

So far efforts have been concentrated on the study of nonidentical nodes, while keeping the coupling function identical. Because the coupling function couples with the network structure of the equations, perturbations in the coupling function can have a drastic impact on the collective dynamics. For instance, if the network has a heterogeneous degree distribution, no perturbations in the coupling function are tolerated. Any perturbation in the large network limit will destabilize the synchronous motion Maia et al. (2015).

IV Methods

IV.1 Inferring coupling functions

Before discussing methods for the inference of coupling functions, we mention earlier discussions and techniques Schiff et al. (1996); Čenys et al. (1991); Schreiber (2000); Stefanovska and Bračič (1999); Arnhold et al. (1999) that paved the way for the subsequent introduction of coupling functions. Although the mathematical and computational facilities of the time did not allow for the full inference of coupling functions, this goal was nonetheless closely approached by different measures that detected the existence of a coupling relationship and characterized its nature.

In a study of this kind, Stefanovska and Bračič (1999) investigated the coupled oscillators of the cardiovascular system from human blood flow signals. Among other methods, wavelet time-frequency analysis was used to detect the instantaneous frequency of the heart, through which the coupling from respiration was assessed. Similarly, many information-based measures were used for assessment of the coupling strength and directionality Schreiber (2000); Paluš and Stefanovska (2003); Barrett and Barnett (2013). Even though these methods are very useful in detecting the net coupling effects, they are essentially directed functional connectivity measures and they are not designed to infer mechanisms.

Another traditional approach for studying interactions is via transfer function analysis Boreman (2001); Sanathanan and Koerner (1963); Saul et al. (1991); Cooke et al. (1998). The transfer function is a mathematical representation that describes the linear relationships between the inputs and outputs of a system considered as a black-box model. Although this approach suffers from some limitations, it has nonetheless been used in the past for understanding interactions, and still is.

The inference of coupling functions involves the analysis of data to reconstruct a dynamical model describing the interactions. 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). The latter has been used historically as means of advanced “filtering”, when one selects and detects the features of interest described by the model – a celebrated example being the Kalman filter Kalman (1960). Fig. 14 presents the main steps in obtaining the coupling function. In short, starting with the data $\mathcal{M}$ from two (or more) interacting dynamical systems, first the appropriate observable data $\tilde{\mathcal{M}}$, like the phase or amplitude, are estimated from the initial data $\mathcal{M}$, so that they can be used by a method that infers a dynamical model from which one can extract the coupling functions.

The data $\mathcal{M}$, are usually represented by amplitude state signals measured dynamically i.e. they describe the time evolution of the system. Very often the raw measurements require pre-processing and pre-estimation procedures. If the systems are of an oscillatory nature, the phase of the periodic signal is extracted; similarly the amplitude can be extracted from the signals. There can be further pre-processing, including filtering within desired intervals, removal of artifacts, noise suppression, removal of common source disturbances, filtering of power supply frequency, etc. The pre-processed signals then act as input for the inference methods.

The inference process aims to reconstruct a model to describe the interacting dynamical systems. It is given by a set of ordinary differential equations (ODEs) or, if there is dynamical noise, by stochastic differential equations (SDEs). The model in Fig. 14 is given in terms of general variable $\chi$, while usually the model uses either the phase or amplitude domain. The dynamics of the system is modeled with a set of base functions, which are usually linearly independent. For example, the set of base functions can be a Fourier series of sine and cosine functions. Base functions can be either linear or nonlinear, and are specified by a set of parameters $c$ that usually act as scaling parameters. When appropriately parameterized by $c$, the base functions then combine to give the all-important coupling function.

The base functions are a part of the model that is assumed to be known beforehand, so the main task of the inference method is to determine the parameters ${\bf c}$ from the data $\tilde{\mathcal{M}}$, given the (SDE or ODE) model $\dot{\chi}$. The choice of the right model can be rather difficult, especially if the dynamical system does not posses some general characteristics. Nevertheless, a number of methods exist for optimal model selection Ljung (1998); Burnham and Anderson (2002); Berger and Pericchi (1996). Given a model and a set of data, one can use different methods to perform the inference. These methods may differ considerably in their characteristics and performance, and we present below some examples of those that are most widely used.

IV.2 Methods for coupling function reconstruction

IV.2.1 Modeling by least-squares fitting

As mentioned above, one of the first works on the reconstruction of coupling functions from data was that developed by Rosenblum and Pikovsky (2001). Their inference of the interaction is based on a least-squares fitting procedure applied to the phase dynamics of the interacting oscillators. The main goal of the method is the detection of coupling and directionality. Nevertheless, part of the results are functions that closely resemble the form of coupling functions.

The technique provides for experimental detection of the directionality of weak coupling between two self-sustained oscillators, from bivariate data. The approach makes use of the well-known fact that weak coupling predominantly affects the phases of the oscillators, not their amplitudes. The principal idea is to investigate and quantify whether the phase dynamics of one oscillator is influenced by the phase of the other. To achieve this, the model of the phase equations (Eq. 5) is fitted to the phase data. From the inferred model and its parameters, one can then quantify the coupling in one or the other direction.

First, for each point in time of the phase time-series, the increments $\Delta_{1,2}(k)=\phi_{1,2}(t_{k}+\tau)-\phi_{1,2}(t_{k})$ are computed, where $\tau$ is a free parameter. These increments $\Delta_{1,2}(k)$ are considered as being generated by some unknown two-dimensional map

[TABLE]

The functions $\mathcal{F}_{1,2}[\phi_{1,2}(k),\phi_{2,1}(k)]$ are decomposed into Fourier series, and their dependences $\Delta_{1,2}(k)$ on $\phi_{1}$ and $\phi_{2}$ are modeled with the least-square fitting procedure. As base functions for the fitting, the Fourier series:

[TABLE]

with $|l|\leq 3$ for $m=0$, $|m|\leq 3$ for $l=0$ and $|l|=|m|=1$ were considered.

It is worth pointing that this notion is close to, though not exactly identical to, the dynamical inference of ODEs: the increments $\Delta_{1,2}(k)$ are close to the Euler method for first order differentiation which would have been $\Delta_{Euler;1,2}(k)=[\phi_{1,2}(t_{k+1})-\phi_{1,2}(t_{k})]/h$, where $h$ is the sampling (differentiation) step. Therefore, the functions $\mathcal{F}_{1,2}[\phi_{1,2}(k),\phi_{2,1}(k)]$ are similar to the coupling functions $q_{1,2}(\phi_{1},\phi_{2})$, i.e. they are close to the form of the genuine coupling functions with close relative but not absolute coupling strength. Despite the differences, these were probably the first extracted functions of oscillatory interactions, and they were of great importance for the work that followed.

The inference itself was performed by least-squares fitting, a widely-used method for finding the best-fitting curve to a given set of points Leon (1980); Lawson and Hanson (1974). The main objective of the fitting consists of adjusting the parameters of a model function to best fit a data set. The result of the fitting process is an estimate of parameters, given the model and the base functions. To obtain the parameter estimates, the least-squares method minimizes the summed square of residuals (often also called offsets). The residual $\wp(k)$ for the $k$-th data point is defined as the difference between the observed response value $\Delta(k)$ and the fitted response value $\tilde{\Delta}(k)$, and is identified as the error associated with the data. The summed square of residuals is then given as:

[TABLE]

The main estimation, aiming to minimize $\mathcal{O}$, involves partial differentiation with respect to each parameter, and setting the result equal to zero Leon (1980); Lawson and Hanson (1974). Such schemes can use linear, weighted or nonlinear fitting. More recent methods for coupling function reconstruction with fitting procedures often involve a kernel density estimation Kralemann et al. (2008). The finally inferred parameters applied to the model base functions, provide explicit definitions of the functions $\mathcal{F}_{1,2}[\phi_{1,2}(k),\phi_{2,1}(k)]$.

To demonstrate the method, a simple example of two coupled phase oscillators subject to white noise was considered:

[TABLE]

where the coupling functions are sines of the phase difference (c.f. Kuramoto Eq. (8)), with frequencies $\omega_{1,2}=1\pm 0.1$, and couplings $\varepsilon_{1}=0.1$ and $\varepsilon_{2}=0.02$ (i.e. weaker than $\varepsilon_{1}$). The noise is assumed to be white Gaussian with $\langle\xi_{1,2}(t)\xi_{1,2}(t^{\prime})\rangle=\delta(t-t^{\prime})2D_{1,2}$, with $D_{1,2}=0.2$. One should note that the least-squares fitting only infers the deterministic part of the ODEs, and the noise here is used to introduce imprecisions only i.e. there is no inference of SDEs. Fig. 15 presents the two functions as reconstructed using least-squares fitting. Note from Fig. 15(a) that the form of the reconstructed function $\mathcal{F}_{1}$ rightly resembles that of the genuine coupling function $\sin(\phi_{2}-\phi_{1})$, i.e. a diagonal form of a wave determined by the phase difference $\phi_{2}-\phi_{1}$; and (b) the strength, or amplitude of $\mathcal{F}_{2}$ is much lower due to the weaker coupling strength.

IV.2.2 Dynamical Bayesian inference

The recently introduced method for the dynamical Bayesian inference of coupling functions Stankovski et al. (2012) relies on a Bayesian framework von Toussaint (2011); Smelyanskiy et al. (2005); Friston (2002); Bayes (1763) and is applied to a stochastic differential model where the deterministic part is allowed to be time-varying.

The method attempts to reconstruct the coupling functions by inferring a model consisting of two weakly-interacting dynamical systems subject to noise. The model to be inferred is described by the stochastic differential equation

[TABLE]

where $i\neq j={1,2}$, and $f({\bf\chi}_{i},{\bf\chi}_{j}|{\bf c})$ are base functions describing the deterministic part of the internal and the interacting dynamics. The parameter vector ${\bf c}$ provides scaling coefficients for the base functions. The noise is assumed to be white, Gaussian, and parameterized by a noise diffusion matrix ${\bf D}$. At this point we speak of ${\bf\chi}_{i}$ in general, but later we will refer separately to the phase or amplitude coupling functions depending on the domain of the application.

Given the $2\times M$ time-series ${\mathcal{X}}=\{{\bf{\bf\chi}}_{n}\equiv{\bf\chi}(t_{n})\}$ ($t_{n}=nh$) provided, and assuming that the model base functions are known, the main task for dynamical inference is to infer the unknown model parameters and the noise diffusion matrix ${\mathcal{P}}=\{{\bf c},{\bf D}\}$. The problem eventually reduces to maximization of the conditional probability of observing the parameters ${\mathcal{P}}$, given the data ${\mathcal{X}}$. For this Bayes’ theorem Bayes (1763) is applied, exploiting the prior density $p_{\mbox{\scriptsize prior}}(\mathcal{P})$ of the parameters and the likelihood function $\ell(\mathcal{X}|\mathcal{P})$ of observing $\mathcal{X}$ given the choice $\mathcal{P}$, in order to determine the posterior density $p_{{\mathcal{X}}}({\mathcal{P}}|{\mathcal{X}})$ of the unknown parameters ${\mathcal{P}}$ conditioned on the observations ${\mathcal{X}}$:

[TABLE]

The next task is to determine the likelihood functions in order to infer the final posterior result. From the time-series the midpoint approximation ${{\bf\chi}}_{n}^{\ast}=({\bf\chi}_{n}+{\bf\chi}_{n+1})/2$ is constructed, followed by the Euler differentiation $\dot{}{\bf\chi}_{n}=({\bf\chi}_{n+1}-{\bf\chi}_{n})/h$. Use of the stochastic integral for noise that is white and independent leads to the likelihood function, which is given by a product over $n$ of the probabilities of observing ${\bf\chi}_{n+1}$ at each time Smelyanskiy et al. (2005). The negative log-likelihood function is then $\mathcal{S}=-\ln\ell({\mathcal{X}}|{\mathcal{P}})$ given as:

[TABLE]

with implicit summation over the repeated index $k$. The likelihood (33) is of quadratic form. Thus if the prior is a multivariate normal distribution, so also will be the posterior. Given such a distribution as a prior for the parameters ${\bf c}$, with mean $\bar{{\bf c}}$, and covariance matrix ${{\bf\Sigma_{\mbox{\scriptsize prior}}\equiv\Xi}^{-1}}_{\mbox{\scriptsize prior}}$, the final stationary point of $\mathcal{S}$ is calculated recursively from:

[TABLE]

where summation over $n=1,\ldots,N$ is assumed and the summation over repeated indices $k$, and $w$ is again implicit. The initial prior can be set to be the non-informative flat normal distribution, ${{\bf\Xi}}_{\text{prior}}=0$ and $\bar{}{\bf c}_{\mbox{\scriptsize prior}}=0$. These four equations (34) are the only ones needed for implementing the method. They are applied to a single block of data ${\mathcal{X}}$ and the resultant posterior multivariate probability ${\mathcal{N}}_{\mathcal{X}}(c|\bar{c},\Xi)$ explicitly defines the probability density of each parameter set of the model (32).

In dynamical Bayesian inference each new prior distribution depends on and uses the previously inferred posterior distribution. In this framework, however, the information propagation is amended in order to allow the method to follow the time-variability of the parameters Stankovski et al. (2012). The new prior covariance matrix becomes $\Sigma_{\text{prior}}^{n+1}=\Sigma_{\text{post}}^{n}+\Sigma_{\text{diff}}^{n}$, where $\Sigma_{\text{diff}}^{n}$ describes how much some part of the dynamics can change with time.

Given the use of Bayesian inference with informative priors, the method is not prone to the overfitting of parameters, and it does not require much data within the windows because, in each new block of data, it only updates the parameters Duggento et al. (2012). For analyses of dynamical oscillators, one can use data windows containing 6 to 10 cycles of the slowest oscillation; care is needed to ensure that the windows are long enough in cases where there is modulation that is slow relative to the eigenfrequencies Clemson et al. (2016). The confidence of the fit is given by the resultant covariance matrix $\Sigma_{\text{post}}$.

The description above is for two interacting oscillators. Nonetheless, the theory also holds for a larger number of oscillators and the dynamical Bayesian inference has been generalized to infer networks of systems with multivariate coupling functions Stankovski et al. (2015).

Fig. 16 shows an application of dynamical Bayesian inference to cardiorespiratory interactions from a resting human subject whose paced respiration was ramped down with decreasing frequency. The inference of the dynamics and the coupling functions were reconstructed from the cardiorespiratory phase dynamics Stankovski et al. (2012). Fig. 16(a) indicates the reconstructed respiration frequency, showing the linearly-decreasing trend. The inferred coupling directionality, defined as $Dirc=(\varepsilon_{r}-\varepsilon_{h})/(\varepsilon_{r}+\varepsilon_{h})$, is also time-varying, with predominant direction of influence from the respiration to the heart. The reconstructed cardiorespiratory coupling functions Fig. 16(c)-(e), are described by complex functions whose form changes qualitatively over time – cf. Fig. 16(c) with Fig. 16(d),(e). This implies that, in contrast to many systems with time-invariant coupling functions, the functional relations for the interactions of an open (biological) system can themselves be time-varying processes. By analyzing consecutive time windows, one can follow the time evolution of the coupling functions.

The time-variability of biological systems and the ability of the method to reconstruct it has implications for the detection of chronotaxic systems, which are a class of nonautonomous self-sustained oscillators able to generate time-varying complex dynamics Suprunenko et al. (2013). Such systems have drive-response subsystems which are inherently connected with appropriate coupling functions. Dynamical Bayesian inference has been applied to reconstruct such chronotaxic systems for cases where the model was known or could be closely approximated Clemson et al. (2014).

IV.2.3 Maximum likelihood estimation: multiple-shooting

The reconstruction of coupling functions has been performed using techniques for maximum likelihood estimation, an approach that was employed for reconstruction of the coupling functions of electrochemical interactions Tokuda et al. (2007, 2013). For this a multiple-shooting method, as one type of maximum likelihood estimation, was used.

The maximum likelihood estimation Aldrich (1997); Myung (2003) is a statistical method of seeking that probability distribution which makes the observed data most likely, which means that one needs to find the value of the parameter vector that maximizes the likelihood function. The procedure of maximization, intuitively describes the “agreement” of the selected model with the observed data, and for discrete random variables it maximizes the probability of the observed data under the resulting distribution. Maximum likelihood estimation gives a unified approach to estimation, which is well-defined in the case of the normal distribution and many other problems. It is of fundamental importance in the theory of inference and provides the basis for many inferential techniques in statistics.

Maximum likelihood estimation is in general different from least-squares fitting (Sec. IV.2.1), as the former seeks the most likely parameters, while the latter is a descriptive tool that seeks the parameters that provide the most accurate description of the data. There is a situation, however, in which the two methods intersect and the same parameters are inferred. This is when observations are independent of one another and are normally distributed with a constant variance Myung (2003).

The multiple-shooting method used for coupling function reconstruction is based on maximum likelihood estimation Baake et al. (1992); Tokuda et al. (2013); Voss et al. (2004). A known approach for inference of the trajectories and the parameters is the so-called initial value approach, where initial guesses for the states $x(t_{1})$ and parameters ${\bf c}$ are chosen and the dynamical equations are solved numerically. However, a problem can appear in such approaches – the inferred trajectory may converge only to a local maximum.

The multiple shooting algorithm provides a possible solution to the problem. In this approach, initial conditions are estimated at several points along the time series, so that the shooting nodes, and thus the estimated trajectory, can be made to stay closer to the true values for a longer time. This task is considered as a multi-point boundary value problem. The interval for fitting $(t_{1},t_{N})$ is partitioned into $m$ subintervals:

[TABLE]

Local initial values $x_{j}=x(\tau_{j})$ are introduced as additional parameters for each subinterval $(\tau_{j},\tau_{j+1})$. In the case of independent Gaussian noise, maximization of the likelihood amounts to minimization of the cost function $\zeta^{2}(x_{1},x_{2}\ldots,x_{m},{\bf c})$, which is the sum of the squared residuals between the data and the model trajectory, weighted by the inverse variance of the noise:

[TABLE]

where $x$ and $y$ are the state and the observed data, respectively, $G$ is a function for the dynamics, and $\sigma_{i}$ gives the noise variance. Thus in the multiple-shooting method the dynamical equations are integrated piecewise and the cost function is evaluated and minimized on the multiple samples from each subinterval.

Assuming that the dynamical parameters ${\bf c}$ are constant over the entire interval, the local initial values are optimized separately in each subinterval. The latter leads to an initially discontinuous trajectory and the final step is to linearize them so as to provide continuous states. This task, called condensation, is often achieved by use of the generalized Gauss-Newton method.

The multiple-shooting method has been used to model the phase dynamics of interacting electrochemical oscillators in order to reconstruct their coupling functions Tokuda et al. (2007), using Fourier series as base functions. The particular application used an electrochemical oscillatory system in which the coupling function had previously been calculated Kiss et al. (2005) (see also Sec. II.5) by applying the perturbation method to a single oscillator, and thus a direct comparison could be made between the two approaches.

Fig. 17 shows the inference of a network of 32 electrochemical oscillators. The effective natural frequencies are well estimated, with slightly higher values than those obtained with for completely uncoupled systems Fig. 17(a). The form of the estimated coupling function is in reasonable agreement with that found by applying the perturbation method to a single isolated electrochemical oscillator Fig. 17(b), with a difference in amplitude of 23.7% between the two. The coupling function is said to be of a form consistent with theoretical predictions for Stuart-Landau oscillators close to a Hopf bifurcation. In a similar way, the technical dependences and conditions, including dependence on the observational noise, network size, number of defects, and data length have also been examined Tokuda et al. (2013).

IV.2.4 Random phase resetting method

The method of random phase resetting can be used for the dynamical inference of interacting systems, and also for the reconstruction of their coupling functions Levnajić and Pikovsky (2011). Initially, the method was designed for the reconstruction of network topology i.e. the inference of coupling strengths; nevertheless, the framework employed dynamical inference and the inferred model allows for the coupling functions to be reconstructed as well.

The main idea relies on repeatedly reinitializing the network dynamics (e.g., by performing random phase resets), in order to produce an ensemble of the initial dynamical data. The quantities obtained by averaging this ensemble reveal the desired details of the network structure and the coupling functions.

The method considers a dynamical model of interacting phase oscillators such as that in Eq. 5, with coupling functions in terms of the phase difference $q_{i}(\phi_{j}-\phi_{i})$. The dynamics starts from a set of initial phases which are denoted as $\bm{\phi}=(\phi_{1},\ldots,\phi_{N})(t=0)$, chosen from a distribution $\hbar(\bm{\phi})>0$ normalized to $(2\pi)^{N}$. The method is based on two assumptions: (i) that one is able to arbitrarily reinitialize the network dynamics $\mathcal{K}$ times, by independently resetting the phases of all nodes to a new state $\bm{\phi}$, and (ii) that one is able to measure all the values $\bm{\phi}_{l}$, and all initial instantaneous frequencies $\dot{\bm{\phi}}_{l}$, each time the dynamics is reinitialized (for $l=1,\ldots,\mathcal{K}$). A $2\pi$-periodic test function $\mathcal{T}=\mathcal{T}(\phi_{i}-\phi_{j})$ with zero mean is given as input and the coupling functions are taken to be represented by Fourier series, to obtain an expression for the index $\Upsilon_{ij}$:

[TABLE]

The dynamical network described by the phase equations can be reconstructed by computing $\Upsilon_{ij}$ for a suitable $\mathcal{T}$ function, e.g. $\mathcal{T}(\phi)=e^{in\phi}$. The practical implementation in terms of the data involves the representation of $\dot{\bm{\phi}}_{l}$ with a kernel smoother Wand and Jones (1994); Kralemann et al. (2008), and appropriate averaging, to get $\Upsilon_{ij}[\mathcal{T}]$:

[TABLE]

To mimic an experimentally feasible situation, $\mathcal{K}$ random phase resets of the network dynamics separated by the time interval $\tau$ are performed. A network of four phase oscillators, coupled as shown in Fig. 18(a), is considered. After applying the inference method, the reconstructed pairwise connections are shown to be in reasonable agreement with the actual coupled values, and especially in identifying the non-couplings, as shown in Fig. 18(b). Once the parameters of the phase model $a^{(1)}_{ij}$ and $b^{(1)}_{ij}$ have been inferred, one can also reconstruct the form of the coupling function. Fig. 18(c) presents an example of the coupling function showing the influence that the fourth oscillator is exerting on the second oscillator.

The approach is related to the methods for reconstruction of phase response curves (see Sec. IV.4.2). Here, however, phase resetting is used somewhat differently i.e. the focus is on the network’s internal interactions, rather than on its response to stimuli. The power of this method lies in a framework that yields both the topology and the coupling functions. Its downside is that it is invasive – requiring one to interfere with the on-going system dynamics (via phase-resets or otherwise), which is often experimentally difficult and sometimes not even feasible.

IV.2.5 Stochastic modeling of effective coupling functions

An important feature of the interacting dynamics in real systems is the presence of noise. Explicit consideration of the stochastic nature of the dynamics can provide a better means of describing the coupling functions, and how they are affected. In Sec. IV.2.2 we discussed the dynamical Bayesian method which treats stochastic dynamics, and is able to infer the deterministic part of the coupling function separately from the random noise. Often when considering noise-induced oscillations, however, one may wish to determine the effective coupling functions including the effect of noise Schwabedal and Pikovsky (2010). When performed on the effective phase dynamics with an invariant phase defined in a nonperturbative way, the phase will depend on the noise intensity, and so will all the corresponding characteristics such as the coupling function.

The authors consider an effective phase model describing periodically-driven, noise-induced oscillations:

[TABLE]

where $\vartheta=\Omega t$ is a $2\pi$-periodic driving phase. The aim is to describe the effective phase dynamics $\mathcal{H}(\theta,\vartheta)$ and the corresponding effective coupling functions. One can express the effective dynamics (35) as: $\dot{\theta}=\mathcal{H}(\theta,\vartheta)=\mathcal{H}_{m}(\theta)+\mathcal{F}(\theta,\vartheta)$, where $\mathcal{H}_{m}(\theta)$ is $\vartheta$-independent marginal effective velocity and $\mathcal{F}(\theta,\vartheta)$ is the effective coupling function. By integrating Eq. (35) over $\vartheta$, the marginal effective velocity $\mathcal{H}_{m}(\theta)$ can be determined. Hence, using $\mathcal{F}(\theta,\vartheta)=\mathcal{H}(\theta,\vartheta)-\mathcal{H}_{m}(\theta)$, one can determine the effective coupling function:

[TABLE]

where $P=P(\theta,\vartheta)$ and $P_{m}=P_{M}(\theta)$ are the probability densities of the full and the marginal dynamics, respectively. The variable $\theta$ can be considered as a protophase and can be further transformed by $\phi=\mathcal{C}(\theta)=2\pi\int_{0}^{\theta}P(\eta)d\eta$ to yield an invariant effective phase dynamics:

[TABLE]

Equation (36) provides the effective phase dynamics of the periodically-driven noise-induced oscillations with an effective coupling function $q$ that depends on the noise intensity.

This theoretical description can be illustrated on a noise-driven FitzHugh-Nagumo model as an example of an excitable system:

[TABLE]

where $a=1.1$, $\epsilon=0.05$ are the parameters of the system, $\xi(t)$ is an additive noise which induces oscillations, and the cosine external function provides the interactions in the system. After estimation of the protophase time-series $\theta$ and its transformation to the phase $\phi$, the effective coupling function $q(\Omega t,\phi)$ can be determined. For this a double Fourier series decomposition was used with least-squares fitting of the model to the data. In this sense, the core of the inference is the same as the least-squares fit discussed in Sec. IV.2.1, even though the difference here is that one reconstructs a stochastic model.

The results of the analysis indicated an increase in the effective coupling for vanishing noise, and masking of the coupling for driven noise-induced oscillations of the FitzHugh-Nagumo model. Fig. 19 presents an unusual case with implications for the interpretation of effective coupling functions. Namely, the effective coupling function was computed with two noise intensities for the same coupling strength. By comparing the two plots in Fig. 19 one can see that the amplitude of $q$ decreases with increasing noise intensity. The change in amplitude may have been related to a more pronounced masking of the coupling induced by the frequency shift, or due to the generic decrease in effective coupling for stronger noise because of flattening of the marginal probability.

IV.2.6 Comparison and overview of the methods

The methods discussed for reconstruction of coupling functions possess some characteristics that are in common, as well other features that differentiate them. The latter eventually lead to different choices of method for use, depending on the circumstances and conditions for the dynamics and the coupling functions to be inferred.

Table 4 summarizes the difference and performance of the methods discussed in the previous sections IV.2.1-IV.2.5. Inference of stochastic dynamics which treated SDEs and the influence of dynamical noise that can cause noise-induced qualitative changes (e.g. phase slips), can be performed with dynamical Bayesian inference and stochastic modelling of the effective coupling functions. The other methods treat ODEs with possible measurement noise that can affect the statistics and precision of the inference.

Often the model for inference has more parameters than the real system. In such cases overfitting of parameters can occur and some methods can infer random error or noise instead of the underlying dynamical relationship. A model that has been overfitted will generally be a poor representation of the real system, as it can exaggerate minor fluctuations in the data. Bayesian inference uses informative prior probabilities and can avoid the problem of overfitting parameters.

The speed of calculation varies between methods, especially as some methods perform additional steps and therefore take longer. Dynamical Bayesian inference has a recursive loop, evaluated for each time-point within a data window until a certain precision is reached; the multiple shooting method requires additional initial conditions, the shooting nodes, to be estimated at several points along the time series; the random phase resetting method uses a large number of additional random initial resetting points. These additional processing steps relative to the initial handling of the time-series, inevitably require more computing power and thus reduce the calculation speed.

The coupling functions are usually evaluated for a sequence of time-series each defined by a certain window-length whose choice determines how well a method is able to follow the time-evolution of the coupling functions. Dynamical Bayesian inference updates the new probabilities within a window of data, based on prior knowledge; the multiple-shooting method exploits the initial shooting nodes; and the random-resetting method also uses resetting points, which are said to require shorter data windows and, in turn, to provide good time-resolution of the inference.

A difficulty in common for all the methods is the reconstruction of coupling functions (and coupling strength in general) when the systems are highly synchronized and coherent in the 1:1 frequency ratio Rosenblum and Pikovsky (2001); Kiss et al. (2005). Namely, in the 1:1 phase synchronous state there is a definite strong relationship between the phases and the trajectory on a $(\phi_{1},\phi_{2})$ torus which is one line; hence $\phi_{1}$ and $\phi_{2}$ are not independent, and the coupling functions of the two variables cannot be separately estimated i.e. one is not able to separate the effect of interaction from the internal dynamics of autonomous systems. In order to obtain information about the coupling one needs to observe deviations from synchrony, either due to dynamical noise or due to the onset of quasiperiodic dynamics outside the synchronization region. Synchronous states with larger $n$:$m$ frequency ratios are favorable, because many revolutions cover the surface of the torus, and the inference is then more successful.

A similar situation applies for the inference of coupling functions of dynamical systems in amplitude states. In such cases, the systems are multidimensional, e.g. two coupled Lorenz chaotic systems, and the inference of the coupling functions is more plausible in a 1:1 generalized synchronization sense Stankovski et al. (2014). Complete synchronization and very strong coupling will again constrict the available space for inference, leading to difficulties in reconstructing the amplitude coupling functions.

IV.3 Towards coupling function analysis

Often one needs to draw conclusions based on comparison and quantitative measures of the coupling functions, after they have been inferred. Such situations occur in experimental studies of some real system interactions, e.g. in biomedicine or chemistry. For example, the biomedical expert would like to have a quantitative measure of the coupling functions to be able to describe or compare different states or behaviours in health and disease.

One needs to quantify some characteristics that describe the coupling functions, and in particular features that are unique to the coupling functions and cannot be obtained from other measures. The form of the function can describe the mechanism of the interactions, so being able to quantify it is of obvious interest.

Quantifying a function is not a trivial task, in general. However, many coupling functions can be decomposed, or are inferred through decomposition into functional components, like for example when the phase coupling functions are decomposed into Fourier series. Therefore, the problem of quantification of the coupling function can be reduced to quantification of its components, and in particular the parameters obtained for the components. In this way, one is left to work with a one-dimensional vector of parameters.

One way to perform such a quantification is through the correlation coefficient and the difference measure evaluated from the inferred coupling functions Kralemann et al. (2013a). The first index $\rho$ measures the similarity of two coupling functions $q_{1}$ and $q_{2}$, irrespectively of their coupling strengths. It is calculated as the correlation coefficient:

[TABLE]

where $\langle\circ\rangle$ denotes spatial averaging over the 2D domain $0\leq\phi_{1},\phi_{2}\leq 2\pi$, $\tilde{q}=q-\langle q\rangle$ and $||q||=\langle qq\rangle^{1/2}$. The similarity index $\rho$ is of great interest as it relates only to the form of the function and is a unique measure of the coupling functions. In a similar way, the difference measure is defined as:

[TABLE]

The difference measure is of less interest than $\rho$ as it relates to the coupling strengths, which can be assessed in different ways through other measures.

Two measures were used to analyze the similarity and difference of cardiorespiratory coupling functions, as shown in Fig. 20. It was found that the functions have a well-pronounced characteristic shape for each of the subjects and that the correlations between the coupling functions obtained in different trials with the same subjects were very high Fig. 20(a). Naturally, the correlation between the functions of different subjects was lower, reflecting the interpersonal variability; however, it is high enough to demonstrate the high similarity of the interactions in the group of subjects. The similarity of the coupling functions, obtained from different observables such as the ECG and the arterial pulse for the cardiac oscillations, support the validity of the use of invariant phase. The similarity index $\rho$ has been also used for quantifying the form of the brain coupling functions Stankovski et al. (2015), quantifying significant differences in the form of the coupling functions when altered by the use of different anaesthetics Stankovski et al. (2016).

IV.4 Connections to other methodological concepts

IV.4.1 Phase reconstruction procedures

When analysing data one needs first to reconstruct the phase, before attempting to detect the underlying phase coupling functions. Various methods exist for phase reconstruction from data, including the marked events method (the marking a particular time event, e.g. a maximum or a zero-crossing, within a cycle of oscillation), the Hilbert transform, and wavelet transform based methods Quiroga et al. (2002); Daubechies et al. (2011); Iatsenko et al. (2016); Pikovsky et al. (2001); Gabor (1946). The effect of the method used can have a direct impact on the form of the reconstructed coupling function. It is therefore important to choose a method to reconstruct a phase that is as genuine as possible. For example, the marked events method reduces the inter-cycle resolution and, despite its limited usefulness in synchronization analysis, it is not appropriate for coupling function analysis. Kralemann et al. (2008) proposed a protophase-to-phase transform that obtains an invariant phase in terms of the genuine, observable-independent phases. This technique can be very useful in checking consistency with the phase estimated by use of the Hilbert transform. Some other phase estimates have also been discussed, noting that use of the synchrosqueezed wavelet transform means that one does not need explicit protophase-phase preprocessing, as it estimates directly the genuine phase Daubechies et al. (2011); Iatsenko et al. (2013). Recently, Schwabedal and Kantz (2016) introduced a method that facilitates a phase description of collective, irregular-oscillatory dynamics from multichannel recordings and they demonstrated it on EEG recordings. Such phase estimates have a potential for the reconstruction of coupling functions from collective dynamics. In any case, one should be very careful when estimating phases for coupling functions, in particular from experiments, as otherwise this can lead to spurious descriptions of the coupling functions.

IV.4.2 Relation to phase response curve in experiments

The phase response curve (PRC), describes how an oscillator responds to an external perturbation Winfree (1980); Kuramoto (1984); Tass (1999). The response of the affected oscillator is manifested as shift of its phase. It has been used in various fields, especially in biological oscillations including the heartbeat, circadian rhythms and neuronal activity Tass (1999); Oprisan et al. (2004); Ermentrout (1996); Czeisler et al. (1989); Preyer and Butera (2005); Ko and Ermentrout (2009); Hannay et al. (2015).

The phase response curve is a function expressed in terms of one phase variable from the affected oscillator (for detailed theoretical description see Sec. III.2.4). In this way, the phase response curve is a similar concept to a coupling function, with that difference that the latter describes the interactions on the whole (two-dimensional) space i.e. depending on the two phase variables. In fact, the phase response curve is a functional component of the coupling function. In terms of the general theory of phase dynamics Kuramoto (1984); Winfree (1980), the coupling function $q_{1}(\phi_{1},\phi_{2})$ can be expressed as the product of two functions:

[TABLE]

where $Z_{1}(\phi_{1})$ is the phase response curve, while $I_{1}(\phi_{2})$ is the perturbation function.

The reconstruction of functional curves from data has been used widely, elucidating the mechanisms underlying the oscillations found in nature Tass (1999); Czeisler et al. (1989); Batista et al. (2012). This approach is widely used in neuroscience Schultheiss et al. (2011); Tateno and Robinson (2007); Gouwens et al. (2010); Ermentrout (1996); Galán et al. (2005). For example, the phase response curve has been estimated with electrophysiological experiments on real neurons from the mouse olfactory bulb Galán et al. (2005). A constant current was injected into the neuron to make it fire at a constant frequency within the beta/gamma frequency band. By following the responses of the neurons to the injected current stimulation, the phase response curve was reconstructed. Fig. 21(a) shows the experimental dots and the fitted phase response curve, which matches well the one from the phase model of the study. The surrogate estimation (shuffled dots) in Fig. 21(b) validates this result. Thus, the method allowed for a simplification of the complex dynamics from a single neuron to a phase model. This study also demonstrates the relationship to the coupling function, which was reconstructed from the convolution of the phase response curve and the perturbation function – an approach used in chemistry as well Kiss et al. (2005).

Going in the opposite direction, the *phase response curve can be estimated by decomposition of the coupling function *Kralemann et al. (2013a). This can be done by decomposing the reconstructed coupling function into a product of two functions Eq. (37) and searching for a minimum in the decomposition error by means of an iterative scheme. In this way, the interactions coming from the second oscillator are used as the perturbation to the first oscillator under consideration, whence there is no need for additional external stimulation – a procedure referred by the authors as in vivo estimation of the phase response curve. This method was applied to the reconstruction of the cardiac phase response curve as perturbed by the respiratory oscillations, as shown in Fig. 22. One can clearly see the interval where the phase response curve is non-zero, so that the cardiac system is susceptible to the respiratory perturbation. Intervals of positive (negative) effective forcing are the intervals where respiration is accelerating (decelerating) the heart rate.

Interesting and relevant parallels could be drawn between coupling functions and amplitude response curves, or phase-amplitude response curves Huguet and de la Llave (2013); Castejón et al. (2013). The latter are similar to phase response curves, with the difference that there is also a response to the amplitude on increasing or decreasing the strength of the oscillations.

IV.4.3 General effective connectivity modeling

Quite generally, methods of modeling dynamical systems from data often contain coupling functions von Toussaint (2011); Smelyanskiy et al. (2005); Friston (2011); Voss et al. (2004). The extent to which these coupling functions resemble the same concept as that discussed in this review can vary, depending on the design of the method and the model itself. For example, there can be a model of one larger system which is different from the interaction of two or many systems, but there can be functions within the model that are coupling certain variables or dimensions.

Similar implications hold for the general description of methods for effective connectivity modeling which exploit a model of differential equations and allow for dynamical mechanisms of connectivity to be inferred from data. Such effective connectivity has particularly wide use in neuroscience, where the methods infer the links on different scales of connectivity and spatially distributed regions within the heavily connected brain network. Although such methods have much in common with coupling function inference methods they do not, however, consider the coupling function as an entity, and nor do they assess or analyze the coupling functions as such.

V Applications and Experiments

In this section we review a number of important applications of the methods for reconstruction of coupling functions and their use for the study and manipulation of the interactions, in various fields.

V.1 Chemistry

The interactions of chemical oscillations have been studied extensively, including in connection with coupling functions Kiss et al. (2005, 2007, 2002); Blaha et al. (2011); Tokuda et al. (2007, 2013); Miyazaki and Kinoshita (2006); Kori et al. (2014). Experiments on chemical, or electrochemical, oscillations provide a convenient way of studying and manipulating interactions and coupling functions under controlled laboratory conditions.

One of the more prominent coupling function applications to chemical oscillators is for engineering complex dynamical structures Kiss et al. (2007). The work exploits the simplicity and analytical tractability of phase models and, in particular, their reconstructed coupling functions in order to design optimal global, delayed, nonlinear feedback for obtaining and tuning the desired behavior. It uses a feedback design methodology capable of creating a large class of structures describable by phase models for general self-organized rhythmic patterns in weakly interacting systems with small heterogeneities. The electrochemical oscillations were achieved with electrode potentials during the electrodissolution of nickel wires in sulfuric acid.

The engineering of the interactions to the desired behavior is achieved in a population of $N$ oscillators through the imposition of nonlinear, time-delayed feedback in the amplitude state. This reduces to a phase model of a population of oscillators with weak, global (all-to-all) coupling described by the Kuramoto model Kuramoto (1984) with a general diffusive coupling function $q(\phi_{j}-\phi_{i})$ (i.e. notation $H(\phi_{j}-\phi_{i})$ in this work). In this way, one can also derive the phase response function $Z(\phi_{i})$ in connection to the feedback function. Given such a feedback function and response function $Z(\phi_{i})$, one can in principle obtain the coupling function $q(\phi_{j}-\phi_{i})$ for use in the phase model. However, in the work discussed, Kiss et al. (2007) proceeded in the opposite manner: they chose a coupling function to produce the desired states, and then designed a feedback loop with optimized feedback parameters to give the desired $q(\phi_{j}-\phi_{i})$.

The method is demonstrated with three interesting and important experiments: (i) tuning the desired arbitrary phase differences between two dissimilar oscillators (see also Fig. 1); (ii) generation of complex patterns that include self-organized switching between unstable dynamical states and clusters; and (iii) the physiologically important problem of desynchronization of oscillators. Below, we devote particular attention to case (ii) involving the generation of sequential states and clusters Ashwin and Timme (2005).

Quadratic feedback to a population of four oscillators is used to reproduce a coupling function proposed for slow switching, Fig. 23(a). The experimental system with feedback that sequentially visits (unstable) two-cluster states with two oscillators in each cluster shows two (saddle type) cluster states in state space, Fig. 23(b). In agreement with the experiments, the phase model predicted a switch between these states due to the existence of heteroclinic orbits connecting them. These switches can be seen as a fluctuation of the system order, as shown in Fig. 23 (c). The engineered feedback produces configurations of two clusters, each containing two elements, connected by heteroclinic orbits. Two types of transitions have been observed: intra-cluster and inter-cluster transitions as presented by the trajectories of the experimental system, and illustrated as phase space plots in Fig. 23(d) and (e).

In a similar way, Kiss et al. (2005) developed a method for reconstruction of coupling functions from electrochemical oscillations, which are then used to predict synchronization, as also discussed in Sec. II.5. Similarly, Tokuda et al. (2007, 2013) used a different technique for inference of the coupling function of multivariate electrochemical oscillations: see Sec. IV.2.3. Also, to capture the whole nature of the interaction of electrochemical oscillations (and not only the synchronization-related ones) the coupling functions were reconstructed in the full two-dimensional ($\phi_{1},\phi_{2}$) space i.e. not only for the one-dimensional diffusive coupling difference $\Delta\phi=\phi_{2}-\phi_{1}$ Blaha et al. (2011).

Of particular interest is a coupling function method that Miyazaki and Kinoshita (2006) applied for studying the interactions of Belousov-Zhabotinsky chemical oscillations. This class of reactions serves as a classical example of non-equilibrium thermodynamics, resulting in the establishment of a nonlinear chemical oscillator Strogatz (2001). The method infers the phase dynamics with diffusive coupling functions from the experimental phase time-series. The coupling function was already discussed in Sec. II.4, as shown in Fig. 8. Here we further review the interpretation and use of such coupling functions.

The inferred coupling function $q(\psi)$ is used to describe the mechanisms of the various synchronous states in two mutually coupled Belousov-Zhabotinsky reactors. The dynamics of the phase difference $\psi$ can be expressed as:

[TABLE]

where $Q(\psi)=q(\psi)-q(-\psi)$, $\Delta_{\omega}=\omega_{2}-\omega_{1}$ and $\varepsilon=\varepsilon_{12}-\varepsilon_{21}$. Then, by varying the frequency $\Delta_{\omega}$ mismatch and coupling strength $\varepsilon$, one can study and predict the occurrence of synchronization. Fig. 24 shows that with increase of $\varepsilon$, a pair of stationary solutions of Eq. (38) are obtained as the intersection points of $Q(\psi)$ and $\Delta_{\omega}/\varepsilon$ (the first horizontal line from the top in Fig. 24). There is one stable solution (solid circle) and in-phase synchronization is realized. If one increases $\varepsilon$ further, a new stable solution appears slightly below $\pi$ in addition to that corresponding to in-phase synchronization – the intersection with the second horizontal line from the top in Fig. 24. This corresponds to out-of-phase synchronization. Thus, a bistability between out-of-phase and in-phase synchronization can appear.

Kori et al. (2014) have performed a comprehensive theoretical analysis and experimental verification of phenomena in electrochemical oscillators, investigating the general occurrence of phase clusters in weakly, globally-coupled oscillators close to a Hopf bifurcation. The amplitude equation with a higher-order correction term, valid near a Hopf bifurcation point, is derived and it is used to calculate analytically the phase coupling function from given limit-cycle oscillator models. Such phase coupling functions, allowed the stability of phase clusters to be analyzed, as demonstrated on the Brusselator model.

Experiments on electrochemical oscillators have demonstrated the existence of three-cluster states near the Hopf bifurcation with negative coupling. Electric potentials were used to control the nature of the oscillations, and they were chosen initially such that the oscillators exhibited smooth oscillations near the Hopf bifurcation. Fig. 25(a) shows the current from one oscillator of each of the three clusters. The nearly balanced three-cluster state with configuration (25:20:19) is shown on a grid of 8 $\times$ 8 circles in Fig. 25(b). Phase response curves (Fig. 25(c)) and coupling functions (Fig. 25(d)) for these oscillators were found experimentally by introducing slight perturbations to the oscillations. The stability of the cluster states was determined, and it was found that the three-cluster state is the most stable, while four- and five-cluster states were observed also at higher potentials. Further increase in the potential resulted in complete desynchronization of the 64 oscillators.

V.2 Cardiorespiratory interactions

The heart and the lungs have physiological functions of great importance for human health and their disfunction may correspond to severe cardiovascular disease. Both organs are characterised by a pronounced oscillatory dynamics, and the cardiorespiratory interactions have been studied intensively using the theory and methods from the nonlinear coupled-oscillators approach Stefanovska et al. (2000); Stefanovska and Bračič (1999); Schäfer et al. (1998); Kenner et al. (1976).

The cardiorespiratory coupling functions are therefore a subject of great interest i.e. the mechanisms through which respiration influences the cardiac period and, in particular, how this relates to different states and diseases. The cardiorespiratory analysis performed with dynamical Bayesian inference Stankovski et al. (2012), as discussed in Sec. IV.2.2 and Fig. 16, revealed the form of the coupling functions in detail. The use of a changing respiration frequency in a linear (ramped) way showed that the form of the reconstructed coupling functions is in itself time-varying. Recently, the method was applied to the study of the effects of general anaesthesia on the cardiorespiratory coupling functions Stankovski et al. (2016). A similar form of the function was reconstructed for the awake measurements as in the previous studies, while its form was more varying and less deterministic for the state of general anaesthesia.

Dynamical Bayesian inference was used to study the effect of ageing on the cardiorespiratory interactions Iatsenko et al. (2013). Analyses were performed on cardiac and respiratory time series recorded from 189 subjects of both genders, aged from 16 to 90 years. By application of the synchrosqueezed wavelet transform, the respiratory and cardiac frequencies and phases were pre-estimated. By applying dynamical Bayesian inference to the phase time-series, measures such as synchronization, coupling directionality and the relative contributions of different mechanisms were then estimated.

The cardiorespiratory coupling function was thus reconstructed, and its time evolution and age dependence were assessed. Fig. 26 shows the time-averaged versions of the cardiorespiratory coupling functions typical of a younger and an older subject. Figs. 26(a) and (b) show the coupling functions of the heart dynamics $q_{h}(\phi_{h},\phi_{r})$. The form of the functions (especially noticeable in Fig. 26(a)) is changing mostly along the respiration phase $\phi_{r}$ axis, while it is nearly constant along the $\phi_{h}$ axis, indicating that this coupling is predominantly defined by the direct influence of respiration on the heart. In physiology, this modulation is known as respiratory sinus arrythmia (RSA). By comparing the coupling functions for the young and old subjects one can see a clear decrease of the RSA amplitude with age. It can also be noted that RSA remains the main stable contribution to the $q_{h}(\phi_{h},\phi_{r})$ coupling function, irrespective of age and that it survives after time-averaging. The coupling function from respiratory dynamics $q_{r}(\phi_{h},\phi_{r})$ shown in Fig. 26(c) and (d) was very low and seemed to be quite irregular and not age-dependent. From analysis of the time-variability of the form of the coupling functions it was observed that, in older people, the heart coupling function $q_{h}(\phi_{h},\phi_{r})$ becomes less stable in time, dominated by the highly time-variable indirect contributions. At the same time, the dynamics of the respiratory coupling function $q_{r}(\phi_{h},\phi_{r})$ did not seem to change with age, being irregular and unstable.

Kralemann et al. (2013a) looked at the cardiorespiratory coupling functions, as an intermediate result, in order to obtain the phase response curve of the heart and the perturbation inserted by respiration. They studied the respiratory and cardiac oscillations of 17 healthy humans while resting in an unperturbed state. The cardiac oscillation was assessed through two different observables – the electrocardiograph (ECG) and the arterial pulse signal. The idea of using two different observables is to demonstrate that an invariant phase can be obtained from each of them, describing a common inherent interaction between respiratory and cardiac oscillations.

By analysing the phase dynamics, first by estimating the protophases and transforming them into genuine phases, the cardiorespiratory phase coupling functions were reconstructed. Fig. 1(b) shows the reconstructed cardiorespiratory coupling functions $Q_{e}$ using the cardiac phase extracted from the ECG signal. The coupling functions estimated from the phases of the arterial pulse signal are shown in Fig. 27. The forms of the functions reveal the detailed mechanism through which respiration influences the cardiac oscillations i.e. the regions with high values of the function mean higher frequencies (acceleration), whereas low regions correspond to lower frequencies (deceleration) of the cardiac oscillations due to the respiratory influence. The existence of such cardiorespiratory coupling functions was tested statistically in respect of inter-subject and intra-subject surrogates.

The high similarity of the cardiorespiratory coupling functions obtained from the phases of the ECG observables Fig. 1(b) and those from the arterial pulse phases Fig. 27(b) demonstrates that the proposed method correctly identified the underlying interaction mechanism. This was achieved because the method was able to transform protophases from different observables into invariant phase dynamics from which a common form of the coupling functions was obtained. The minor differences in the form of the functions, as compared to the previous one with ageing (Fig. 27), are related to possible inter-subject variations and the different phase estimation approaches.

The similarity of the cardiorespiratory coupling functions among different subjects and between the two observable phases was further quantified with the similarity indices Kralemann et al. (2013a), as discussed in Sec. IV.3 and shown in Fig. 20. The similarity index that quantifies the correlation between the form of the functions, also proved very useful in assessment of the state of general anaesthesia Stankovski et al. (2016). It was found that the intersubject correlation similarity of the cardiorespiratory coupling functions, in comparison to the awake state, decreased with the onset of propofol-induced general anaesthesia, and to an even greater extent when sevoflurane was used.

V.3 Neural coupling functions

Neural states often manifest themselves as changes in brain electrophysiological activity, which emanates from the dynamics of large-scale cell ensembles oscillating synchronously within characteristic frequency intervals. Individual ensembles communicate to integrate their local information flows into a common brain network. One way to describe such an integration or communication is through cross-frequency coupling, an approach that has led to numerous studies elucidating the respective roles of cognition, attention, memory and anaesthesia Canolty et al. (2006); Jirsa and Müller (2013); Lakatos et al. (2008); Stefanovska (2007). Unlike these cross-frequency coupling strength approaches, the methods discussed here can assess the neural states through the computation of the coupling functions, hence describing the functional forms and mechanisms of individual cross-frequency interactions. In this way, one infers *neural cross-frequency coupling functions *Stankovski et al. (2017).

The methods for the reconstruction of coupling functions have been applied to electroencephalogram (EEG) recordings. The brainwave intervals, including the $\delta$ (0.8-4 Hz), $\theta$ (4-7.5 Hz), $\alpha$ (7.5-14 Hz), $\beta$ (14-22 Hz) and $\gamma$ (22-80 Hz), were first extracted from a single EEG channel recording. The phase was then extracted from each filtered time series, using for example the Hilbert transform or the synchrosqueezed wavelet transform. During this preprocessing procedure, particular care was taken to minimise overlap between the spectra of the intervals Lehnertz et al. (2014): overlaps of consecutive frequency intervals would result in overestimation of the corresponding phase-to-phase coupling. Dynamical Bayesian inference was then used to reconstruct the coupling functions from the multivariate five-phase oscillators. In a similar manner, dynamical Bayesian inference was applied to a study of neural interactions during epileptic seizures Wilting and Lehnertz (2015), though not for cross-frequency coupling.

As the brain is a highly complex system that can mediate a variety of functions from fixed structure Park and Friston (2013), the coupling relationships between the brainwaves can be different. One important coupling relation is the $\delta,\alpha\rightarrow\alpha$, as it has been found that the $\delta$-waves typical of deep sleep in adults can influence the $\alpha$-activity, which is related to the processing of information Feinberg et al. (1987); Jirsa and Müller (2013). Fig. 29(a) shows how the form of the $\delta,\alpha\rightarrow\alpha$ coupling functions varies in relation to their spatial locations on the head. It can be seen that the tridimensional waves propagate mostly in the $\delta$ dimension. This tendency can be seen better in Fig. 29(b), which shows the averaged coupling function. Its form depends predominantly on the direct delta oscillation, changing mostly along the $\phi_{\delta}$-axis. This reveals how and when within one cycle the $\delta$ oscillations accelerate and decelerate the $\alpha$ oscillations. Other coupling relationships could include for example the pairwise $\theta,\gamma\rightarrow\gamma$ and $\alpha,\gamma\rightarrow\gamma$, or the multivariate triplet $\theta,\alpha,\gamma\rightarrow\gamma$, as shown previously in Fig. 7.

Neural cross-frequency coupling functions were used recently to elucidate the mechanisms of general anaesthesia Stankovski et al. (2016). In fact, the analyses also included the cardiac and respiratory oscillations (in a sense integrating Secs. V.2 and V.3). Here, we review an important finding based on the neural $\delta$-$\alpha$ coupling functions. The study included 25 awake and 29 anaesthetized healthy subjects, of which 14 subjects were anaesthetized with the intravenous anaesthetic propofol and 15 subjects with the inhalational anaesthetic sevoflurane. The aim of the study was to determine if there are any differences in the interaction mechanisms in respect of the three states: awake; and anaesthetized with either propofol or sevoflurane.

Fig. 28 shows the group $\delta$-$\alpha$ coupling functions for the three states. The coupling functions for the awake resting, propofol, and sevoflurane states are evidently quite different from each other, both in the forms and strengths of the couplings. The $\delta$-$\alpha$ coupling function for the awake state has a relatively complex and varying form, and low amplitude. The coupling functions for propofol and sevoflurane are similar to each other and they look significantly different from those for the awake state. The sevoflurane coupling function has the largest coupling amplitude. Careful surrogate testing showed that the propofol and sevoflurane coupling functions are statistically significantly different from the corresponding surrogates. The qualitative form of the $\delta$-$\alpha$ coupling function has a sine-like form along the $\phi_{\delta}$-axis, while remaining nearly constant along the $\phi_{\alpha}$-axis. This implies that much of the $\delta$-$\alpha$ coupling comes from the direct contribution of the delta oscillation. The specific form of the function (e.g. Fig. 28(c)) reveals the underlying coupling mechanism, i.e. it shows that, when the delta oscillations are between $\pi$ and $2\pi$, the sine-wave coupling function is higher and the delta activity accelerates the alpha oscillations; similarly, when the delta oscillations are between [math] and $\pi$, the coupling function is decreased and delta decelerates the alpha oscillations.

The delta-alpha coupling has been linked to the coding mechanism of feedback valence information Cohen et al. (2009), non-REM sleep Bashan et al. (2012) and the eyes-closed state Jirsa and Müller (2013). The findings with anæsthesia are consistent with, and have further extended, these findings. The form of the $\delta$-$\alpha$ coupling functions (e.g. Fig. 28 (c)) indicates that the influence is direct modulation from delta to alpha, where the couplings are significantly stronger in anæsthesia than when awake. This showed that, once the subject is anaesthetized, delta activity influences the alpha oscillations by contributing to the reduction of information processing and integration.

V.4 Social sciences

In a recent application in social sciences, Ranganathan et al. (2014); Spaiser et al. (2014) have identified the coupling functions that capture interactions between social variables, employing a Bayes factor to decide how many interaction terms should be included in the model.

The work presents an interesting study of the relationship between democracy and economic growth, identifying nonlinear relationships between them. Economic growth is assessed through the GDP per capita (from the World Bank), while the level of democracy is gauged from the democracy index (from Freedom House) Ranganathan et al. (2014). It is well known that the GDP per capita and democracy are highly correlated: higher GDP implies more democracy. The linear Pearson correlation coefficient between the two variables is 0.571 ($p<0.01$). However, by use of coupling functions one can try to determine a more precise, causal relationship between the variables, revealing the underlying mechanism.

In its general form, the model considered is:

[TABLE]

where $\mathcal{D}$ denotes the democracy and $\mathcal{G}$ is GDP per capita. In this way the change of the variables $\mathcal{D}$ and $\mathcal{G}$ is represented with ordinary differential equations, even though the original data are discrete and one should really use difference rather than differential equations. Nevertheless, this approximation was used for mathematical simplicity. Further, the model can have some of the functions from a set of seventeen base functions of polynomial form, including reciprocal, quadratic and cubic terms. The main idea of the method is to select the optimal base functions thereby reducing the number of terms in the model.

The inference itself consists of two main steps. The first is an inferential fitting to obtain a model from the data, based on a maximum likelihood procedure, and involving multiple linear regression (similar to that discussed in Sec. IV.2.3). The second step uses a Bayesian Berger and Pericchi (1996) model selection procedure222We note that the Bayes factor uses the Bayesian probability theory too; however, as it is purely statistical procedure, it differs from the DBI as discussed in Sec. IV.2.2.. Here, the method decides how many interaction terms should be included in the model, i.e. it selects a subgroup of base functions of the seventeen polynomials available, after trying all possible combinations among them. Thus, the method punishes overly complicated models and identifies the models with the most explanatory power. Such procedures could benefit greatly if a surrogate testing procedure (see Sec. I.2.2) were used to determine whether the finally selected model is genuinely reliable.

The method was applied to model the interaction of democracy and GDP per capita for the years 1981-2006 for 74 countries. The resulting coupling functions for two selected models are shown in Fig. 30. The simplest model shown in Fig. 30(a) includes a coupling function with two terms $\dot{D}=0.11G^{3}-0.067D/G$, while the best fit five-function model Fig. 30(b) was given as $\dot{D}=0.77G^{3}+1.9D-0.85D/G-0.96DG-0.14D^{2}$. For the middle GDP, both of the coupling functions show dependences that closely relate to the linear dependence determined with the simple correlation coefficient. However, there were some nonlinear deviations from this, especially for very low and very high GDP. In particular, the threshold for very high GDP indicated that there is no significant improvement in democracy with further GDP growth.

Similarly, the best model for $G$ was inferred to be $\dot{G}=0.014+0.0064DG-0.02G$, which shows primarily that the GDP is growing at a constant rate, but in addition demonstrates that it is positively affected by democracy interacting with GDP, and that the growth is self-limiting at high levels of GDP.

Finally, we point out that the method was further applied to investigate the interactions of other social variables, such as the case of interactions between democracy, development and cultural values Spaiser et al. (2014). These works could benefit significantly from further coupling function assessment and analysis.

V.5 Mechanical coupling functions

Mechanical clocks and oscillators provide an important cornerstone in the study of interactions and synchronization phenomena, starting from the earliest observations of the phenomenon in pendulum clocks by Huygens (1673), up to the more comprehensive and detailed studies based on current methods Kapitaniak et al. (2012).

Kralemann et al. (2008) describe an experiment using two coupled mechanical metronomes for the analysis of coupling functions. The metronomes were placed on a rigid base and the coupling through which they interact and influence each other was achieved by connecting them with an elastic rubber band – Fig. 31(top). A digital camera was used for acquiring the data, from which the oscillatory signals were extracted. Coupling functions were determined for three different experimental conditions, when: (i) the pendulums of the metronomes were linked by a rubber band; (ii) the pendulums were linked by two rubber bands; and (iii) the metronomes were uncoupled.

From the extracted signals, the Hilbert transform protophases were first estimated, and then transformed to genuine phases. The coupling functions were then reconstructed with a fitting procedure based on kernel smoothing. Fig. 31 shows the results for the three cases. By comparison of the coupling functions in the case of one rubber coupling Fig. 31(a),(b) with the coupling of two rubber coupling Fig. 31(c),(d), one can see that the form is very similar, while the coupling strength is slightly higher for the case of two-band coupling. The coupling functions is of a complex form, changing along both axes, thus reflecting the bidirectional influence and contribution within the couplings. The methodology correctly detects extremely weak coupling functions Fig. 31(e),(f) for the case of no explicit coupling. Future developments of this work could benefit from comparison of the extracted coupling functions with the actual mechanics of the coupled metronomes, as well as from validation of the weak coupling regime for better justification of the use of phases.

V.6 Secure communications

The findings that the cardiorespiratory coupling function can be decomposed into a number of independent functions, and that the latter can have a time-varying nature Stankovski et al. (2012) and Sec. V.2, inspired the creation of a new class of secure communications characterized by high efficiency and modularity Stankovski et al. (2014).

The protocol (Fig. 32) starts with a number of information signals coming from different channels or communications devices (e.g. mobile phone, sensor networks, or wireless broadband) needing to be transmitted simultaneously. Each of the signals $s_{i}$ is encrypted in an amplitude coupling function; i.e. they serve as scaling parameters in the nonlinear coupling functions between two self-sustained systems in the transmitter. The coupling functions constitute the private key and, in principle, have an unbounded continuum of possible combinations. Two signals, one from each system, are transmitted through the public channel. At the receiving end, two similar systems are enslaved, i.e. completely synchronized, by the two transmitted signals. Finally, by applying time-evolving dynamical Bayesian inference (as discussed in Sec. IV.2.2) to the reconstructed systems, one can infer the model parameters and decrypt the information signals $s_{i}$.

This application is similar to that where amplitude coupling functions are reconstructed from data. The coupled systems are multidimensional and may be, e.g. chaotic Lorenz or Rössler systems. The great advantage is that the model of the coupled systems is known exactly on the side of the receiver where the inference is performed. Thus the problem of not knowing the amplitude model and its dimensionality does not exist. The main task of the decryption lies in inferring the time-evolution of the parameters.

The protocol can encrypt multiple signals simultaneously as time-evolving parameters. Each of them scales one of the coupling functions, which are nonlinear and mutually linearly-independent. Thus the method inherently allows for multiplexing i.e. simultaneous transmission of multiple signals. Another property of the protocol is that it is highly resistent to channel noise. This is because the dynamical Bayesian inference is performed for stochastic dynamics, so that the method is able very effectively to separate the unwanted noise from the deterministic dynamics carrying the messages.

The method is demonstrated on two bidirectionally coupled chaotic Lorenz systems. Ten information signals $s_{1},\ldots,s_{10}$ are encrypted with ten coupling functions, as indicated on the ordinate axis in Fig. 33. The choice of the particular forms of the coupling functions prescribe the private key. After the transmission, the systems are reconstructed on receiver side and, by application of Bayesian inference, the information signals are reconstructed. Fig. 33 shows good agreement of the original with the decrypted bits.

VI Outlook and Conclusion

VI.1 Future directions and open questions

Coupling functions have been studied as early as some of the first theoretical works on interactions and they remain a very active field of research that is attracting increasing interest from the scientific community. They bring a certain complexity in understanding, but at the same time they also illuminate, and provide deeper insight into, the interaction mechanisms. As such, coupling functions pose many open questions and there still remain many related aspects that are not well understood. Below, we discuss some of the open questions and current possibilities for further developments related to coupling functions.

VI.1.1 Theory

The theoretical development of coupling functions will lead to a better understanding of the mechanisms responsible for the resulting overall dynamics, and they may help to incorporate seemingly different models into a general overall framework. Future theoretical studies need to identify the classes of coupling functions that lead to particular physical effects. In doing so, one needs to determine if there are some classes of functions which demonstrate unique characteristics, and more importantly whether some particular functions lead to common effects. As a consequence, these can then lead to classes of functions to be used in engineering, for controlling or predicting the outcome of the interactions.

The experimental results suggest two important directions for the theoretical development. First, coupling functions can be nonautonomous. Secondly, coupling functions can lead to the coexistence of attractors.

The theory of nonautonomous dynamical systems has gained recent interest, mainly in relation to finite-time bifurcations Kloeden and Rasmussen (2011). These mathematical developments will play a major role in the theory of coupling functions. The theoretical studies should include systematic and comprehensive descriptions of the different classes of coupling function, including the nonautonomous case.

The coexistence of attractors has gained considerable attention. As experimental results show, coupling functions are important, especially in relation to network structure and the effect on the basin of attraction, e.g. the basin of attraction of synchronization. Recent results have shown that the roles of coupling function and network structure can be nontrivial Menck et al. (2013). We still need new theoretical methods to tackle problems associated with the involvement of coupling functions in the coexistence of attractors. These questions are intimately related to the stability of the system and will play a role in important practical applications, such as to electrical power-grids.

VI.1.2 Methods

The future development of coupling function methods needs to take into account all the advantages and pitfalls of current methods, e.g. as outlined in the critical comparison Sec. IV.2.6. So far, all the coupling functions methods have been applied to pairs of coupled systems, or small to medium size networks. New methods should allow applications to the more prevalent large-scale networks. In line with this, they should aim to achieve faster calculations so as to facilitate the ever-growing demand for extensive computation.

There is also a room for improvements of the amplitude coupling function methods. The search for more generally applicable amplitude models remains open. These should be as general as possible, or at least general enough for specific sub-systems.

There is a need to overcome the problem of inferring coupling functions from systems that are highly synchronized and coherent in 1:1 frequency ratio. Currently, this is a common deficiency in all of the available methods for coupling (function) reconstruction. A possible direction for solving this issue could lie in the use of perturbation, for example starting from different initial conditions or with some other form of temporary deviation from the highly coherent state.

To enhance the inference of more general interactions, efforts are needed to develop and design robust methods for distinguishing direct from indirect couplings, better surrogates and null models for determining significant couplings and interpretation of the couplings in high-dimensional networks.

VI.1.3 Analysis

We also point to the need for further development in coupling function analysis: although the basic coupling function theory and methods are relatively well developed, there is a pressing need for measures able to exploit the computed coupling functions to better effect. The current tools, including similarity analysis and coupling function decomposition, are very useful, but there is a clear need for development of even more analysis tools. The task here is to find better and more systematic ways of quantifying and describing the form of the coupling functions and the other functional characteristics unique to coupling functions.

The development of such methods for analysing and characterizing coupling functions could be linked to the mathematical theory of functional analysis. To date, this theory including the main concepts in vector spaces, and measures of mappings between the functions, have not yet been fully exploited in relation to coupling functions.

VI.1.4 Integration theory-applications

We emphasize that further interplay between theory and applications for the development of coupling functions is still needed. Although the applications usually take into account theoretical developments, recent experimental findings have not yet been properly addressed by the theory. For example, the theoretically most studied form of coupling functions is that for diffusive coupling, which includes the state or phase difference as an argument. The latter is mostly used because it provides convenient solvable solutions (see e.g. history, Sec. II.2). On the other hand, the coupling decomposition experiments have shown that it is the direct coupling function that often predominates in reality, especially in biological oscillatory interactions. Hence, further theoretical studies are needed to establish the phenomena and the nature of interactions for direct coupling functions. Such theoretical investigations can usefully be performed numerically in cases where the relevant model cannot be solved analytically.

VI.1.5 Applications

Coupling functions have universal implications for all interactions between (dynamical) systems. As such, they can describe mechanisms operating between systems that are seemingly of very different natures. We have reviewed a number of important applications, including for example chemical, biomedical, mechanical, social and secure communications; however, the unique features of coupling functions promise even further application in these and in other fields. We outline below some foreseeable directions for new applications, notwithstanding that many others are also possible.

Recently, there has been a significant interest and developments in the study of interactions and synchronization in power grids Rubido (2015); Rohden et al. (2012). To ensure a reliable distribution of power, the network should be highly controllable and synchronized. It is therefore very important that the state of synchronization should be highly stable (i.e. deep in the Arnold tongue), so as to ensure that small disruptions and glitches will not interrupt the function of the network. Coupling functions should be investigated in order to establish how to design and engineer a persistent Pereira et al. (2014) and very stable state of synchronization.

Similar problems occur with the control, synchronization and optimization in transport grids Rodrigue (1999); Albrecht (2004), for example in a rail network. In such cases, of vital importance are the dynamical and the time-varying events. The developed methods and theory for time-varying coupling functions could be of great use in these applications.

Increasing the scale of the networks often leads to higher-level organization, including networks of networks and multilayer networks Stern (2013); Kivelä et al. (2014). In such high-dimensional spaces, a variety of different physical effects can be observed, e.g. synchronization, chimeras, and clustering. The coupling functions of different levels and layers could provide deeper insight into the functions or subfunction integration of the networks.

Coupling functions have been found very useful in studying the interactions between macroscopic physiological systems, such as those between the cardiorespiratory and neural systems reviewed above. Further coupling function investigations will probably be developed between different oscillations in integrated network physiology Bashan et al. (2012); Stefanovska (2007). In a similar way, coupling functions between microscopic physiological organizations could be developed. The latter could explore coupling functions between cells including, for example, the oscillations of neurons or stem cells Murthy and Fetz (1996); Eytan and Marom (2006); Méndez-Ferrer et al. (2008); Jackson et al. (2001).

VI.2 Conclusion

In recent years, the investigation of coupling functions has developed into a very active and rapidly evolving field. Their study and use have brought huge progress in the understanding of the mechanisms underlying the diverse interactions seen in nature. The enterprise has now reached a critical mass, offering increased potential for new and important discoveries, and in this way the topic has attained a substance and unity justifying the present review.

The concept of the function in the coupling functions is perhaps its most important characteristic. Yet, precisely because of being a function, it is inevitably harder to interpret, assess and compare than is the case for quantitative measures such as the coupling strength. In attempting to integrate and pull together existing knowledge about coupling functions, therefore, we have tried to organize, explain and, as far as possible, to standardize their description in the hope of making them more generally accessible and useful.

Interactions underlie many important phenomena and functions of the systems found in nature, and it is of great importance to be able to describe and understand the mechanisms through which the interactions occur. Coupling functions are opening up new perspectives on these interactions and we envisage that they will catalyze increased research activity on coupled dynamical systems and their interactions in the future.

Acknowledgements.

We gratefully acknowledge valuable and continuing discussions with Ralph Andrzejak, Peter Ashwin, Murilo Baptista, Miroslav Barabash, Chris Bick, Philip Clemson, Andreas Daffertshofer, Andrea Duggento, Jaap Eldering, Deniz Eroglu, Dmytro Iatsenko, Viktor Jirsa, Peter Kloeden, Ljupco Kocarev, Thomas Kreuz, Juergen Kurths, Jeroen Lamb, Gemma Lancaster, Klaus Lehnertz, Zoran Levnajić, Dmitry Luchinsky, Hiroya Nakao, Milan Paluš, Thomas Peron, Spase Petkoski, Arkady Pikovsky, Antonio Politi, Alberto Porta, Martin Rasmussen, Francisco Rodrigues, Michael Rosenblum, Rajarshi Roy, Björn Schelter, Ivana Stankovska, Sinisa Stojanovski, Yevhen Suprunenko, Peter Tass, Valentina Ticcinelli, and Dmitry Turaev. The work was supported by the Engineering and Physical Sciences Research Council UK (grant numbers EP/100999X1, EP/M015831/1 and EP/M006298/1), by the ITN COSMOS programme (funded by the EU Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie Grant Agreement No. 642563), by Action Medical Research UK (grant number GN1963), by Lancaster University Department of Physics, by the Fundação de Amparo à Pesquisa do Estado de São Paulo Brazil (grant number FAPESP 2013/07375-0), by the Institute of Pathophysiology and Nuclear Medicine, Faculty of Medicine, UKIM Skopje, Macedonia and by the Slovenian Research Agency (Program No. P20232).