## Abstract

A mild, nonlinear, time-delayed feedback signal was applied to two heterogeneous oscillators in order to synchronize their frequencies with an arbitrary and controllable phase difference. The feedback was designed using phase models constructed from experimental measurements of the intrinsic dynamical properties of the oscillators. The feedback signal produced an interaction function that corresponds to the desired collective behaviour. The synchronized phase difference between the elements can be tuned to any value on the interval 0 and 2*π* by shifting the phase of the interaction function using the feedback delay. Numerical simulations were conducted and experiments carried out with electrochemical oscillators.

## 1. Introduction

The collective dynamical behaviour of a rhythmic population is dependent on the interactions between individual elements. Coupling among rhythmic elements can lead to synchronization, where a number of rhythmic elements organize into a single group with a uniform phase and frequency. Such synchronization can lead to coherent light emission from laser systems (Petrov *et al.* 1997; Oliva & Strogatz 2001), neuronal clustering in the suprachiasmatic nucleus (Iglesia *et al.* 2000; Yamaguchi *et al.* 2003), oscillations in chemical systems (Epstein & Pojman 1998), epileptic events (Traub & Wong 1982) and Parkinsonian tremors (Tass 1999).

By manipulating the interactions between individual elements, it is possible to steer the collective behaviour of the system to a desired state. A number of methods have been developed to control the collective behaviour of rhythmic systems (Battogtokh & Mikhailov 1996; Popovych *et al.* 2006). Proportional–integral–derivative controller-based feedback has been used to create phase-locked states (Di Donato *et al.* 2007), feedback based on mutual information has been shown effective in controlling the characteristic time scales of oscillators (Belykh *et al.* 2005) and tunable heterogeneities have been used to steer the behaviour of excitable media (Mikhailov & Showalter 2006).

We have developed a general engineering methodology based on phase models that allows precise control over the collective dynamical behaviour of a rhythmic system using global nonlinear time-delayed feedback. A key advantage of this method is that it can be tailored to the unique dynamical behaviour of a physical system. This is accomplished by using experimental measurements to calibrate the phase model to the target system, eliminating the need for physico-chemical models. Mild feedback is used to gently steer the collective behaviour of the target system towards the desired state, ensuring that the rhythmic behaviour of the individual elements of the system is preserved. This method has been used to generate dynamical behaviours such as phase synchronization, phase desynchronization, clustering and itinerant clustering within a population of rhythmic electrochemical elements (Kiss *et al.* 2007; Kori *et al.* 2008).

In this paper, we use our methodology to address the general question of how to tune the stationary phase difference between two dissimilar oscillators without prior knowledge of their underlying dynamic behaviour. This two-oscillator system represents a canonical model for a large class of rhythmic biological systems (Winfree 1980; Iglesia *et al.* 2000) and incorporates many of the non-idealities that are present in typical experimental systems, yet has a simple analytical solution. Heterogeneities are common in experimental and natural rhythmic systems where the intrinsic oscillator frequency is not under experimental control. We demonstrate how to choose and obtain an interaction function that produces an arbitrary phase difference between two rhythmic elements, and we demonstrate the application of the method using an experimental chemical system.

## 2. General methodology

The phase behaviour of a population of oscillating elements can be approximated in the limit of weak global interactions as (Kuramoto 1984)
2.1
where *ϕ*_{i} and *ω*_{i} are the phase and inherent frequency of element *i*, *K* is the interaction strength, *N* is the number of elements in the population and *H*(Δ*ϕ*) is an interaction function that describes the effect one element has on all other elements in the system. The interaction function can be derived from measurable quantities of the experimental system as
2.2
where *Z*(*ϕ*) is the response function, *x*(*ϕ*) is the waveform and *h*(*x*) is the coupling function between the elements. The response function, *Z*(*ϕ*), represents the sensitivity of a single element to perturbations along its limit cycle. This function is an intrinsic property of the oscillator and can be experimentally measured (Galan *et al.* 2005; Kiss *et al.* 2005; Tsubo *et al.* 2007). The coupling function, *h*(*x*), represents the form of the physical connection between the elements in the system. For the simulations and experiments described in this work, the interactions among elements are produced by global nonlinear time-delayed feedback of the form
2.3
where *t* is the time, *S* is the polynomial feedback order, *k*_{n} is the *n*th order polynomial feedback coefficient, *τ*_{n} is the *n*th order polynomial feedback delay and Δ*τ* is a feedback delay that is common to all polynomial terms. By manipulating the shape of the interaction function, it is possible to control the collective behaviour of the system. In general, this can be accomplished using a three-step process: determine a target interaction function that will produce the desired global behaviour, measure the response function of the rhythmic elements to be used and numerically optimize the shape of the coupling function to achieve the target interaction function.

## 3. Tuning the phase difference between two rhythmic elements

### (a) Determining the required interaction function

Applying the Kuramoto phase approximation (equation (2.1)) to a system of two non-identical elements yields equations for the two phase variables {*ϕ*_{1},*ϕ*_{2}} which when subtracted give a relationship for the phase difference, Δ*ϕ*=*ϕ*_{1}−*ϕ*_{2}, between the elements
3.1
where Δ*ω* is the difference between the natural frequencies of the two oscillators. Equation (3.1) can be simplified to
3.2
where *H*^{−} is the odd part of the interaction function, defined as
3.3
From equation (3.2), it is seen that stationary solutions will occur at phase differences that satisfy
3.4
For a homogeneous system (identical elements, Δ*ω*=0), the stationary solutions correspond to the roots of *H*^{−}(Δ*ϕ*). The behaviour of a slightly heterogeneous system (Δ*ω*≠0) can be brought close to that of a homogeneous system by increasing the interaction strength between the elements. In the limit as , the fixed point solutions approach the roots of *H*^{−}(Δ*ϕ*). However, caution must be used with this approach since an excessively large interaction strength will disrupt the waveform and invalidate the phase approximation.

For a given interaction function *H*(Δ*ϕ*), the stationary phase difference, Δ*ϕ**, is obtained from the roots of its odd part, *H*^{−}(Δ*ϕ**)
3.5
By shifting the phase of the interaction function by an amount Δ*τ*, Δ*ϕ** can be tuned to take on values between 0 and 1 rad/2*π*. (As shall be seen below, the phase shift, Δ*τ*, of *H*(Δ*ϕ*) is equal to the common feedback delay in the experiments to be discussed.) The stationary phase differences between the two oscillating elements are then given by the solutions of
3.6
By definition, *H*(Δ*ϕ*) is a 2*π* periodic function; therefore equation (3.6) must always have roots at 0 and 0.5 rad/2*π*. These roots represent the trivial in-phase and anti-phase synchronization states, respectively. The non-trivial roots of equation (3.6) are governed by the higher order terms of *H*(Δ*ϕ*). An interaction function with second-order harmonics is sufficient to obtain a non-trivial, stationary phase difference between the two oscillators. A convenient form for this purpose is the function
3.7
Figure 1*a* illustrates the interaction function defined in equation (3.7) setting *R*=0.1 and Δ*τ*=0, while figure 1*b* illustrates its associated synchronized states as a function of Δ*τ*. The non-trivial bifurcation branches illustrated in figure 1*b* are relatively narrow, indicating that small changes in the parameter Δ*τ* will cause large changes in the synchronized state. To reduce this sensitivity, the parameter *R* was selected to maximize the width of the non-trivial bifurcation branches. To simplify the analysis, the parameters of equation (3.7) were constrained such that only supercritical bifurcations occur.

To determine the optimum value of *R*, a linear stability analysis of equation (3.2) was performed, indicating that stable stationary states occur when d*H*^{−}(Δ*ϕ*)/d Δ*ϕ*>0. Therefore, the phase-locked solution, Δ*ϕ**, is stable if
3.8
From this, it can be seen that the in-phase synchronization state (Δ*ϕ**=0) is stable if
3.9
while the anti-phase state is stable if
3.10
Equations (3.9) and (3.10) indicate that the trivial bifurcation branches correspond to the region where . Since the trivial and non-trivial bifurcation branches are mutually exclusive, the non-trivial branch must correspond to the region where . Therefore, a wide non-trivial bifurcation branch will occur when *H*(Δ*ϕ*) has a large domain of negative slope.

The value for the parameter *R* of equation (3.7) was selected to maximize the width of the non-trivial branch of the bifurcation diagram, that is where . The value of *R* was constrained to the interval 0≤*R*≤0.5 owing to the presence of subcritical bifurcations outside this interval. Numerical analysis determined that the width of the region where monotonically increased with the value of *R* in this interval. As a result, the widest non-trivial bifurcation branch occurred when *R*=0.5. The target interaction function can be seen in figure 1*c* for Δ*τ*=0. With these parameters, equation (3.7) can be analytically solved for *H*^{−}(Δ*ϕ**)=0, yielding Δ*ϕ**=0,*π*, arccos. The stable and unstable stationary states are shown as a function of Δ*τ* in figure 1*d*.

### (b) Numerical simulation of phase difference tuning

Before conducting experiments on the experimental system, the methodology was verified with numerical simulations using the Brusselator oscillator, a simple two-variable ODE system that exhibits self-sustained oscillations (Glansdorff & Prigogine 1971). The phase difference between two interacting Brusselator oscillators is tuned to values between 0 and 2*π*. The dynamical equations for an interacting homogeneous Brusselator population under global feedback are
3.11
where *h*(*x*) is the feedback function defined in equation (2.3) that is constructed from and applied to the variables *x*_{i}. The variables *x*_{i} and *y*_{i} are transformed such that the fixed point is shifted to (*x*,*y*)=(0,0). For a single uncoupled oscillator, a Hopf bifurcation occurs at *B*=*B*_{c}≡1+*A*^{2}. The parameters of equation (3.11) were chosen to be *A*=1.0 (so that *B*_{c}=2.0) and *B*=2.3. The waveform *x*(*ϕ*) and the response functions *Z*(*ϕ*) along the *x*-direction are displayed in figure 2. The response function was calculated using the XPPAUT program (Ermentrout 2002).

A feedback parameter set {*k*_{n},*τ*_{n}} is chosen to create a tunable non-trivial solution such that an arbitrary stationary phase difference Δ*ϕ** can be achieved. The target interaction function was selected to be equation (3.7), where *R*=0.5. The phase shift of the interaction function can be directly controlled using the common feedback delay parameter Δ*τ*. For the purpose of determining the feedback parameters, the common feedback delay can be set to zero. Second-order feedback is used in order to produce the required first and second harmonics in the target interaction function. The constant term in equation (2.3) is arbitrary and can be safely neglected (*k*_{0}=0). Although both the linear and quadratic terms of the feedback function contain weak third (and higher) order harmonics, their effect on the interaction function is negligible since the third (and higher) order harmonics of the response function are small. The feedback parameters are found by numerically solving for *h*(*x*) (Kori *et al.* 2008), with one solution being
3.12
As expected, the magnitudes of the higher harmonics are very small,
3.13
for the parameter set (equation 3.12).

The parameter set in equation (3.12) is verified by direct numerical simulations of two coupled Brusselators (equation (3.11) with *N*=2), starting from random initial conditions. The feedback parameters of equation (2.3) were set to be (*k*_{1},*k*_{2},*τ*_{1},*τ*_{2},Δ*τ*)=(2.56,4.68,5.34,1.98,*ω*^{−1}*α*), where *ω* is the natural frequency of the oscillator and *α* is a variable control parameter. After a short transient period, the phase difference between the elements achieves a nearly constant value. Figure 3 displays Δ*ϕ** for 0<Δ*ϕ**≤*π* and 2*π*−Δ*ϕ** for *π*<Δ*ϕ**≤2*π*, which shows excellent agreement with theoretical predictions. Other parameter sets found numerically from equation (2.2) give quantitatively similar results (data not shown).

## 4. Experimental demonstration of phase difference tuning

### (a) Experimental setup

Experiments were conducted using an apparatus consisting of an electrochemical cell, a set of zero resistance ammeters (ZRAs) and a real-time data acquisition system (figure 4). The cell consisted of two Ni electrodes (99.99% pure, 1 mm diameter) in a 3 M H_{2}SO_{4} solution, a Pt mesh counter electrode and an Hg/Hg_{2}SO_{4}/K_{2}SO_{4} (sat) reference electrode. The cell was enclosed in a jacketed glass vessel held at a constant temperature of 11^{°}C. An EG&G potentiostat was used to adjust the circuit potential, *V* _{0}, of the cell, causing the nickel electrodes to undergo transpassive dissolution. The dissolution current of each electrode, *I*_{m}(*t*), was measured by a ZRA. A 650 Ω resistor (*R*_{p}) was attached to each channel to induce oscillations in the electrode potential (Zhai *et al.* 2004). A Labview-based real-time data acquisition computer was used to acquire data from the ZRA, send the data to a host machine, calculate a global feedback signal from these data and apply the feedback signal to the potentiostat at 250 Hz. The measurements of dissolution current were scaled such that the mean value of each channel () was removed and then normalized by the amplitude of its oscillation (*A*_{m}) relative to the mean amplitude of the population (*A*_{mean}):
4.1
The host machine was used to continuously determine the offset and amplitude of each rhythmic element in the population over time. The feedback signal, *δV* , fed back to the potentiostat was calculated using the equation
4.2
where *K* is the overall feedback gain (i.e. interaction strength), *h*(*x*) is equation (2.3) and *x*_{m}(*t*) is the potential drop across the double layer for the *m*th electrode, calculated as
4.3
where *V* (*t*) is the applied voltage and *R*_{p} is the channel resistance.

Before the nickel electrode array was placed into the system, it was polished with a wet rotary polisher to remove any initial oxide layer that may have been present. Six polishing discs were used, decreasing in roughness from 180 to 4000 grit. This ensured that the electrodes started from roughly identical conditions and that these initial conditions were approximately identical for each experiment. After polishing, the electrodes were placed in the acid solution and the potential was ramped from −0.68 to 1.25 V and back to 0 V without any resistors present. This caused a thin passive oxide layer to form on each electrode. After this cycle was complete, the resistors were reconnected and the system was brought up to the desired operating voltage. The system was allowed to line out at this voltage for approximately 1.5–2.0 h, as this was found to reduce the amount of drift in the inherent frequency of the oscillators. The electrode potential waveform of both elements can be seen in figure 5*a*. In the absence of feedback, the phase difference between the two elements increases over time (figure 5*b*).

### (b) Determining the response function

The response function can be determined from equation (2.2), provided that the interaction function and feedback function are both known. Work by Miyazaki & Kinoshita (2006) was extended to provide a method for measuring the interaction function of a rhythmic system composed of two non-identical oscillators under weak global feedback. Physically, the interaction function represents the change in the frequency of a rhythmic element as a function of its phase relative to all other interacting elements in the system. Therefore, the interaction function can be determined by recording the frequency (or period) of two interacting oscillators as a function of their phase difference. Two electrochemical oscillators were created, with different inherent frequencies, and a weak feedback signal was applied to the system. The magnitude of the feedback was selected such that the system was unable to synchronize, allowing the phase difference between the elements to grow over time (figure 6*a*). The weak interaction caused the period of the two oscillators to fluctuate as the phase difference between the elements changed (figure 6*b*). By integrating the instantaneous frequency of one of the oscillators over the observed period of a single cycle and equating this to 2*π* rad (Miyazaki & Kinoshita 2006), the interaction function can be determined using the equation
4.4
where *P*_{base} is the intrinsic period of the oscillator and *K* is the overall feedback gain.

Measuring interaction functions under different feedback conditions allowed the response function (figure 6*d*) to be calculated from equation (2.2). These sets of data were used in a multi-objective optimization to find the single best response function that could reproduce all the measured interaction functions simultaneously. This fitting was necessary to reduce the effects of experimental noise in the interaction function measurements. The Fourier coefficients of the response function were used as the optimization parameters. The number of Fourier terms to be optimized was determined by the number of higher harmonics present in the waveform and the measured interaction functions. Close to the Hopf bifurcation point (approx. 1.105 V), only one to two terms were necessary since the system was largely sinusoidal. At higher circuit potentials (approx. 1.20 V), the system became relaxational and the response function required approximately seven terms. At the experimental circuit potential of 1.165 V, four Fourier terms were used since the system was at a slightly higher voltage than the Hopf bifurcation point, but still exhibited relatively smooth oscillations. Additionally, it was found that increasing the number of coefficients did not significantly alter the shape of the response function. A simplex optimization algorithm was used to determine the value of the Fourier coefficients by minimizing an objective function of the form
4.5
where is the *i*th data point of the *n*th measured interaction function, *D* is the interaction function calculated using the optimized response function, *N* is the number of measured interaction function datasets and Pts is the number of measured points in the function. The initial conditions were taken to be the Fourier coefficients of a sine wave, since the response function approaches a sine wave as the system approaches the Hopf bifurcation point. On average, the optimization required approximately 1200 iterations to converge, corresponding to approximately 20 min.

### (c) Calculating the feedback signal

The target interaction function (figure 7*a*) was engineered into the experimental system using the nonlinear time-delayed feedback defined in equation (2.3). The feedback parameters *k*_{n} and *τ*_{n} can be calculated using equation (2.2), provided that the response function of the system and the target interaction function are known. The common feedback delay, Δ*τ*, was set to zero for this calculation. A simplex optimization algorithm was used to minimize the error between the calculated and target interaction functions by manipulating the feedback parameters. The objective function for this optimization was
4.6
where *T*_{i} is the *i*th data point of the target interaction function, *H* is the optimized interaction function and Pts is the number of data points in the functions *T* and *H*. Second-order feedback was used, since this is the highest order harmonic present in the target interaction function. The initial conditions for the optimization were taken to be *k*_{n}=10^{n−1} and *τ*_{n}=0.01 rad/2*π*. This allows each polynomial feedback term to be of the same order of magnitude, since |*x*|<1. The optimization successfully found a feedback parameter set that produced the target interaction function with sufficient accuracy:
4.7

Figure 7*a* shows the optimized interaction function (solid line) compared with the target interaction function (dashed line). There is good agreement between the target interaction function and the optimized interaction function. Since the optimized function has a large domain of negative slope, it will produce a large non-trivial synchronization region.

The feedback parameters were experimentally validated by measuring the interaction function produced by the optimized feedback conditions. Figure 7*b* illustrates the agreement between the predicted interaction function and the measured interaction function, indicating that the optimized parameters were successful in producing the desired interaction function.

### (d) Experimental results

The feedback parameters (equation (4.7)) were applied to the experimental system and the common feedback delay, Δ*τ*, was adjusted to achieve the desired phase difference between the two rhythmic elements. The phase of an element was linearly interpolated between adjacent peaks in the observed waveform; the phases of the peaks were defined as 0 and 2*π* rad, respectively. When Δ*τ*=0 rad/2*π*, the elements phase-synchronized with Δ*ϕ**=0 rad/2*π* (figure 8*a*,*d*). Under these conditions, the system exhibited only one stable stationary state, corresponding to the root of *H*^{−}(Δ*ϕ*) located at Δ*ϕ*=0 (figure 8*g*). Increasing Δ*τ* to 0.23 rad/2*π* caused the interaction function to shift, changing the synchronized phase difference to Δ*ϕ**=0.73 rad/2*π* (figure 8*b*,*e*). In this case, the system exhibited bi-stability, in which two stable stationary states (Δ*ϕ**=0.73 and Δ*ϕ**=0.28 rad/2*π*) coexist simultaneously (figure 8*h*). Further increasing Δ*τ* to 0.5 rad/2*π* caused the elements to synchronize in an anti-phase configuration, where Δ*ϕ**=0.5 rad/2*π* (figure 8*c*,*f*). This anti-phase configuration corresponded to the root of *H*^{−}(Δ*ϕ*) located at Δ*ϕ*=0.5 rad/2*π* (figure 8*i*), which was previously unstable.

A more quantitative comparison between the phase model predictions and experimental results can be seen in figure 9. The bifurcation diagram in figure 9*a* was generated by tracking the roots of *H*^{−}(Δ*ϕ*) (and their associated derivatives) as a function of Δ*τ*. Therefore, the diagram represents all the possible stationary states (both stable and unstable) for the experimental system at different values of Δ*τ*. To determine the accuracy of these predictions, the stationary phase difference between two electrochemical elements was experimentally measured with different values of Δ*τ*. These measurements were made by first synchronizing the system on the upper non-trivial branch of the bifurcation diagram using the appropriate delay and initial conditions. When the initial phase difference between the elements was 0.5<Δ*ϕ*<1 rad/2*π*, the system approached the upper bifurcation branch, otherwise it approached the lower branch. Once on the upper branch, Δ*τ* was adjusted to probe the extent of the branch. The global feedback was not removed during this process, ensuring that the system remained on the upper branch during these experiments. Once the upper branch was mapped, the feedback was removed and the system was allowed to drift into the range of initial conditions required to reach the lower branch (0<Δ*ϕ*<0.5 rad/2*π*). Upon reaching these conditions, the feedback was reapplied and the lower branch was mapped using the same process.

The unstable stationary states were experimentally probed using the negative feedback; the overall feedback gain, *K*, was multiplied by −1, while all the other feedback parameters remained constant. According to the linear stability analysis of equation (3.2), changing the sign of *K* inverts the stability of the fixed points. Therefore, the unstable states under positive feedback become stable under negative feedback, allowing for experimental observation. Excellent agreement was found between experimental measurements and theoretical predictions for both the stable and unstable stationary states (figure 9*b*).

## 5. Conclusions

We have demonstrated the application of a synchronization methodology by controlling the stationary phase difference between two oscillators with nonlinear feedback; we demonstrate the method in both numerical simulations and chemical experiments. Our methodology is based on the fact that the existence and stability conditions of dynamical states in weakly coupled oscillators are characterized by the interaction function in the phase model. Thus, by knowing the interaction function that results in a desired collective behaviour, the only remaining issue is how to construct the physical interaction that yields the phase interaction function. An interaction function can be constructed using the proposed feedback function (equation (2.3)). The choice of the specific form of the feedback function was motivated by the flexible application of the imposed interaction function for synchronization engineering. It has been shown that the *n*th order term of the feedback signal effectively enhances the *n*th Fourier components of the interaction function. The time delay of the *n*th term is used to adjust the balance of even and odd parts of the Fourier components. As a result, all the harmonics in the interaction function can be tuned by adjusting the appropriate polynomial gain and time delay (Kori *et al.* 2008).

The proposed methodology is applied to tune the phase difference of uncoupled oscillators through feedback perturbations. The methodology is still applicable when the oscillators are relatively weakly coupled (even with time varying strengths) because the feedback signal can dominate the relatively weak inherent coupling effects within the limitations of the phase model description; a similar approach was previously used (Kiss *et al.* 2007) for the problem of desynchronization of a population of coupled oscillators. Extension of the proposed method is required for the application to strongly interacting oscillators and to a population of oscillators, especially when there is a range of inherent coupling strengths. Nonetheless, the general framework proposed here should serve as a guideline for tackling these important problems that are expected to arise in biology.

A major advantage of this methodology is that the feedback, which results in a target interaction function, can be designed through the knowledge of the macroscopic observable of an isolated oscillator, that is, the waveform and the phase response function. This point is beneficial when applications to biological systems are considered; it is usually a formidable task to construct an appropriate detailed mathematical model of a biological system; however, the investigation of the phase response function is often possible.

## Acknowledgements

This work was supported in part by the National Science Foundation through grant CBET-0730597. I.Z.K. thanks the organizers for financial aid to attend ECC10. H.K. acknowledges financial support from the Grants-in-Aid for Young Scientists (no. 19800001) and the Sumitomo Foundation (no. 071019).

## Footnotes

One contribution of 10 to a Theme Issue ‘Experiments in complex and excitable dynamical systems’.

- © 2010 The Royal Society