Coupling functions: Universal insights into dynamical interaction mechanisms

Coupling functions: Universal insights into dynamical interaction mechanisms

Tomislav Stankovski    Tiago Pereira    Peter V. E. McClintock    Aneta Stefanovska Faculty of Medicine, Ss Cyril and Methodius University, 50 Divizija 6, Skopje 1000, Macedonia Department of Physics, Lancaster University, Lancaster, LA1 4YB, United Kingdom Department of Mathematics, Imperial College London, London SW7 2AZ, United Kingdom Institute of Mathematical and Computer Sciences, University of São Paulo, São Carlos 13566-590, Brazil
July 19, 2019
Abstract

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

pacs:
05.45.Xt, 05.45.-a, 05.45.Tp, 02.50.Tt
Contents

I Introduction

i.1 Coupling functions, their nature and uses

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

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

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

Given these properties, it is hardly surprising that coupling functions have recently attracted considerable attention within the scientific community. They have mediated applications, not only in different subfields of physics, but also beyond physics, predicated by the development of powerful methods enabling the reconstruction of coupling functions from measured data. The reconstruction within these methods is based on a variety of inference techniques, e.g. least squares and kernel smoothing fits Rosenblum and Pikovsky (2001); Kralemann et al. (2013a), dynamical Bayesian inference Stankovski et al. (2012), maximum likelihood (multiple-shooting) methods Tokuda et al. (2007), stochastic modeling Schwabedal and Pikovsky (2010) and the phase resetting Galán et al. (2005); Levnajić and Pikovsky (2011); Timme (2007).

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

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

Figure 1: (color online). Examples of coupling functions used in chemistry, cardiorespiratory physiology and secure communications, to demonstrate their diversity of applications. (a) Coupling functions used for controlling and engineering the interactions of two (left) and four (right) non-identical electrochemical oscillations. (b) Human cardiorespiratory coupling function reconstructed from the phase dynamics the heart and respiration phases. (c) Schematic description of the coupling function encryption protocol. Multiple information signals are encrypted by modulating the parameters of linearly-independent coupling functions between (chaotic) dynamical systems at the transmitter. These applications are discussed in detail in Sec. V. Fig. 1(a) is from Kiss et al. (2007), (b) from Kralemann et al. (2013a) and (c) from Stankovski et al. (2014).

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

i.2 Significance for interacting systems more generally

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

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

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

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

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

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

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

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

i.2.2 Coupling strength and directionality

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

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

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

i.2.3 Coupling functions in general interactions

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

Ii Basic Concept of Coupling Functions

ii.1 Principle meaning

ii.1.1 Generic form of coupled systems

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

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

Figure 4: (color online). The state of synchronization described through phase difference dynamics, versus . Depending on the existence of stable equilibria, the oscillators can be synchronized (a),(c) or unsynchronized (b). Stable points are shown with white circles, while unstable with black circles. Adapted from Kuramoto (1984).

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

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

ii.1.3 Example of coupling function and synchronization

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

(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 stable-unstable 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 saddle-node bifurcation, leading to a stable synchronization.

The coupling functions of real systems are often more complex than the simple sine function presented in Fig. 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 stable-unstable 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 limit-cycle oscillators theoretically. Notably, he considered the phase dynamics of interacting oscillators, where the coupling function was a product of two periodic functions of the form:

(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 epoch-making book Kuramoto (1984). Today this model is the cornerstone for many studies and applications Acebrón et al. (2005); Strogatz (2000), including neuroscience Breakspear et al. (2010); Cumin and Unsworth (2007); Cabral et al. (2014), Josephson-junction arrays Wiesenfeld et al. (1996, 1998); Filatrella et al. (2000), power grids Dorfler and Bullo (2012); Filatrella et al. (2008), glassy states Iatsenko et al. (2014) and laser arrays Vladimirov et al. (2003). The model reduces the full oscillatory dynamics of the oscillators to their phase dynamics, i.e. to so-called phase oscillators, and it studies synchronization phenomena in a large population of such oscillators Kuramoto (1984). By setting out a mean-field description for the interactions, the model provides an exact analytic solution.

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

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

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

(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 limit-cycle, one can apply phase reduction procedures (see Sec. III.2 for further theoretical details) which systematically approximate the high-dimensional dynamical equation of a perturbed limit cycle oscillator with a one-dimensional reduced-phase equation, with just a single phase variable representing the oscillator state Nakao (2015). In uncoupled or weakly-coupled contexts, the phases are associated with zero Lyapunov stability, which means that they are susceptible to tiny perturbations. In this case, one loses the amplitude dynamics, but gains simplicity in terms of the single dimension phase dynamics, which is often sufficient to treat certain effects of the interactions, e.g. phase synchronization. Thus phase connectivity is defined by the connection and influence between such phase systems.

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

(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.

Figure 5: (color online). Schematic illustration of a phase dynamics coupling function. The first oscillator influences the second oscillator unidirectionally, as indicated by the directional diagram on the left of the figure. (a) Amplitude signal during one cycle of period . (b) Coupling function in space. (c) -averaged projection of the coupling function . (d) Amplitude signal of the second driven oscillator , during one cycle of the first oscillator. From Stankovski et al. (2015).

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 influence-sensitivity product function of Winfree model is given by a more complex form spread differently along the two-dimensional 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.

Figure 6: (color online). Two characteristic coupling functions in the phase domain. (a) The coupling function is of sinusoidal form for the phase difference, as used in the Kuramoto model. (b) The coupling function is a product of the influence and sensitivity functions, as used in the Winfree model.

ii.3.2 Amplitude coupling functions

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

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

To present the basics of amplitude coupling functions, we discuss a simple example of two interacting Poincaré limit-cycle oscillators. In the autonomous case, each of them is given by the polar (radial and angular ) coordinates as: and . In this way, a Poincaré oscillator is given by a circular limit-cycle and monotonically growing (isochronous) phase defined by the frequency parameter. In our example, we transform the polar variables to Cartesian (state space) coordinates , , 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 .

Figure 7: (color online). Schematic illustration of an amplitude dynamics coupling function. The first oscillator Eqs. (10) is influencing the second oscillator Eqs. (11) unidirectionally, as indicated by the directional diagram on the left of the figure. (a) Amplitude state signal during one cycle of period . (b) Coupling function in space during one period of each of the oscillations. (c) -averaged projection of the coupling function . (d) Amplitude signal of the second (driven) oscillator , during one cycle of the first oscillator.

Figure 3: (color online). The state of synchronization described through phase difference dynamics, versus . Depending on the existence of stable equilibria, the oscillators can be synchronized (a),(c) or unsynchronized (b). Stable points are shown with white circles, while unstable with black circles. Adapted from Kuramoto (1984).
Figure 10: (color online). Inference of multivariate interactions. True (structural) configurations (left), and the reconstructed phase model (right). Middle: the table shows the corresponding inferred coupling strengths. Note the multivariate triplet link – the arrows from the centres of the diagrams. From Kralemann et al. (2011).

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 pre-estimated most probable models Berger and Pericchi (1996), but without including explicit assessment of the coupling functions. Due to the multi-dimensionality and the lack of a general property in a dynamical system (like for example the periodicity in phase dynamics), there are countless possibilities for generalization of the coupling function. In a sense, this lack of general models is a deficiency in relation to the wider treatment of amplitude coupling functions. There are open questions here and much room for further work on generalising such models, in terms both of theory and methods, taking into account the amplitude properties of subgroups of dynamical systems, including for example the chaotic, oscillatory, or reaction-diffusion nature of the systems.

ii.3.3 Multivariate coupling functions

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

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

Multivariate coupling functions have been studied by inference of small-scale networks where the structural coupling can differ from the inferred effective coupling Kralemann et al. (2011). The authors considered a network of three van der Pol oscillators where, in addition to pairwise couplings, there was also a joint multivariate cross-coupling function, for example of the form 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 cross-frequency coupling in neuroscience Jensen and Colgin (2007). Recently, neural cross-frequency coupling functions were extracted from multivariate networks Stankovski et al. (2015) (see also Sec. V.3). The network interactions between the five brainwave oscillations , , , 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 .

Figure 11: (color online). Multivariate triplet coupling functions between neural oscillations. The phase coupling function shows the influence that and jointly insert on the cortical oscillations. From Stankovski et al. (2015).

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

ii.3.4 Generality of coupling functions

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

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

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)
Table 1: Different examples of coupling functions . These pairwise coupling functions (CFs) are considered in relation to the system: .

In the literature, authors often speak of the commonly-used coupling functions including, but not limited to, those listed in Table 1. Note that reactive and diffusive coupling have functionally the same form, the difference being that the reactive case includes complex amplitudes. This results in a phase difference between the coupling and the dynamics. Also in the literature, a diffusive coupling function 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 limit-cycle shape parameters. By definition, the mechanism involves some kind of function or process leading to a change in the affected system. Its significance is that it may lead to qualitative transitions and induce or reduce physical effects, including synchronization, instability, amplitude death, or oscillation death.

But why is the mechanism important, and how it can be used? The first and foremost use of the coupling function mechanism is to illuminate the nature of the interactions themselves. For example, the coupling function of the Belousov-Zhabotinsky chemical oscillator has been reconstructed Miyazaki and Kinoshita (2006) with the help of a method for the inference of phase dynamics. Fig. 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 spike-like i.e. they resemble delta functions. Consequently, neuronal coupling functions (which are convolution of phase response curves and perturbation functions) then depend only, or mainly, on the phase response curves. So the interaction mechanism is defined by the phase response curves: quite a lot of work has been done in this direction Schultheiss et al. (2011); Tateno and Robinson (2007); Gouwens et al. (2010); Ermentrout (1996); see also Sec. IV.4.2. For example, Tateno and Robinson (2007) and Gouwens et al. (2010) reconstructed experimentally the phase response curves for different types of interneurons in rat cortex, in order to better understand the mechanisms of neural synchronization.

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

After learning the details of the reconstructed coupling function, one can use this knowledge to study or detect the physical effects of the interactions. In this way, the synchronous behavior of the two coupled Belousov-Zhabotinsky reactors can be explained in terms of the coupling function as illustrated by the examples given in Fig. 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.

Figure 12: Coupling function determined from the phase dynamics of two interacting chemical Belousov-Zhabotinsky oscillators. The coupling function is reconstructed in terms of the phase difference . Points obtained from reactors 1 and 2 are plotted with open circles and triangles, respectively. The full curves represent smooth interpolations. From Miyazaki and Kinoshita (2006).

ii.5 Synchronization prediction with coupling functions

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

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

Figure 13: Experimental coupling function from electrochemical oscillators, used for the prediction of synchronization. (a)-(c) Coupling function evaluated in respect of the phase difference shown on the left panel and its odd part shown on the right panel – for the case of (a) smooth oscillator, and (b) and (c) for relaxation oscillator with slightly different parameters. on the plots is equivalent to the notation used in the current review. From Kiss et al. (2005).

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 in-phase synchronization at the phase difference of . The third case of relaxation oscillations Fig. 13(c) has an inverted sinusoidal form , leading to stable anti-phase 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 in-phase () and anti-phase () entrainments are stable, in which case the actual state attained will depend on the initial conditions.

Next, the knowledge obtained from experiments with a single oscillator was applied to predict the onset of synchronization in experiments with 64 globally coupled oscillators. The experiments confirmed that for smooth oscillators the interactions converge to a single cluster, and for relaxational oscillators they converge to a two-cluster synchronized state. Experiments in a parameter region between these states, in which bistability is predicted, are shown in Fig. 14. A small perturbation of the stable one-cluster state (left panel of Fig. 14) yields a stable two-cluster 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.

Figure 14: Mutual entrainment and stable single-cluster (left panel) and two-cluster (right panel) states of a population of 64 globally-coupled electrochemical relaxation oscillators under the same experimental conditions. The two-cluster state was obtained from the one-cluster state by a small perturbation acting as a different initial condition for the population. From Kiss et al. (2005).

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

ii.6 Unifying nomenclature

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

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

Iii Theory

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

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

iii.1 Strong interaction

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

iii.1.1 Two coupled oscillators

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

Diffusion driven oscillations

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

Assume that the system

(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 diffusion-driven oscillators. Further oscillations in originally passive systems have been reported in spatially extended systems Gomez-Marin et al. (2007). In diffusively coupled membranes, collective oscillation in a group of nonoscillatory cells can also occur as a result of spatially inhomogeneous activation factor Ma and Yoshikawa (2009). These ideas of diffusion leading to chemical differentiation have also been observed experimentally and generalized by including heterogeneity in the model Tompkins et al. (2014).

Oscillation death

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

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

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.

Figure 15: (color online) Illustration of the coordinates parallel and transverse to synchronization. In the left panel we also show a trajectory converging to the synchronization subspace implying that . Once the coupled systems reach synchronization, their amplitudes will evolve together in time, but the evolution can be chaotic as illustrated in the right panel. The dynamics along the synchronization subspace is the Lorenz attractor.

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.

  1. 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).

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

  3. 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 Taylor-expand 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 Taylor-expand 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: Comparison between classes of coupling function and the techniques to obtain synchronization. Dashes indicate that typically we are unable to construct the Lyapunov function in such cases.

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

iii.2 Weak regime

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

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

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

iii.2.1 Stable Periodic Orbit and its phase

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

(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 re-parameterize 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.

Figure 16: (color online) Periodic orbits are shown as filled (black) circles and isochrons as (blue) lines. Every point in an isochron has the same value of asymptotic phase. Moreover, the distance between two points in the same isochron tends to zero exponentially fast, as illustrated by the red (dark) and blue (light) points. In the lower figure, we show the effect on the phase dynamics of a small perturbation. The point is initially at phase zero. The perturbation moves the system from its initial point to another isochron, thereby advancing the phase. The periodic orbit and the isochrons are for Eq. 25.

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.

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

iii.2.3 Synchronization with external forcing

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

(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 phase-lock 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 so-called phase response curve or phase resetting curve (PRC). For an oscillator to be able to adjust its rhythm and synchronize, it must respond differently to the perturbations at different phases . 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 so-called phase response curve.

In the theory of weakly coupled oscillators, we use the concept of an infinitesimal PRC (iPRC). It is equivalent to the gradient of the phase , 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, pulse-coupled oscillators are an important class of models. Here the interactions happen in instantaneous pulses of communication. The collective dynamics of such models are of great interest Mirollo and Strogatz (1990). The relationship between pulse-coupled oscillators and the phase reduction presented above has recently been elucidated by Politi and Rosenblum (2015) who showed that the models are equivalent.

iii.2.5 Examples of the phase sensitivity function

Because of the works of Winfree (1967, 1980), Kuramoto (1984), Ermentrout (1996); Stiefel et al. (2008), we now have a good understanding of the phase sensitivity 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.

Figure 17: (color online) Three typical bifurcations appearing as the result of changing a single parameter. As the parameter changes for the SNIPER bifurcation, two fixed points collapse to a saddle node on the circle, and then the system oscillates. In the Hopf bifurcation, a periodic orbit appears after destabilization of the fixed point. In the homoclinic bifurcation, the stable and unstable manifolds of the saddle point join to form a homoclinic orbit, as the parameter changes; with further parameter change, the homoclinic orbit is destroyed and a periodic orbit appears.

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, 111We 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
Table 3: The phase sensitiveness for various models and their bifurcation after Brown et al. (2004). Izhikevich (2000, 2007) obtained the phase sensitivity for relaxation oscillators. They are discontinuous and are not shown here.

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 best-studied coupling function, and it has offered deep insights into collective properties for both globally coupled oscillators Acebrón et al. (2005) and complex networks Arenas et al. (2008); Rodrigues et al. (2016).

First, we consider . 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 mean-field parameters)

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

When the distribution of the natural frequencies is an even, unimodal, and non-increasing function, and the coupling is weak, the incoherent state is neutral Strogatz and Mirollo (1991), but the order parameter 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 bi-modal distributions (keeping the coupling function purely harmonic) has been given by Acebrón et al. (2005).

iii.3.1 Coupling functions leading to multistability

As we have seen, if the oscillators are close to a Hopf bifurcation as is typical, then the corresponding phase sensitivity . This means that the coupling function will have a phase shift

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

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

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

Higher harmonics. In the previous discussion, the coupling function contained one harmonic. . Depending on the underlying bifurcation we must now include further, higher-order, Fourier components,