## Abstract

We used phase models to describe and tune complex dynamic structures to desired states; weak, nondestructive signals are used to alter interactions among nonlinear rhythmic elements. Experiments on electrochemical reactions on electrode arrays were used to demonstrate the power of mild model-engineered feedback to achieve a desired response. Applications are made to the generation of sequentially visited dynamic cluster patterns similar to reproducible sequences seen in biological systems and to the design of a nonlinear antipacemaker for the destruction of pathological synchronization of a population of interacting oscillators.

Complex system responses can emerge from interactions among nonlinear rhythmic components (*1*, *2*). Examples abound in biology (*3*, *4*), communications (*5*), population dynamics (*6*), and chemical reaction systems (*7*, *8*). Inherent feedback is often an integral component; for example, circadian rhythms are controlled by the interactions within the multi-cellular master circadian clock in the brain, entrainment from sunlight, and feedback from other brain parts and locomotive activities (*9*).

External feedback can be used to control the behavior of complex rhythms, both to tune essential behavior, such as by heart pacemakers (*10*), or to alter pathological behavior, such as by deep-brain “antipacemakers” in tremors or Parkinson's disease (*11*). In such applications, a mild control is desired so that the system can be tuned to a desired behavior without destroying its fundamental nature. The efficient description and design of complex dynamic structure is a formidable task that requires simple yet accurate models, incorporating integrative experimental and mathematical approaches that can handle hierarchical complexities and predict emergent, system-level properties. Such approaches include phase models (*3*, *12*, *13*) and pulse-coupled models (*14*) that have been used to describe mutual entrainment of weakly interacting neuronal assemblies (*15*–*17*).

Here we present a methodology for the design of complex dynamical structure that does not require a detailed chemical or biological description of the individual units. The simplicity and analytical tractability of phase models (*3*, *12*, *18*, *19*) are exploited to design optimal global, delayed, nonlinear feedback for obtaining and tuning desired behavior. The feedback design methodology is capable of creating a large class of structures describable by phase models for general self-organized rhythmic patterns in weakly interacting systems with small heterogeneities. The method is demonstrated in three experiments: the tuning of desired arbitrary phase differences between two dissimilar oscillators, the generation of complex patterns that include self-organized switching between unstable dynamical states, and the physiologically important problem of desynchronization of oscillators.

We engineer a desired behavior of a population of *N* oscillators through the imposition of a nonlinear, time-delayed feedback. A time-dependent system parameter perturbation [δ*p*(*t*)], where *p* is a system parameter and *t* is time, is chosen to be a nonlinear function of the measured variables *x*_{k}(*t*) summed over the population (1) where *K* is the overall gain; here, we choose the nonlinear feedback function (*h*) to be a polynomial (2) where *k*_{n} and τ_{n} are the gain and the delay of the *n*th-order feedback, respectively, and *S* is the overall order of the feedback.

The challenge is obtaining the best form of the feedback: that is, obtaining the order and time delays best suited for the desired complex structure. We exploit the simplicity and flexibility of phase models in describing collective behavior of complex rhythms. Phase models, which describe each rhythmic unit by a single variable, the phase, provide an efficient tool for the analysis of the collective behavior of coupled oscillators (*3*, *12*). Because of their simplicity and mathematical tractability (as compared with full ordinary differential equation descriptions), phase models can often be constructed to yield typical dynamical states in populations of interacting rhythms including, for example, synchronized or desynchronized behavior, stable or intermittent clustering, or bistability between synchronized and nonsynchronized states (*3*, *12*, *18*–*21*).

A phase model is constructed that reproduces the desired state. A population of oscillators with weak, global (all-to-all) coupling can be described by (*3*, *12*) (3) where φ_{i} and ω_{i} are the phase and the natural frequency of the *i*th oscillator, *K* is the global coupling strength, and *H* is the interaction function. [Such a reduction is possible for a weakly heterogeneous population where the heterogeneities are small as compared to the intensity of coupling (*12*).] Equation 3 shows that the phase of an element increases at a rate equal to its inherent frequency (ω_{i}), slightly modified by slowing down or speeding up resulting from interactions with other elements. The interaction function *H*(Δφ) characterizes the extent of phase advance or delay as a result of the interaction between oscillators. For a desired target state (i.e., time variation of the phases of the oscillators), an optimal target interaction function *H*(Δφ) is determined through analytical and numerical investigations of the phase model (Eq. 3). The manipulation of the harmonics in the interaction function provides flexibility in the development of desired states.

After obtaining an appropriate interaction function *H*(Δφ) for the phase model, the feedback parameters for use in the experiments (Eqs. 1 and 2) can be obtained from the relation (*12*) (4) where the response function *Z*(φ), proportional to the phase-response curve widely used to interpret external entrainment in circadian rhythms, shows the phase advance per unit perturbation as a function of the phase of the oscillator. Consequently, we obtain *Z*(φ) (*22*), and thus *H*(Δφ), from direct experiments on one (*13*, *15*, *23*) or two (*24*) oscillators.

Given a feedback δ*p*(*t*) and a response function *Z*(φ), we could obtain the interaction function *H*(Δφ) for use in the phase model. However, we proceed in the opposite manner and choose an interaction function to produce desired states and then design a feedback loop δ*p*(*t*) with optimized feedback gains *k*_{n} and delays τ_{n} to give the desired *H*(Δφ). The parameters *k*_{n} and τ_{n} are found with standard optimization techniques (*22*). It can be shown analytically that, in weakly nonlinear oscillators, the order of feedback enhances the corresponding harmonic in the interaction function, and the delay time produces an offset in the phase difference (*22*). Thus, if we need an interaction function with predominantly first- and second-order harmonics, linear and quadratic feedback shall be applied and the delays of the feedback used to tune the ratio of the cosine and sine terms of *H*. The optimized feedback is then expected to produce the target dynamics through imposing the proper interaction function in the phase-model description.

We first demonstrate the method with a seemingly simple, but nevertheless nontrivial, example: tuning the (phase-locked) phase difference between two electrochemical oscillators with different inherent frequencies. The problem may arise in the design of oscillator arrays for communications and radar applications (*25*) and also in the dephasing of the timing of spikes in neurons (*20*). The electrode potential *E*(*t*) in the experimental chemical system is oscillatory (Fig. 1A). The phase difference between two (non-interacting) electrodes of somewhat different frequencies changes approximately linearly (Fig. 1B). For illustration, we choose the special case of an out-of-phase entrained target state with phase difference Δφ*= π/2 that can be obtained with an interaction function with first- and second-order harmonics *H*(Δφ) = cos(Δφ) - sin(2Δφ). [A general target phase difference can be obtained with an interaction function having an odd part *H*(Δφ)=sin(Δφ)+ *R*sin(2Δφ), where *R* is a control parameter (*22*).] Because the target interaction function is composed of first- and second-order harmonics, a feedback composed of linear and quadratic terms is chosen; the design requires the waveform of the oscillators and the response function (Fig. 1C). The optimized and target interaction functions and odd part are shown in Fig. 1, D and E. With feedback, the desired locked phase difference of Δφ*= π/2 (or –π/2) (Fig. 1, B and F) was achieved.

We now consider the generation of sequential dynamical states (*26*, *27*). Sequential patterns can play a role in information processing and in the functioning of memory in neural systems (*26*, *27*); scent cues processed by the olfactory system are encoded in complex spatial and temporal patterns of firing neurons (*28*). The mathematical concept of slow switching (*21*) predicts an alternation between synchronized cluster states in a population of (at least) four oscillators with, for example, *H*(Δφ) = sin(Δφ – 1.32) - 0.25 sin(2Δφ) in Eq. 3. Because heteroclinic orbits connect the unstable dynamic states and these orbits are typically not robust against heterogeneities and noise caused by their structural instabilities, their demonstration in experimental systems is a challenging task.

We optimized a quadratic feedback to a population of four oscillators that reproduced an interaction function proposed for slow switching (Fig. 2A). The experimental system with feedback sequentially visits (unstable) two-cluster states with two oscillators in each cluster; Fig. 2B shows two (saddle-type) cluster states in state space. The phase model predicts a switch between these states as a result of the existence of heteroclinic orbits connecting the states. In the experiments, we observed switching between the cluster states (red line in Fig. 2B) along the theoretically predicted orbit (black line in Fig. 2B). We observed many switches along the heteroclinic orbits in a long time series. These switches can be seen as a fluctuation of the system order (Fig. 2C). The time scale of the cluster switching is 60 s, much greater than the 2.2-s period of the individual oscillating elements.

The engineered feedback produces configurations of two clusters, each containing two elements, connected by heteroclinic orbits; the particular elements in each cluster and the switching dynamics depend on initial conditions as well as heterogeneities and noise. Figure 2D shows that, in a long experiment, the system iterates among predicted two-cluster configurations; when the orbits approach a cluster state, the (color-coded) phase velocity slows down, indicating that the states are of saddle types. During the iteration among orbits, the system exhibits a complicated pattern with occasional jumps from one heteroclinic orbit to another. Two types of transitions have been experimentally observed. Intracluster transitions occur when the elements of a single cluster reverse themselves (compare *t* = 60s and *t* = 180 s in Fig. 2C). An intercluster transition occurs when an element from one cluster pairs with an element from the adjacent cluster [fig. S4 (*22*)]. Both types of transitions are seen in the trajectory of the experimental system illustrated in phase-space plots in Fig. 2, D and E. Because of inherent heterogeneities (about 0.2% frequency differences) between the oscillators, a motion in state space close to but not on the heteroclinic orbits predicted by the (homogeneous) phase model is seen both in the experiments and in phase-model simulations with heterogeneous oscillators (red surface in Fig. 2E).

Synchronization in biological systems can be pathological. Deep-brain stimulation with high-frequency, high-amplitude electrical signals is being used to destroy synchronized rhythms in Parkinson's disease and essential tremor and has potential application in epilepsy antipacemakers (*11*). Low side-effect antipacemakers require the development of advanced desynchronization methods (*29*, *30*).

Our proposed phase-model methodology provides an efficient design of mild nonlinear feedback antipacemakers for weakly interacting systems. A system of 64 weakly relaxation oscillators, synchronized with global coupling through a common resistor (*13*), has the interaction function shown in Fig. 3B obtained from the response function and waveform shown in Fig. 3A. The interaction function exhibits a positive slope at Δφ = 0, resulting in a synchronized (one-cluster) state with a large order parameter (left side of Fig. 3C). Although the one-cluster state can be broken with a linear feedback technique (*29*), instead of a desynchronized state, synchronized cluster states can appear (here as a three-cluster configuration as seen in Fig. 3C). The occurrence of such clusters is caused by the higher harmonics in the overall interaction function, including the aggregate effects of coupling and linear feedback.

A desynchronized state without stable, ordered cluster states can be obtained with nonlinear feedback. Many nonlinear feedbacks can produce a nonsynchronized state without clusters. Mild, effective desynchronization can be achieved by minimizing the power of the feedback signal under the condition that the feedback produces a family of target interaction functions with negative odd components: for example, , where *M* is the largest harmonics considered and the ϵ_{k} are small numbers. A linear programming optimization (*22*) resulted in a mild, second-order feedback that produces an interaction function with negative odd harmonics (up to *M* = 5) (Fig. 3D). This quadratic feedback can successfully desynchronize the system as seen in Fig. 3E. The initially synchronized state (*t* < 60s) desynchronizes with the designed nonlinear feedback (60 s < *t* < 440 s); the elements almost uniformly populate the cycle, and all order parameters drop to low values. When the feedback is turned off (*t* = 440 s), the system returns to the synchronized state.

We have experimentally demonstrated an effective method of designing complex dynamic structure and tuning a wide spectrum of emergent collective behavior in systems composed of rhythmic elements. The method is precise enough to engineer delicate synchronization features of nonlinear systems and can be applied to both small sets and large populations. Because the method does not require a priori detailed physical, chemical, and biological models, it should find applications in pacemaker and anti-pacemaker design in systems where there is a need for tuning complex dynamical rhythmic structures but where such detailed models are difficult to obtain.

**Supporting Online Material**

www.sciencemag.org/cgi/content/full/1140858/DC1

Materials and Methods

Figs. S1 to S4

References