Coupling functions: Universal insights into dynamical interaction mechanisms
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 globallycoupled 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.
pacs:
05.45.Xt, 05.45.a, 05.45.Tp, 02.50.TtContents
 I Introduction
 II Basic Concept of Coupling Functions
 III Theory
 IV Methods
 V Applications and Experiments
 VI Outlook and 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 coevolution 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 manmade 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 (multipleshooting) 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 timevariability and for studying the evolution of the cardiorespiratory coupling functions with age; in neuroscience for revealing the crossfrequency 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 nonidentical 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 lowdimensional 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 quasisynchronized 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 dimensionallyreduced 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 BarEli (1985); Mirollo and Strogatz (1990); Prasad (2005); SuárezVargas 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 informationtheoretic 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 weaklycoupled systems, the methods will usually detect nonzero 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 “zerolevel” 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 nonoscillatory, nondynamical, 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:
where is a differentiable vector field with 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:
(1)  
(2) 
where are vector fields describing the isolated dynamics (perhaps with different dimensions) and 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:
(3)  
(4) 
where the functions 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 wellstudied 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):
(5) 
where are the phase variables of the oscillators, are their natural frequencies, 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 . In this case, from and Eqs. (5) we can express the interaction in terms of as:
(6) 
Synchronization will then occur if the phase difference is bounded, i.e. if Eq. (6) has at least one stableunstable pair of solutions Kuramoto (1984). Depending on the form of the coupling function, in this case the sine form , 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 .
Fig. 4 illustrates schematically the connection between the coupling function and synchronization. An example of a synchronized state is sketched in Fig. 4(a). The resultant coupling strength has larger values of the frequency difference at certain points within the oscillation cycle. As the condition is fulfilled, there is a pair of stable and unstable equilibria, and synchronization exists between the oscillators. Fig. 4(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. 4(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 saddlenode bifurcation, leading to a stable synchronization.
The coupling functions of real systems are often more complex than the simple sine function presented in Fig. 4(a) and (b). For example, Fig. 4(c) also shows a synchronized state, but with an arbitrary form of coupling function that has two pairs of stableunstable points. As a result, there could be two critical coupling strengths ( and ) and either one, or both, of them can be larger than the frequency difference , 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. 4(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 limitcycle 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:
(7) 
Here, is the influence function through which the second oscillator affects the first, while the sensitivity function 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 , , 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 and 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 epochmaking 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), Josephsonjunction 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 socalled phase oscillators, and it studies synchronization phenomena in a large population of such oscillators Kuramoto (1984). By setting out a meanfield 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 alltoall diffusive coupling rather than local coupling, took the meanfield limit, introduced a random frequency distribution, and assumed that a limitcycle orbit is strongly attractive Kuramoto (1975). As already mentioned, Kuramoto’s coupling function was a sinusoidal function of the phase difference:
(8) 
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. 4).
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 limitcycle, one can apply phase reduction procedures (see Sec. III.2 for further theoretical details) which systematically approximate the highdimensional dynamical equation of a perturbed limit cycle oscillator with a onedimensional reducedphase equation, with just a single phase variable representing the oscillator state Nakao (2015). In uncoupled or weaklycoupled 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 phasecoupled:
(9) 
Our aim is to describe the effect of the coupling function through which the first oscillator influences the second one. From the expression for in Eq. (9) one can appreciate the fundamental role of the coupling function: is added to the frequency . Thus changes in the magnitude of will contribute to the overall change of the frequency of the second oscillator. Hence, depending on the value of , 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. 5. Because in real situations one measures the amplitude state of signals, we explain how the amplitude signals (Fig. 5(a) and (d)) are affected depending on the specific phase coupling function (Fig. 5(b) and (c)). In all plots, time is scaled relative to the period of the amplitude of the signal originating from the first oscillator (e.g. ). For convenient visualisation of the effects we set the second oscillator to be fifteen times slower than the first oscillator: . The particular coupling function presented on a grid (Fig. 5(b)) resembles a shifted cosine wave, which changes only along the axis, like a direct coupling component. Because all the changes occur along the axis, and for easier comparison, we also present in Fig. 5(c) a averaged projection of .
Finally, Fig. 5(d) shows how the second oscillator is affected by the first oscillator in time in relation to the phase of the coupling function: when the coupling function is increasing, the second oscillator accelerates; similarly, when decreases, decelerates. Thus the form of the coupling function 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 (). 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. 6. The sinusoidal function of the phase difference from the Kuramoto model exhibits a diagonal form in Fig. 6(a), while the influencesensitivity product function of Winfree model is given by a more complex form spread differently along the twodimensional space in Fig. 6(b). Although these two functions differ from those in the previous example (Fig. 5), 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) higherdimensional space are necessarily observable through the lowdimensional space of the measurements. Frequently, a delay embedding theorem Takens (1981) is used to reconstruct the multidimensional dynamical system from data. In real application with nonautonomous and nonstationary 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é limitcycle oscillators. In the autonomous case, each of them is given by the polar (radial and angular ) coordinates as: and . In this way, a Poincaré oscillator is given by a circular limitcycle and monotonically growing (isochronous) phase defined by the frequency parameter. In our example, we transform the polar variables to Cartesian (state space) coordinates , , and we set unidirectional coupling, such that the first (autonomous) oscillator:
(10) 
is influencing the state of the second oscillator through the quadratic coupling function :
(11) 
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 , and we set a relatively high coupling strength .
The description of the amplitude coupling function is illustrated schematically in Fig. 7. In theory, the coupling function has four variables, but for better visual illustration, and because the dependence is only on , we show it only in respect to the two variables and i.e. . The form of the coupling function is quadratic, and it changes only along the axis, as shown in Figs. 7(b) and (c). Finally, Fig. 7(d) shows how the second oscillator is affected by the first oscillator in time via the coupling function: when the quadratic coupling function is increasing, the amplitude of the second oscillator increases; similarly, when decreases, decreases as well.
The particular example chosen for presentation used a quadratic function ; other examples include a direct linear coupling function e.g. , or a diffusive coupling e.g. 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 preestimated most probable models Berger and Pericchi (1996), but without including explicit assessment of the coupling functions. Due to the multidimensionality 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 reactiondiffusion 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 , 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 smallscale 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 crosscoupling function, for example of the form in the dynamics of the first oscillator . 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. 10. Comparing the true (Fig. 10 left) and the inferred effective (Fig. 10 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 and which is tested in respect of the condition , for 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 crossfrequency coupling in neuroscience Jensen and Colgin (2007). Recently, neural crossfrequency coupling functions were extracted from multivariate networks Stankovski et al. (2015) (see also Sec. V.3). The network interactions between the five brainwave oscillations , , , and were analysed by reconstruction of the multivariate phase dynamics, including the inference of triplet and quadruplet coupling functions. Fig. 11 shows a triplet coupling function of how the and influence 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 to .
Very recently, Bick et al. (2016) have shown theoretically that symmetricallycoupled phase oscillators with multivariate (or nonpairwise) coupling functions can yield rich dynamics, including the emergence of chaos. This was observed even for as few as oscillators. In contrast to the KuramotoSakaguchi equations, the additional multivariate coupling functions mean that one can find attracting chaos for a range of normalform 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 welldefined 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 nonuniqueness of the estimate. That is, we project the function (which is infinitedimensional) onto a finitedimensional 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.
Type of CF  Model  Meaning  Reference 

Direct  unidirectional influence  Aronson et al. (1990)  
Diffusive  dependence on state difference  Kuramoto (1984)  
Reactive  complex coupling strength  Cross et al. (2006)  
Conjugate  permutes the variables  Karnatak et al. (2007)  
Chemical synapse  is a sigmoidal  Cosenza and Parravano (2001)  
Environmental  given by a differential equation  Resmi et al. (2011) 
In the literature, authors often speak of the commonlyused 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 satisfying a local condition 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 the coupling is called repulsive Hens et al. (2013). Chemical synapses are an important form of coupling where the influences of and 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 for , so that the variables are considered as external fields driving the equation. Its solution is taken as the coupling function and, for , 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 limitcycle 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 BelousovZhabotinsky chemical oscillator has been reconstructed Miyazaki and Kinoshita (2006) with the help of a method for the inference of phase dynamics. Fig. 12 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 and abruptly increases at a larger , with its minimum and maximum at around and , 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 spikelike 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 oneaxis), 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 selfcoupling describes the inner dynamics of an oscillator which results from the interactions and has little physical meaning. Directcoupling describes the influence of the direct (unidirectional) driving that one oscillator exerts on the other. The last component, indirectcoupling, often called commoncoupling, 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 BelousovZhabotinsky reactors can be explained in terms of the coupling function as illustrated by the examples given in Fig. 12 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 orderdisorder 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 of two coupled oscillators are reconstructed for three characteristic cases, as shown in Fig. 13(a)(c), left panels. The right panels in Fig. 13 show the corresponding odd (antisymmetric) part of the coupling functions , which is important for determination of the synchronization. The coupling functions Fig. 13 (a)(c) have predominantly positive values, so the interactions contribute to the acceleration of the affected oscillators. The first coupling function Fig. 13(a) for smooth oscillations has a sinusoidal which can lead to inphase synchronization at the phase difference of . The third case of relaxation oscillations Fig. 13(c) has an inverted sinusoidal form , leading to stable antiphase synchronization at . The most peculiar case is the second one Fig. 13(b) of relaxation oscillations, where the odd coupling function takes the form of a second harmonic () and both the inphase () and antiphase () 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 twocluster synchronized state. Experiments in a parameter region between these states, in which bistability is predicted, are shown in Fig. 14. A small perturbation of the stable onecluster state (left panel of Fig. 14) yields a stable twocluster state (right panel of Fig. 14). 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 phaseresponse curves (PRCs) or spiketime responsecurves (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 spiketime responsecurves of individual neuronal cells. Results from these singlecell experiments were then used to predict the multicell network behaviors, which were found to be compatible with previous modelbased 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 fullycoupled oscillators whereas in the weak coupling we can reduce the analysis to lowerdimensional 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
(12) 
where 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
(13)  
The problem proposed by Smale was to find (if possible) a coupling function (positive definite matrix) 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 . The following example illustrates the main ideas. Consider
(14) 
where . Note that all the eigenvalues of have negative real parts. So the origin of the system Eq. (14) is exponentially attracting.
Consider the coupling function to be the identity
For the origin is globally attracting; the uniform attraction persists when is very small, and so the origin is still globally attracting. However, for large values of the coupling 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 diffusiondriven oscillators. Further oscillations in originally passive systems have been reported in spatially extended systems GomezMarin 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 BarEli (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
where
So, each isolated system has a limit cycle of amplitude and a frequency of . Note that the origin 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 . The equations have the same form as Eq. (13) with being the identity.
The effect can be better understood in terms of phase and amplitude variables. Let , be the amplitudes and the phases of and , respectively. We consider which captures the main causes of the effect, as well as the phase difference . Then the equations in these variables can be well approximated as
(15)  
(16) 
The conditions for oscillation death are a stable fixed point at 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. 4. There is a fixed point if and . We will assume that which implies that, when the fixed point exists, .
Next, we analyse the stability of the fixed point . This is determined by the linear part of Eq. 15. Hence, the condition for stability is
Using the equation for the fixed point we have . Replacing this in the stability condition we obtain . 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, , then there are minimum and maximum coupling strengths for oscillation death
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
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 , and discuss the key mechanism for synchronization. Note that there are natural coordinates to analyze synchronization
These coordinates have a natural meaning. If the system synchronizes, and with . Hence, we refer to as the coordinate parallel to the synchronization subspace , and to as the coordinate transverse to the synchronization subspace, as illustrated in Fig. 15.
The synchronization analysis follows two steps: (i) Obtaining a governing equation for the modes 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 is small. Then we can obtain a linear equation for by neglecting the high order terms proportional to . Noting that , using Eqs. (13) for and , and expanding in a Taylor series we obtain
(17) 
where is the Jacobian evaluated along a solution of .
(ii) Coupling function to provide damping. The term coming from the coupling now plays the role of a damping term. So we expect that the coupling will win the competition with and will force the solutions of to decay exponentially fast to zero. To see this, we observe that the first term
(18) 
depends on the dynamics of alone. Typically, for . Now, to obtain a bound on the solution of Eq. (17), we consider the ansatz and notice that differentiating we obtain Eq. (17). Hence
from the growth behaviour of the disturbance 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 decay to zero. For coupling strengths
the oscillators synchronize.
Meaning of : This corresponds to chaotic behaviour in the synchronization manifold. If then, will be the Jacobian along a solution of So depends on the dynamics on the synchronization manifold. If and the solutions are bounded, the dynamics of the synchronized system is chaotic. Roughly speaking, 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. 15 for the chaotic Lorenz attractor. The number 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.
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 . This is the case, when the coupling function eigenvalues have positive real parts. Pereira et al. (2014).

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.

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 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
(19) 
Since is positive definite we obtain where is the smallest eigenvalue of . The global stability of the system can be obtained by constructing a Lyapunov function . The system will be stable if is positive and its derivative is negative. This system admits a quadratic Lyapunov function Indeed, taking the derivative
Hence all solutions of Eq. (19) will converge to zero exponentially fast. Next, consider the coupled system
where by the mean value theorem we obtain
Because we did not Taylorexpand the vector fields, the equation is globally valid. Assuming that the Jacobian is bounded by a constant , we obtain .
Computing again the Lyapunov function for the coupled system (including the vector fields) we obtain
(21) 
The system will synchronize if is negative. So synchronization is attained if
Again the critical coupling has the same form as before. The coupling function came into play via the constant , and instead of we have . Typically, is much larger than . 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 Taylorexpand the system to obtain
(22) 
in just the same form as before. Note however that the trick we used previously, by defining , is no longer applicable. Indeed, we use the ansatz to obtain and, since and do not commute,
the ansatz cannot be used. Thus we need a better way forward. So in the same way as we calculated the expansion rate for , 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 to obtain a function
This function is called the Master Stability Function (MSF). We will extract the synchronization properties from . As we already discussed, the solutions of Eq. (19) will behave as
Now note that (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 such that
and for which will become negative, and will converge to zero.
Global versus local results. In the global analysis, the critical coupling depends on , 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 . The critical coupling was then obtained by analysing the sign of the Lyapunov exponent. Generically, . 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 can be nonlinear. In this case, we consider the Jacobian
Moreover, we say that belongs to the Lyapunov class if there are positive matrices and such that
Whenever the matrix is in the Lyapunov class we can construct the Lyapunov function algorithmically.
Coupling Function  Class  Technique  Global  Persistence 

+ve definite  Lyapunov  LF  Yes  Yes 
+ve definite  Lyapunov  LF  No  Yes 
differentiable  generic  LF  
differentiable  generic  MSF  No  No 
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 weakcoupling 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 welldeveloped. 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 reactiondiffusion 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
(23) 
where , has a uniformly exponentially attracting periodic orbit with period , that is, . 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 , . We can also reparameterize time such that that phase increases uniformly along the orbit . That is, the phase is uniform frequency equal to unity. By the chain rule we then have
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 ). To this end we need to extend the phase 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 is defined only along the orbit . Our first step is to extend to a neighbourhood of . Since the periodic orbit is exponentially and uniformly attracting, it will attract an open neighbourhood of . We call this set the basin of attraction of the periodic orbit. Note that every initial point in the basin of attraction of the orbit will converge to the orbit. Hence, we have a such that
where is the solution of the system with initial condition . For each initial point in the basin of attraction of we can assign a unique point in the orbit . This is called the asymptotic phase.
Isochron. For each value of phase in the orbit 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 . So points in the isochron have the same value of phase and will move at the same speed. See Fig. 16 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 being perturbed
where is the phase of the external influence, and the influence is a periodic on . The weak coupling implies .
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 , there is a new attracting periodic orbit close to the orbit , the difference between the orbits being of order . Moreover, 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 . This neighborhood contains the new orbit as we consider small . So, calculating the phase along , by the chain rule we obtain
But by construction in a neighborhood of . Because is distant from we can expand both and and evaluate them along at the expense of a perturbation of order . It is standard to denote . In this setting we have
and we have successfully reduced the problem to the phase of the unperturbed orbit . In general terms, we study problems of the following type
(24) 
The main insight was to obtain the coupling function in terms of how the phase of the unperturbed orbit behaves near the orbit . We performed the following steps:
 (1) Phase sensitivity of the unperturbed system.

Once we have the asymptotic phase, we can use it as the new phase variable , 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
where the gradient is evaluated along the orbit .
 (2) Obtain coupling function by .

For this step, we need to take the inner product of with the perturbation . When studying collective phenomena will contain fast variables and slow variables. Typically, only the slow variable are of interest, so we will average over the fast variables.
Meaning of . 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 and the net effect will be the sum of and . Notice, that the linearity is only in terms of the perturbations. The equation itself is nonlinear in the phase variable . The linearity with respect to perturbations is because we have discarded all nonlinear terms and terms of order (by computing 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.
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
(25) 
By inspection, it is clear that is an unstable point and the system has an attracting periodic orbit of radius 1. This can be better seen by changing to polar coordinates and using which, we obtain
(26) 
The orbit corresponds to and is shown in Fig. 16.
Asymptotic phase: The phase as defined in the orbit has a constant frequency equal to unity. Along the orbit , we therefore have (by inspection of the equations). For points outside the orbit, however, this is no longer true. The asymptotic phase will fix this issue because the points then move at the same speed as the corresponding points in the orbit, so that for points outside the orbit.
Because of the symmetry ( does not depend on ) we can use the ansatz
where we aim to find the function . Differentiating we obtain and, using the isochron’s properties together with the equations for and , we obtain so Since we want to extend the phase continuously from the orbit, if then . We choose the constant Therefore,
and we can define the isochron with asymptotic phase In Fig. 16(top) we show four level sets of the asymptotic phase corresponding to and .
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 we obtain
(27) 
which is valid, not only along the orbit 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 .
We can now readily take the gradient (in polar coordinates), yielding Along the unperturbed orbit we have so that the phase sensitivity functions are
Next we obtain the coupling function.
Obtaining the coupling function of external forcing. Now we consider the system being forced at frequency:
(28) 
where . We obtain the coupling function through the isochron. We now justify in detail why we have discarded the corrections in .
We compute the equation for the phase dynamics. Note that, by chain rule . Using Eq. (27) and evaluating the gradient along the orbit , we obtain
For small , we know that the difference between and is of order , so we can replace the gradient along the perturbed orbit and unperturbed orbit with corrections of order (because is already multiplying the function). Hence,
Synchronization and coupling function. The main idea is that the coupling function can help in adjusting the frequency of the system to the frequency of the forcing. As we discussed above, we will neglect the terms . Introducing the phase difference
and considering we obtain
If is of order then the dynamics of will be slow in comparison with the dynamics of . Roughly speaking, for each cycle of we have cycles of . Because, the dynamics of is faster than that of , we use the averaging method to obtain the coupling function
where is the period of as a function of . Note that, for our result, is sinusoidal so that by integrating over while keeping fixed, we obtain . Hence we obtain the dynamics in terms of the phase difference
(29) 
which is exactly the equation shown in Fig. 4.
We are now ready to study collective phenomena between the driving and the system. For instance, the system will phaselock with the driving when . In this case, the oscillators will have the same frequency. Moreover, because is a periodic function, the fixed point will exist only when .
Higher order n:m phase locking. Our assumption is that , so that Eq. (29) for the phase difference is a slow variable. It may happen that gives rise to a slow variable. In such cases, we perform the same analysis for and further information on the higher order phase locking can be obtained Ermentrout (1981).
iii.2.4 Phase response curve
The phase sensitivity function 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 socalled 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 . 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 is . Every point on the orbit can be uniquely described by the phase . A small perturbation applied at a phase can cause the phase to complete its full cycle at a time . The normalized phase difference between the cycles is
Note that the PRC depends on the phase at which the small perturbation was applied, that is, . This is the socalled 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 , and it is defined as the PRC normalized by the amplitude of the perturbation
Indeed, the isochrons and the PRC are closely related. If a point moving along the orbit is instantaneous and the perturbation is small, the point will land on an isochron, which tells us the new phase 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, pulsecoupled 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 pulsecoupled 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 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). is a vector and, in our example in Sec. III.2.3, the norm of was proportional to . This is typical of Hopf bifurcations. In Fig. 17 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 . So we only need to compute the derivative with respect to . Thus, ^{1}^{1}1We are abusing the notation by using 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 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. 17. 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).
Bifurcation  

SNIPER  
Hopf  
Homoclinic  
Integrate and Fire  
Leaky Integrate and Fire 
iii.3 Globally coupled oscillators
Now suppose that we have coupled oscillators
(30) 
We assume that, when they are uncoupled , each system has an exponentially attracting periodic orbit. So the dynamics of the uncoupled system occurs on a torus that is exponentially attracting. Moreover, we also assume that is close to . 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
where each oscillator has its own period , and is the coupling function describing the influence of the th oscillator on the th oscillator. Here
Note that is independent of the index because we assumed that ’s are all close to . 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
This is by far the beststudied 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 . The model is then written as
If the oscillators are identical then any small coupling leads to synchronization (the phases will converge to the same value). If the distribution of natural frequencies is broad, then at a critical coupling 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
and to rewrite the equations in terms of the parameters and (which are now meanfield parameters)
Taking the limit , we can write the model in terms of selfconsistent equations Kuramoto (1984).
When the distribution of the natural frequencies is an even, unimodal, and nonincreasing function, and the coupling is weak, the incoherent state is neutral Strogatz and Mirollo (1991), but the order parameter vanishes (at a polynomial rate) if is smooth Fernandez et al. (2014). Moreover, on increasing the coupling, the incoherent solution bifurcates for
A systematic review of the critical coupling, including bimodal 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 . This means that the coupling function will have a phase shift
This coupling function is called the SakaguchiKuramoto coupling Sakaguchi and Kuramoto (1986). The slight modification can lead to nonmonotonic behaviour of synchronization Omel’chenko and Wolfrum (2012). For certain unimodal frequency distributions , 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 SakaguchiKuramoto coupling
This coupling function generalizes the standard SakaguchiKuramoto model as it allows for different contributions of oscillators to the mean field, on account of the phase shifts and and coupling factors . In turn, the mean field acts on each oscillator differently. This scenario is tractable in terms of the selfconsistency 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 powergrids.
Higher harmonics. In the previous discussion, the coupling function contained one harmonic. . Depending on the underlying bifurcation we must now include further, higherorder, Fourier components,