Synchronization transitions caused by time-varying coupling functions

Interacting dynamical systems are widespread in nature. The influence that one such system exerts on another is described by a coupling function; and the coupling functions extracted from the time-series of interacting dynamical systems are often found to be time-varying. Although much effort has been devoted to the analysis of coupling functions, the influence of time-variability on the associated dynamics remains largely unexplored. Motivated especially by coupling functions in biology, including the cardiorespiratory and neural delta-alpha coupling functions, this paper offers a contribution to the understanding of effects due to time-varying interactions. Through both numerics and mathematically rigorous theoretical consideration, we show that for time-variable coupling functions with time-independent net coupling strength, transitions into and out of phase- synchronization can occur, even though the frozen coupling functions determine phase-synchronization solely by virtue of their net coupling strength. Thus the information about interactions provided by the shape of coupling functions plays a greater role in determining behaviour when these coupling functions are time-variable. This article is part of the theme issue ‘Coupling functions: dynamical interaction mechanisms in the physical, biological and social sciences’.


SW7 2AZ, UK
TS, 0000-0003-2974-2836; AS, 0000-0001-6952-8370 Interacting dynamical systems are widespread in nature. The influence that one such system exerts on another is described by a coupling function; and the coupling functions extracted from the time-series of interacting dynamical systems are often found to be time-varying. Although much effort has been devoted to the analysis of coupling functions, the influence of time-variability on the associated dynamics remains largely unexplored. Motivated especially by coupling functions in biology, including the cardiorespiratory and neural delta-alpha coupling functions, this paper offers a contribution to the understanding of effects due to time-varying interactions. Through both numerics and mathematically rigorous theoretical consideration, we show that for time-variable coupling functions with time-independent net coupling strength, transitions into and out of phasesynchronization can occur, even though the frozen coupling functions determine phasesynchronization solely by virtue of their net coupling strength. Thus the information about interactions 2019 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/ by/4.0/, which permits unrestricted use, provided the original author and source are credited.
In a differential equation or stochastic differential equation describing a system of interacting components, the terms on the right-hand side arising from the interactions between the components are referred to as coupling functions. Coupling functions can be reconstructed from time series recorded from the interacting components, as a result of which one can obtain information about the interactions. For example, from cardiac and respiratory time series one can obtain the phase at which the cardiac beat is most susceptible to respiratory drive, from which one can extract the respiratory-related component of heart rate variability [26]. Another example is that general anaesthesia can lead to important changes in the forms of coupling function between brain waves [11]. Several studies show that the possession of time-varying coupling functions typifies the behaviour of interacting systems in real situations [27][28][29], and that such time-variations can play a major role in the dynamics [30][31][32].
The coupling function can be described by its net coupling strength, and its form [33]. By the net coupling strength, we mean the norm of the coupling function; its form, on the other hand, defines the functional law specifying the interactions, and it thereby introduces a new dimension and perspective [28,29,33]. Thus the net coupling strength quantifies only one aspect of the coupling. Many recent studies of interactions are designed for, and focus exclusively on, the effect of the net coupling strength of interacting systems. This approach is often found in informationtheoretic methods for the detection of directionality and causality of influence between time-series including, for example, methods for Granger causality, transfer entropy, mutual information and symbolic transfer entropy [34][35][36][37]. From this perspective, time-variability of coupling strength was observed in cardiorespiratory interaction in [36].
In this paper, we extend [23] by studying theoretically the effects of time-varying coupling functions that induce a transition to synchronization, while keeping the net coupling strength constant. In particular, we will analyse numerically a model of two coupled oscillators with timevarying coupling functions while maintaining a constant net coupling strength. We will thus obtain information about synchronization epochs and phase slips in terms of the dynamics of the phase difference of the oscillators. Then we will set out theoretical considerations and theorems for quantifying these phenomena. For all of this, we will consider unidirectional coupling, i.e. the master-slave configuration. We will demonstrate that fixing a constant value for the net coupling strength does not necessarily enable one to predict whether or not synchronization transitions will occur, even though for the same set-up in the absence of time-variability, the value of the net coupling strength would have been sufficient to determine synchronization.
The paper is organized as follows. To set the context and to provide some motivation for the study, §2 describes an experiment on the time-variability of biological interactions. Section 3 presents a typical model of two coupled phase oscillators. Some basic concepts regarding coupling functions are presented in §4. Section

Motivation from time-variability of biological interactions
Much of the motivation for studying time-varying coupling functions has come from biology. In particular, numerous situations arise where there are inherent time-variations, not only in the internal parameters and the quantitative characteristics of components of the system but also in the physical laws and functions that define the interactions between these components.
For example, the coupling functions of cardiorespiratory interactions were found to vary in time in [27], where it was shown that not only the coupling strength but also the form of the coupling function varies over time. Similarly, human delta-alpha neural coupling functions in the resting state also vary in time [28]. To illustrate, we now consider an example of time-varying cardiorespiratory and delta-alpha neural coupling functions, calculated using simultaneous recordings from the same subject. More specifically, the cardiorespiratory interactions were analysed from an ECG signal and a respiration-belt signal, while the delta and alpha brainwaves were extracted from an EEG signal, measured on the left forehead (equivalent to an FP1 electrode in the 10-20 EEG standard). These particular relationships were chosen for analysis, as they had previously been found significant when tested against surrogates. During recordings, the subject was resting in a supine position. The data are drawn from an earlier study of the effects of general anaesthesia on physiological oscillations [11]. Here, in order to demonstrate the time-variability properties, we consider data from only one of the earlier subjects.
Coupling functions were extracted from the phase dynamics using the dynamical Bayesian inference method [27,38,39], which involves the use of a sliding time-window to reveal the timedependence. It was assumed that within each time-window, the phases φ 1 (t) and φ 2 (t) of the two interacting oscillators are governed by coupled differential equations of the general form (3.3, see §3 below) with the addition of Gaussian white noise, i.e.
where ω i and φ i are the natural frequency and phase of oscillator i, ξ i (t) is Gaussian white noise, and q i is the coupling function describing the influence of oscillator j on the phase of oscillator i. First, from the recorded time-series the phases of the cardiac φ h , the respiratory φ r , the delta brainwave φ δ and the alpha brainwave φ α oscillations were extracted using the Hilbert protophase-to-phase procedure [40]. The coupling functions were then inferred within each timewindow up to a second order of Fourier expansion on the 2-torus. In this way, we inferred the time-evolving cardiorespiratory coupling function q h (φ r , φ h ) for the influence of the respiration on the heart, and the time-evolving neural cross-frequency coupling function q α (φ δ , φ α ) for the influence of δ brainwaves on α brainwaves. Figure 1 shows the results for cardiorespiratory interactions in (a-c) and for interactions in the neural delta and alpha waves in (d-f ). We show the time-variations in quantitative measures of the coupling functions (figure 1a), namely the coupling strength ε(t), and the similarity index ρ(t) ∈ [−1, 1] between the form of the coupling function and the form of the time-averaged coupling function. The coupling strength ε(t) is the norm (with respect to L 2 ) of the coupling function at time t, as described in §4a; the similarity index ρ(t) [26,29,33,41] is the cosine similarity (again with respect to L 2 ) between the coupling function at time t and the time-averaged coupling function. Note that both ε(t) and ρ(t) are time-varying, but they often vary differently, reinforcing the argument that the strength and form of the coupling represent two separate dimensions of the coupling function, often carrying different information about the interactions.    Comparing the cardiorespiratory and neural interactions shows that the former exhibits the more stable and invariant form, in that it varies less between different time windows, arguably implying greater determinism as described in the time-averaged coupling function. On the other hand, the neural coupling functions vary more and are individually less similar to the timeaveraged coupling function. (These findings are consistent with previous results for the resting state in the multi-subject studies [11]). Common to both of the interactions is that they are time-varying, to a lesser or greater extent, and that the coupling strength and the form of the function often vary over time quite differently from each other. Hence, these characteristics can have correspondingly different effects on the outcome and the possible transitions caused by the interactions-a phenomenon worth exploring further theoretically.

Model
Rhythmic phenomena can be described by a stable periodic dynamics. Consider two weakly coupled oscillators described by the dynamical system where f 1,2 represent the unperturbed dynamics of the first and second oscillators, g 1,2 represent the effects of one oscillator on the other and ε is a small parameter which represents the strength of the perturbations. We assume that f i has an exponentially stable limit cycle γ i of period T i for i = 1, 2; and that the phase of the ith oscillator φ i is defined on the limit cycle γ i in such a way that it grows monotonically with time, satisfying the phase equation for i = 1, 2, where ω i = 2π/T i is the natural frequency of the ith oscillator. One can apply the phase reduction method to reduce the dynamics of the high-dimensional system equation (3.1) to lower-dimensional phase equations [42][43][44][45][46]. Small external perturbations to each oscillator x i , such as its interaction with the other oscillator x j , may force x i off the limit cycle γ i of f i . Therefore, for phase reduction of equation (3.1), one needs to define the phases outside the limit cycles. Using the concept of 'isochrons' [43,47,48] which are the level sets of phases, one can extend the definitions of the phases φ 1 and φ 2 of the oscillators to the whole basins of the corresponding limit cycles γ 1 and γ 2 in such a way that they rotate uniformly according to equation (3.2), not only on the cycle, but also in their corresponding neighbourhoods. With this, the model (3.2) of interacting oscillators can be reduced to a phase model taking the form and where the functions q 1 (φ 1 , φ 2 ) and q 2 (φ 1 , φ 2 ) are 2π -periodic with respect to their arguments. We refer to these functions q 1 and q 2 as the coupling functions for the phase-reduced model (3.3). When considering time-variable interactions, the coupling functions will be time-dependent functions q 1 t (φ 1 , φ 2 ) and q 2 t (φ 1 , φ 2 ). In what follows, we will specifically consider unidirectional coupling, meaning that q 1 = 0. But first, before we carry out our numerical and theoretical analysis of the model, we will in the next section introduce the concepts of net coupling strength and synchronization transitions.

Basic concepts (a) Net coupling strength
Because each coupling function q is 2π -periodic in both arguments, we can identify it on a torus of the form T 2 = S 1 × S 1 , where S 1 ∼ = R/2π Z. Hence, the torus T 2 can be identified by the square [−π , π ] 2 or [0, 2π ] 2 . We define the inner product in L 2 (T 2 , R) [49] by If the coupling function q is smooth, then we can decompose q into a Fourier series on the square [−π , π ] 2 , using Parseval's identity [49] to compute the norm in terms of Fourier coefficients. In applications, the coupling functions can be well approximated by using a finite number of Fourier terms. Each coupling function is typically of the diffusive type where c 1 (t) and c 2 (t) are time-varying parameters. With the second expression for q t in equation (4.1) being a finite sum of Fourier components, Parseval's identity gives This norm provides a quantitative measure of the net coupling strength between the linked systems. For any specified (time-dependent) form of coupling function, this norm represents the scaling parameter of the coupling function.

(b) Intermittent synchronization
A mark of the level of adjustment of rhythmic behaviour due to interaction is whether it is sufficient to cause synchronization of the oscillators [42]. For time-variable coupling functions, the time-variability can cause transitions into or out of synchronization, i.e. it can cause there to be both epochs of synchrony during which the phase difference remains nearly constant and epochs of phase slipping during which the phase difference changes rapidly. An example to illustrate these behaviours is shown in figures 2 and 3 which plot time series of the phase difference between two coupled oscillators. We will sometimes refer to this behaviour as intermittent synchronization. We emphasize that these transitions into and out of synchronization are not the same as the intermittency of apparent synchronization often occurring in autonomous systems with parameters close to or on the boundary of a region of synchronization or chaos [42]. Transitions into and out of synchronization occurring due to time-variability of the coupling functions in equation (3.3) correspond physically to dynamical consequences of the openness of a coupled-oscillator system to time-variable external influences. An example of intermittent synchronization induced by such non-autonomicity of models of open systems is the focus of this present paper. It is also known that for closed systems, transitions into or out of synchronization can occur due to the time-evolution of slow variables within the system representing internal adaptation based on the current or previous state of the oscillators, as exemplified in [50] for a large network with slowly adaptive coupling via time-delays. The intermittency in this latter case is similar to intermittency induced by non-autonomicity, but represents a very different physical cause of the exhibited dynamical behaviour.

Numerics
Our main goal is to establish numerically the effect of time-variation in the coupling functions between two-phase oscillators. Our focus will be on the case where the coupling exists in one  Specifically, c 1 (t) = √ 2α cos(f (t)t) and c 2 (t) = √ 2α sin(f (t)t) as in equation (5.3), where f (t) is the periodic function defined in equation (5.4). In red is shown the phase difference ψ(t) = φ 1 (t) − φ 2 (t) as governed by equation (5.2), and in blue is shown f (t). The parameters ε and k were set to ε = 0.01 rad s −1 , k = 100 rad s −1 , and the net coupling strength was set to α = 1.55/ √ 2 s −1 . The inset shows the transition to synchronization. The dynamics of the phase difference is shown to alternate between synchrony states and phase slips (indicated by bold arrows in the plot of ψ(t)), due to the time-variability of the coupling function q t in equation (4.1) via the parameters c 1 (t) and c 2 (t) while the net coupling strength remains constant. direction only, sometimes known as unidirectional coupling, or the master-slave configuration. We choose this configuration because it provides the clearest case where time-varying coupling functions can lead to synchrony and phase slips, while the net coupling strength remains constant.
Consider the master-slave configuration where the coupling function q t (φ 1 , φ 2 ) is equal to the expression in equation (4.1) and ω 1 , ω 2 are the natural frequencies of the oscillators. The presence of the coupling term q t (φ 1 , φ 2 ) could cause the fundamental frequency of the driven oscillator (whose phase is represented by φ 2 (t)) to become different from its natural frequency ω 2 , and to become time-dependent as q t varies over time.
From the previous section, the net coupling strength of the master-slave configuration q t in equation (5.1) can be defined as where c 1 (t) and c 2 (t) are the time-varying coupling parameters of the coupling functions [  y (t) |ω 1 − ω 2 | then the oscillators will synchronize [23]. In the present study, unlike the autonomous case, the oscillators are not guaranteed to be synchronized all of the time.
To analyse whether the two oscillators synchronize or not, we consider the phase difference between them. From equation (5.1) with q t as in equation (4.1), the phase difference ψ(t) obeys the equation where Ω = ω 1 − ω 2 is the natural frequency difference (sometimes called the frequency mismatch, or detuning) between the oscillators. For numerical investigation of synchrony in our coupled oscillator model, we simulate the differential equation (5.2), taking the time-dependent coupling parameters c 1 (t) and c 2 (t) as where f (t) is a T-periodic function defined by with k > ε ≥ 0, and with the values of ε, (T/2 − T 1 )/T and (T − T 2 )/T being small. The expression for the function in equation (5.4) has been chosen to exhibit the existence of synchrony epochs and phase slips in the dynamics of the phase difference. Note that, from equation (5.3), we obtain showing that the net coupling strength is constant for all time.
In all our simulations, the length of the simulated time series is set to 5000 s and the sampling step is set as h = 0.001 s, and we use the following parameter values: Ω = 1.08 rad s −1 , T = 1000 s, T 1 = 490 s, T 2 = 990 s. We will consider different possible values for the remaining parameters k, ε and α. We simulate equation (5.2) via its simplified form equation (6.1) with ϕ(t) := f (t)t, taking the initial value of the phase difference as ψ(0) = 0. Figure 2 presents a time series of the phase difference ψ(t) (red) and the function f (t) (blue), and the zoomed figure (top-left inset) shows the transition to synchrony. Further numerical simulations were carried out using different parameter values, for which results are shown in figure 3. For both figures, the net coupling strength α is constant over time, with √ 2α being greater than the magnitude of the natural frequency difference Ω.
In both figures 2 and 3a-d, we see intermittency of synchronization, similar to that observed experimentally and numerically in [7,20,24,27,51]: there is an alternation between synchronized epochs (plateaux) and phase slips (rapid increases) in ψ(t). While ϕ(t) := f (t)t is slowly varyingi.e. on the intervals [nT, nT + T 1 ] where f (t) has small magnitude-we have phase synchrony. But while ϕ(t) has rapid angular velocity-i.e. on the intervals [(n + 1 2 )T, nT + T 2 ] where f (t) is large-we have unbounded slipping of the phase difference.
Given that the net coupling strength between the systems was invariant, it is evident that the continuing alternation between synchronization epochs and phase slips was just due to timevariation in the coupling function. This shows that the net coupling strength does not in itself give us enough information to characterize the interactions of the oscillators. Based on the numerics implemented in this section, we can generalize the choices of the coupling parameters c 1 (t) and c 2 (t) in order to analyse the effects of the time-varying coupling functions.

Explanation and generalization of numerical findings
Consider equation (5.3) with a general function ϕ(t) in place of f (t)t. Based on its behaviour, we determine the dynamics of the phase difference ψ(t) of the oscillators. Whether they exhibit synchronous and/or asynchronous states hinges on three considerations: (i) If ϕ(t) is a slowly varying function, i.e. |φ| is small, then the phase difference ψ(t) corresponds to a synchronous state for a finite time T 1 : for arbitrarily large T 1 we can take |φ| sufficiently small that the solution will remain in an arbitrarily small ball. See theorem 6.1 and corollary 6.2 for the proof. (ii) If the phase shift ϕ(t) is fast-winding round the circle (i.e. its unwrapped phase is rapidly growing) then, using an averaging argument, we can show that the dynamics induces phase slips. See theorem 6.3. The averaging argument essentially means that due to the fast timescale of ϕ, the oscillators φ 1 and φ 2 feel no interactions and they rotate independently at their own natural angular velocities. (iii) If the function ϕ(t) has both the behaviours stated in (i) and (ii), then the coupled oscillators will synchronize when ϕ proceeds on the slow timescale, but they will exhibit phase drifts when ϕ proceeds on the fast timescale. Hence, the dynamics of the phase difference in equation (6.1) subject to a function ϕ(t) with both slow variation and fast winding induces both synchronous and asynchronous states.
We now provide rigorous theorems to justify these three considerations. Before doing so, however, we briefly discuss the setup. The phase difference ψ(t) = φ 1 (t) − φ 2 (t) obeys equation (5.2). We take the parameters c 1 (t) and c 2 (t) in equation (5.2) to be of the form c 1 (t) = √ 2α cos(ϕ(t)) and c 2 (t) = √ 2α sin(ϕ(t)), for some constant α > 0 and C 1 function ϕ : R → S 1 ∼ = R/2π Z. The constant α corresponds to the net coupling strength. Equation (5.2) can be written as Note that equation (6.1) is a non-autonomous system, and we can consider the effect of both slow variation and fast variation of ϕ(t). Introducing a new variable η(t) = ψ(t) + ϕ(t), the dynamics of equation (6.1) is then equivalent to the dynamics of the equation In the case that √ 2α > |Ω|, let η * and η * * be respectively the stable and unstable fixed points of the equation 3) The following theorems address the effect of the time-varying coupling functions, while the net coupling strength is invariant. First, in Theorem 6.1 and its subsequent corollary, we consider the synchronization exhibited for slowly varying ϕ.
Combining these gives the result.
Proof. Take ε 0 and T 0 as in Theorem 6.1, and then for T 1 > T 0 take ε = min(ε 0 ,ε 2T 1 ). So if |φ| ≤ ε 2ε for all t ∈ [0, T 1 ], and so we have the desired result. Now in Theorem 6.3, we consider the unbounded phase slips exhibited for fast-winding ϕ. The following theorem may be regarded as a kind of non-autonomous averaging principle. For a detailed exposition of averaging principles, see [52]. Proof. As in the proof of Theorem 6.1, let F(η) := Ω − √ 2α sin(η). First take arbitrary K 1 > 0 and K 2 > √ 2π, and suppose thatφ(s) > max(K 1 , and therefore Now fix any s ∈ (t i , t i+1 ). We have that s − t i < 2π/K 1 and so Also, by Taylor's theorem, we have that for some ξ 1 (s) ∈ (t i , s). But by the mean value theorem we have that for some ξ 2 (s) ∈ (t i , ξ 1 (s)), Combining equations (6.4), (6.5) and (6.6), we have that Hence, for each i ≤ N, Now for any s ∈ [0, t], taking the largest i with t i ≤ s, we have and so .
Combining theorems 6.1 and 6.3, we have that if ϕ has alternating epochs of slow variation and rapid oscillation, this can lead to the system (5.1) having alternating epochs of synchrony and asynchrony.

The results in context
There are different ways in which time-variability can enter a system, and so other forms of nonautonomous driving that give rise to intermittent synchronization have also been considered. Of particular interest has been the case where time-variability enters through modulation of the natural frequency of the driving oscillator [22,24,53]. As shown in [24,53], for fixed-frequency driving we have synchronization either all of the time or none of the time, but when the driving frequency is allowed to vary then we can have intermittent synchronization. Through the finitetime dynamical considerations of [24] to account for a free shape of variation over time, it was shown that this intermittent synchronization results overall in stability of the driven oscillations, as indicated by negativity of Lyapunov exponents; and in this manner, greater time-variability of the driving frequency increases the region of stability in parameter space. Similar results were also observed numerically for higher-dimensional oscillatory systems.
The theoretical and numerical considerations of [24] for a unidirectionally coupled pair of phase oscillators have been extended to phase-oscillator networks in [54]. The considerations of this present paper can be generalized to networks such as those considered in [55]: network structure can also have an impact on the dynamics exhibited, such as synchronization [30,56]. Network synchronization induced by sufficiently fast non-autonomous driving has been considered in [57], where the non-autonomous driving consists of rapid switching between the existence and non-existence of links between given nodes, while the dynamical interaction along existing nodes takes a specific time-independent form.
As intermittent synchronization in [24,54] led to stability, we likewise expect that the intermittent synchronization obtained in our present study leads to stability of the driven oscillations, as would be indicated by a stability analysis of the non-autonomous one-dimensional differential equation where φ(t) is the phase of the driven oscillator and φ 1(0) + ω 1 t is the phase of the driving oscillator.

Conclusion
Interacting dynamical systems can have time-varying coupling functions, where the net coupling strength depends in a number of different ways on the time, sometimes resulting in synchronization transitions. The experiments on the cardiorespiratory and neural delta-alpha coupling functions whose results are shown in figure 1, for example, illustrate the existence of time-varying functional relationships that can cause synchronization transitions. A model of two coupled oscillators with time-evolving coupling functions has been shown to exhibit transitions to/from synchronization even when the net coupling strength remains constant. The analysis was carried out in terms of the phase difference between the oscillators. The corresponding numerical simulations show that, in the case of time-varying coupling functions, one can have sequential epochs of synchrony and asynchrony while the net coupling strength remains unchanged. Thus, by itself, the net coupling strength does not provide enough information to describe the dynamics of the interacting systems. To generalize the results, based on the model considered, we discussed three main ideas arising from the periodic function f (t). The first of these was that, when ϕ(t) := f (t)t varies slowly with time, the dynamics of the two coupled oscillators induces synchrony over the slow timescale. The second was that when, by contrast, ϕ(t) has rapid angular velocity, the oscillators do not synchronize. The third observation was the combined effect: the oscillators can exhibit sudden changes between synchrony and drifting phase difference occurring at transitions between slow variation and fast winding of ϕ(t). So we have transitions in exhibited behaviour due to time-variability of the coupling functions despite constant net coupling strength. This confirms that, in the time-variable setting, the net coupling strength does not give sufficient information about the interaction of the oscillators to predict their behaviour.
Note that this insufficiency of the net coupling strength as a criterion carries implications for information-theoretic methods that assess the statistical mutual dependences of signals from interacting systems [10,[34][35][36]58,59]. The latter are statistical measures that can determine a causal relationship and the predominant direction of influence, thus measuring a directed functional connectivity. In this way, however, they usually reveal only the net coupling strength and direction, but are unable to detect variations in sub-coupling components like those used in the present study.
Data accessibility. Data used for analysis of coupling functions shown in figure 1 were collected during the project BRACCIA https://cordis.europa.eu/project/rcn/74753/factsheet/en and are available upon request to A.S. Codes used to generate data in figures 2 and 3 are also available.