Coupling functions: Universal insights into dynamical interaction mechanisms

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. A coherent and comprehensive review is presented 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.

Figure 1

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) nonidentical electrochemical oscillations. From [154]. (b) Human cardiorespiratory coupling function $Q_e$ reconstructed from the phase dynamics of the heart ${\phi }_e$ and respiration ${\phi }_r$ phases. From [167]. (c) Schematic description of the coupling function encryption protocol. From [285]. 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. 5.

Figure 2

The state of synchronization described through phase difference dynamics $\dot{\psi }$ vs $\psi$. Depending on the existence of stable equilibria, the oscillators can be (a), (c) synchronized or (b) unsynchronized. Stable points are shown as white circles, while unstable are black circles. Adapted from [172].

Figure 3

Schematic illustration of a phase dynamics coupling function. The first oscillator $x_1$ influences the second oscillator $x_2$ unidirectionally, as indicated by the directional diagram on the left. (a) Amplitude signal $x_1(t)$ during one cycle of period $T_1$. (b) Coupling function $q_2({\varphi }_1,{\varphi }_2)$ in $\{{\varphi }_1,{\varphi }_2\}$ space. (c) ${\varphi }_2$-averaged projection of the coupling function $q_2({\varphi }_1,{\varphi }_2)$. (d) Amplitude signal of the second driven oscillator $x_2(t)$, during one cycle of the first oscillator. From [287].

Figure 4

Two characteristic coupling functions in the phase domain. (a) The coupling function $q({\varphi }_1,{\varphi }_2)$ is of sinusoidal form for the phase difference, as used in the Kuramoto model. (b) The coupling function $q({\varphi }_1,{\varphi }_2)$ is a product of the influence and sensitivity functions, as used in the Winfree model.

Figure 5

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. (a) Amplitude state signal $x_1(t)$ during one cycle of period $T_1$. (b) Coupling function $q_2(x_1,x_2)$ in $\{x_1,x_2\}$ space during one period of each of the oscillations. (c) $x_2$-averaged projection of the coupling function $q_2(x_1,x_2)$. (d) Amplitude signal of the second (driven) oscillator $x_2(t)$, during one cycle of the first oscillator.

Figure 6

Inference of multivariate interactions. (Left) True (structural) configurations, and the (right) reconstructed phase model. Middle: The table shows the corresponding inferred coupling strengths. Note the multivariate triplet link, the arrows from the centers of the diagrams. From [168].

Figure 7

Multivariate triplet coupling functions between neural oscillations. The phase coupling function $q_{\gamma }({\varphi }_{\theta },{\varphi }_{\alpha })$ shows the influence that $\theta$ and $\alpha$ jointly insert on the $\gamma$ cortical oscillations. From [287].

Figure 8

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 $\psi ={\varphi }_2-{\varphi }_1$. Points obtained from reactors 1 and 2 are plotted with open circles and triangles, respectively. The full curves represent smooth interpolations. From [206].

Figure 9

Experimental coupling function from electrochemical oscillators used for the prediction of synchronization. (a)–(c) Coupling function $q(\psi )$ evaluated with respect to the phase difference $\psi ={\varphi }_2-{\varphi }_1$ shown in the left panel and its odd part $q_{-}(\mathrm{\Delta }\varphi )$ shown in the right panel—for the case of (a) a smooth oscillator, and (b), (c) for relaxation oscillators with slightly different parameters. $H(\mathrm{\Delta }\varphi )$ on the plots is equivalent to the $q(\psi )$ notation used in the current review. From [156].

Figure 10

Mutual entrainment and stable (left panel) single-cluster and (right panel) two-cluster 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 [156].

Figure 11

Illustration of the coordinates parallel $y$ and transverse $z$ to synchronization. In the left panel we also show a trajectory converging to the synchronization subspace implying that $z\to 0$. 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.

Figure 12

Periodic orbits are shown as solid (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 $\mathrm{\Delta }x$ moves the system from its initial point to another isochron, thereby advancing the phase. The periodic orbit $\gamma$ and the isochrons are for Eq. (25).

Figure 13

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.

Figure 14

Schematic illustration of the procedure for the inference of coupling functions. From left to right: measurement data $\mathcal{M}$, preestimation procedure where the phase or amplitude $\tilde{\mathcal{M}}$ is estimated from those data, the inference of a dynamical model from the $\tilde{\mathcal{M}}$ data, and the coupling function emerging as the end result of the procedure.

Figure 15

The reconstructed functions of the phase interactions. (a) The function ${\mathcal{F}}_1({\varphi }_1,{\varphi }_2)$ for the influence of the second on the first oscillator, and (b) the function ${\mathcal{F}}_2({\varphi }_1,{\varphi }_2)$ for the influence of the first on the second oscillator. From [252].

Figure 16

Application of dynamical Bayesian inference to cardiorespiratory interactions when the (paced) respiration is time varying. (a) The inferred time-varying respiration frequency. (b) The coupling directionality between the heart and respiration (denoted as $h$ and $r$, respectively). (c), (d), and (e) The cardiorespiratory coupling function evaluated for the three time windows whose positions are indicated by the gray arrows. From [284].

Figure 17

(a) The estimated natural frequencies (vertical axis) of 32 electrochemical oscillators vs their measured natural frequencies (horizontal axis). (b) The coupling function estimated by the multiple-shooting method (dotted line), compared with that estimated by application of the perturbation to a single isolated electrochemical oscillator (solid curves). From [313].

Figure 18

Inference of interacting dynamics by application of the random phase resetting method. (a) Four-node network of interacting phase oscillators. (b) Reconstruction of the four-node network. Circles are the actual parameter values and crosses are the inferred values; left $a_{ij}^{(1)}$, right $b_{ij}^{(1)}$, for each pair $i\to j$. (c) Coupling function with respect to the phase difference ${\psi }_{42}={\varphi }_2-{\varphi }_4$ from the reconstructed parameters $a_{42}^{(1)}$ and $b_{42}^{(1)}$. From [184].

Figure 19

Coupling functions for noise-induced oscillations in the FitzHugh-Nagumo model with $b=0.1$ and two different values of the noise intensity $D$: for (left panel) $D=0.08$, the mean frequency is $\omega \approx 0.62$; and for (right panel) $D=0.11$, $\omega \approx 0.95$. From [270].

Figure 20

Box plots illustrating the similarity of cardiorespiratory coupling functions. (a) The correlation coefficient $\rho$ and (b) the difference measure $\eta$, for all available pairs of functions (high similarity corresponds to large $\rho$ and small $\eta$). ES: the similarity between the respiration-ECG coupling functions of the same subject, obtained from two trials. EG: the same relation similarity between different subjects in the group demonstrates low interpersonal variability. PS and PG: intra- and interpersonal similarities, respectively, for the respiration-arterial pulse coupling function. EPS and EPG: intra- and interpersonal similarities, respectively, between between the two types of coupling functions. From [167].

Figure 21

Experimental estimation of neuronal phase response curves (PRCs). (a) Raw estimation of the PRC (dots) and smoothing over a $2\pi /3$ interval (gray line) compared with the estimated PRC (black line) from the approach in [104]. Both curves match, which indicates that the raw data are consistent with a phase model. (b) The same as (a) but after shuffling the raw data. The PRC is roughly flat and yields inconsistent results with the smoothing, implying that the shuffled data cannot be described by a phase model. From [104].

Figure 22

Phase response curve and effective forcing for the cardiorespiratory interactions. (a) Individual PRCs $Z$ and (b) effective forcing $I$ for all ECG-based coupling functions (gray curves). In both main panels the thick (blue) lines show the average over all the individual (gray) curves. The thick (red) lines are obtained by decomposition of the averaged coupling function. The small panel on the top in (a) shows for comparison the average ECG cycle as a function of its phase. The small panel on the top in (b) shows the average respiratory cycle as a function of its phase, with the epochs of inspiration and expiration marked (approximately). From [167].

Figure 23

Engineering a system of four nonidentical oscillators using a specific coupling function to generate sequential cluster patterns. (a) The target [solid line, $H(\mathrm{\Delta }\varphi )=\sin (\mathrm{\Delta }\varphi 1.32)0.25\sin (2\mathrm{\Delta }\varphi )$] and optimized coupling function with feedback (dashed line). (b) Theoretical and experimentally observed heteroclinic orbits and their associated unstable cluster states. (c) Time series of the order parameter [$R1=\sum _{j=1}^N\exp (i{\varphi }_j)$] along with some cluster configurations. (d), (e) Trajectories in state space during slow switching. The black lines represent calculated heteroclinic connections between cluster states (fixed points). The (red) surface in (e) is the set of trajectories traced out by a heterogeneous phase model. $H(\mathrm{\Delta }\varphi )$ on the plots is equivalent to the $q(\psi )$ notation used in the current review. From [154].

Figure 24

The mechanisms of synchronization for Belousov-Zhabotinsky chemical oscillations as determined by their coupling functions. The solid curve shows $Q(\psi )=q(\psi )-q(-\psi )$ estimated from the coupling functions $q(\psi )$. $q(\psi )$ and $q(-\psi )$ are presented as dashed and dot-dashed lines, respectively. Stable and unstable solutions of Eq. (38) are shown as solid and open circles, respectively. $\mathrm{\Delta }\omega$ on the plots is equivalent to the ${\mathrm{\Delta }}_{\omega }$ notation used in this review. From [206].

Figure 25

Experiments on a three-cluster state close to a Hopf bifurcation with negative global coupling of 64 electrochemical oscillators. (a) Current time series and the three-cluster configuration. Solid, dashed, and dotted curves represent the currents from the three clusters. (b) Cluster configuration. White, black, and gray circles represent the three clusters. (c) Response function and wave form (inset) of the electrode potential from a current of single oscillator. (d) Phase coupling function. $\mathrm{\Gamma }(\mathrm{\Delta }\varphi )$ on (d) is equivalent to the $q(\psi )$ notation used in this review. From [163].

Figure 26

Cardiorespiratory coupling functions for the study of human aging. Typical time-averaged coupling functions for (a), (c) a young subject aged 21 yr and (b), (d) an old subject aged 71 yr. (a), (b) are from the cardiac, while (c), (d) are from the respiration phase dynamics. From [132].

Figure 27

Coupling functions for the human cardiorespiratory system. The reconstructed functions specify the dependence of the instantaneous cardiac frequency, measured in radians per second, on the cardiac and respiratory phases. The functions $Q_p({\varphi }_r,{\varphi }_p)$ are computed from the arterial pulse and respiration. Results from the subject who had the lowest levels of determinism and similarity to the coupling functions obtained from ECG phases $Q_e$ are shown in (a), and those for the subject with the highest determinism and similarity in (b). (c) The averaged coupling function over all measurements for all subjects. From [167].

Figure 28

Neural cross-frequency coupling functions between $\delta$ and $\alpha$ oscillations in general anesthesia. (a)–(c) The average coupling functions from all subjects within the group. Note that, for realistic comparison, the vertical scale of the coupling amplitude is the same in each case. Here “Awake” refers to the state when the subject is conscious and resting, and “Propofol” and “Sevoflurane” to states when the subject is anesthetized with propofol or sevoflurane, respectively. From [286].

Figure 29

Examples of neural cross-frequency coupling functions. (a) Spatial distribution of the $\delta -\alpha$ coupling functions over the head, based on the different probe locations. (b) Average coupling function along all the probes for the $\delta -\alpha$ coupling relation. Each $\delta -\alpha$ coupling function $q_{\alpha }({\varphi }_{\delta },{\varphi }_{\alpha })$ is evaluated from the $\alpha$ dynamics and depends on the bivariate $({\varphi }_{\delta },{\varphi }_{\alpha })$ phases. From [287].

Figure 30

Coupling functions of the change in democracy $q_1(D,G)$, from the interaction between GDP and democracy. (a) The two term model, and (b) the five term model. The black lines are the solutions $\dot{D}=0$. The strengths of the coupling functions are encoded by the color bars shown on the side of each figure. From [244].

Figure 31

Coupling functions for the two coupled mechanical metronomes. (Top) Experimental apparatus with the two metronomes, placed on a rigid support. (a), (b) The coupling functions in each direction from the case of coupling with one rubber band, (c), (d) with two rubber bands, and (e), (f) coupling functions for the “uncoupled” without any rubber bands. The vertical scales are the same so that one can clearly see the reduction of the coupling function in the uncoupled case. From [166].

Figure 32

Schematic diagram showing the communication protocol based on coupling functions. Messages $s_1,\dots ,s_n$ are encrypted by modulation of the coupling functions connecting two dynamical systems at the transmitter. Only two signals are transmitted through the public domain. The receiver consists of two systems of the same kind with the same coupling functions (forming the private key) and uses dynamical Bayesian inference to reconstruct $s_1,\dots ,s_n$. From [285].

Figure 33

Transmission of ten pseudorandom binary signals encrypted in different coupling functions. The high values (binary “1”) at the transmitter are indicated by gray shading. The received signals, after decrypting, are shown by thick (red) lines, each of which (a)–(j) represents one information signal $s_i(t)$. The particular amplitude coupling functions that were used for encrypting each signal are indicated on the ordinate axis. The bit words are indicated by $m_1-m_4$ on the top of the figure. From [285].