Abstract
The ability to organize and finely manipulate the hierarchy and timing of dynamic processes is important for understanding and influencing brain functions, sleep and metabolic cycles, and many other natural phenomena. However, establishing spatiotemporal structures in biological oscillator ensembles is a challenging task that requires controlling large collections of complex nonlinear dynamical units. In this report, we present a method to design entrainment signals that create stable phase patterns in ensembles of heterogeneous nonlinear oscillators without using state feedback information. We demonstrate the approach using experiments with electrochemical reactions on multielectrode arrays, in which we selectively assign ensemble subgroups into spatiotemporal patterns with multiple phase clusters. The experimentally confirmed mechanism elucidates the connection between the phases and natural frequencies of a collection of dynamical elements, the spatial and temporal information that is encoded within this ensemble, and how external signals can be used to retrieve this information.
Introduction
Complex interactions among nonlinear periodic phenomena emerge in many natural and engineered systems^{1,2}. Numerous instances appear in chemical reactions^{3,4} and biological systems^{5,6}, which exhibit endogenous and emergent multiscale oscillations^{7}. There is significant interest in characterizing synchronization in oscillators interconnected in networks^{8,9}, which is especially important for understanding the highly complex dynamics of manmade systems such as electric power grids^{10}, and elucidating the functions of neural systems^{11,12,13}. Understanding entrainment of oscillating systems to an exogenous forcing signal is crucial to modelling circadian timekeeping^{14}, dynamic neural regulation^{15} and for the design of synchronizing or desynchronizing treatments of cardiac arrhythmias^{16}, Parkinson’s disease^{17}, epilepsy^{18} and movement disorders^{19}.
Even the simplest models of interacting oscillators can exhibit highly complex behaviour^{20}, and individual oscillating units may themselves possess complicated dynamics^{21}. These factors are aggravated in practice by model and parameter uncertainty and the impracticality of obtaining feedback information, such as for in vivo biological applications, and pose challenges to manipulating or controlling oscillating ensembles. Such tasks require tractable yet accurate simplifications of the complex dynamic interactions involved, and demand suitable mathematical approaches that characterize ensemblelevel properties while withstanding experimental uncertainties.
Controltheoretic techniques have been applied to control a single oscillator^{22,23,24}. In contrast, finely manipulating individual subsystems in underactuated ensembles, such as thousands of neurons in the brain affected by one electrode, rather than activating them homogeneously remains a fundamental challenge. Synchronization has been engineered in collections of oscillators using feedback^{25,26,27}, or tuning coupling strengths^{4,28,29}. Such approaches require certain coupling structures, exact model specification, state feedback information, or precise knowledge of initial conditions, but still are not able to produce a prescribed phase pattern corresponding to frequency clusters of the oscillators.
Versatile openloop control techniques were developed for simultaneous control of ensembles of quantum spin systems, which motivated the field of ensemble control^{30}. Inspired by selective pulse design in nuclear magnetic resonance (NMR)^{31}, which enabled revolutionary applications including functional magnetic resonance imaging (fMRI), we develop a method for selectively manipulating the subunits of oscillator ensembles using periodic inputs that are robust to parameter uncertainty and disturbances. Specifically, we exploit the slight heterogeneity and high nonlinearity of an ensemble of structurally similar oscillators far past the Hopf bifurcation, rather than relying on a known coupling structure, state feedback or initial condition information.
In this manuscript, we present a methodology for constructing weak, globally applied, openloop control inputs that synchronize a collection of structurally similar yet heterogeneous nonlinear oscillators while selectively assigning their relative phases on the periodic orbit. Using the technique, the synchronization structure of an oscillating ensemble can be manipulated among diverse phase patterns, seen in relative positions on the limit cycle. Our theory is developed specifically to overcome the challenges of experimental implementation when feedback information is unavailable, initial conditions are unknown and the oscillators are subject to uncertainty in subsystem parameters and stochastic disturbances. The control inputs create and maintain such phase patterns when the coupling between oscillators is weak, while preserving the intrinsic nature of the ensemble to enable nondestructive application to fragile biological and chemical systems. The dynamics of the oscillators may be arbitrary, as long as all are structurally similar and exhibit sufficient nonlinear relaxation for the control design to be realizable. A coherent structure may be established and robustly maintained indefinitely by a single periodic waveform, which can be altered to switch between patterns. We demonstrate the theoretical methodology in practice with experiments to control complex electrochemical ensembles whose dynamics are variable and unknown, and for which state information is unavailable^{26,32}.
Results
Phase model approximation
We approximate the effect of inputs on periodic dynamical systems using phase models^{33}, which can be computed for systems with known dynamics^{34} or experimentally inferred in practice when the dynamics are unknown^{35}. These models are used to characterize circadian cycles^{36}, cardiac rhythms^{16} and phenomena in neural and chemical systems^{11,37}, and their simplicity has enabled control design for neuron models^{22} given initial conditions and exact parameters. Control techniques have recently been developed for individual nonlinear oscillators and finite collections described by phase models that require exactly known initial conditions, parameters, and dynamics^{22,38,39}. Many studies on synchronization focus on the network structure of couplings between oscillators, and the nonlinearity in the phase response of each unit is simplified to sinusoidal couplings with its neighbours^{9}. However, for the manipulation and desynchronization of electrochemical and neural systems^{11,40,41}, complex, hierarchical interactions must be established or broken in large collections of nonlinear systems. The dynamics, parameters, and interconnections of these systems are typically problematic to infer, may be noisy, variable, or uncertain, and state observations may be incomplete or unavailable. Such conditions elude tractable formulation, and require an approach where synchronization properties of the systems, that is, asymptotic phase structure, are manipulated rather than steering the system states directly^{7,42}.
Entrainment of an ensemble
Our method relies on entrainment, which refers to the dynamic synchronization of oscillating systems to periodic inputs. Each system in an ensemble of N structurally similar units exhibits endogenous oscillation along an attractive periodic orbit with period T_{j}, and is represented by the Winfree phase model
where ω_{j}=2π/T_{j} is the natural frequency and Z is the phase response curve (PRC), which quantifies how a weak perturbation u advances or delays the phase (refs 2, 33). A value of (equivalently 2π) corresponds to a measurement of the jth system reaching its maximum. More details on phase coordinate transformation are given in Supplementary Note 1. We demonstrate our phaseselective entrainment techniques using experiments in which nickel electrodes undergo anodic dissolution in sulfuric acid and exhibit electrochemical oscillations^{32}, for which the experimental apparatus is described in section ‘Experimental apparatus’. Before these experiments, the PRCs of the ensemble elements, which are shown in Fig. 1a, were estimated and averaged for use as the nominal PRC in equation (1). This was done using a pulse perturbation procedure for system identification that was previously used for electrochemical oscillators^{43}, and is described in Supplementary Note 2.
To synchronize oscillation of ensemble elements, each subsystem receives the same weak, periodic forcing input of frequency Ω of the form u(t)=v(Ωt), where v is 2πperiodic. When the forcing frequency is near the natural frequencies in the ensemble, averaging theory^{44} states that the mean difference between each phase and the forcing phase θ=Ωt follows the timeindependent dynamic equation
where Δω_{j}=ω_{j}−Ω is called the frequency detuning, and
is a 2πperiodic interaction function that characterizes the average effect of the periodic input on the oscillation phase^{43}. Such ergodic averaging is discussed in more detail in Supplementary Note 3 and illustrated in Supplementary Fig. 1. If equation (2) has a unique attractive fixed point that satisfies , then the phase of the jth oscillator becomes entrained to the forcing phase θ with an average offset of . This analysis is widely applied to determine the interaction function resulting from a measured PRC and an input waveform, and equation (2) is used to infer the entrained system’s stability^{27}. We reverse this approach by choosing the desired asymptotic behaviour, constructing a suitable interaction function and using the PRC to obtain the input by circular deconvolution of equation (3).
Interaction functions for phase selection
We design the input v(Ωt) so that each system in equation (1) is entrained to a forcing frequency Ω, e.g., the mean of natural frequencies ω_{1}<ω_{2}<…<ω_{N}, and such that the jth oscillator cycles its orbit with a phase offset of relative to the forcing phase θ. The set of pairs constitutes a pattern for selective entrainment. We require that eventually holds for each oscillator, that is, equation (2) has an attractive fixed point at for all j at which the slope of the interaction function Λ_{v} is negative^{43}. The function that best satisfies these ideal conditions has steep decreases at phase values where entrainment must occur, and crosses (from above) horizontal lines at frequency detuning values −Δω_{j}. This creates the desired attractive fixed points for equation (2). Because the interaction function is periodic, it must then increase so that Λ_{v}(2π)=Λ_{v}(0) holds. Crossings of −Δω_{j} from below are unstable fixed points, and do not affect convergence.
The concept is illustrated in Fig. 1, which describes an experiment where a phase difference is assigned between two entrained oscillators. In Fig. 1c we desire inphase synchronization at phase offsets of , so the ideal interaction function has one steep decrease at that intersects horizontal lines through −Δω_{1} (blue) and −Δω_{2} (red) at π radians. Figure 1d illustrates antiphase synchronization with , where the interaction function has two steep decreases at and that intersect horizontal lines at −Δω_{1} (blue) and −Δω_{2} (red), respectively. The best achievable interaction function (solid line) and the PRC estimate yield the input from equation (3). The right columns of Fig. 1c,d show the observed current of two oscillators entrained in inphase and antiphase arrangements by the input waveform (shown above). These configurations are achieved regardless of initial oscillator phases, because the interaction function crosses the line −Δω_{j} only once from above, so each system has a globally attractive fixed point. For the electrochemical system, phase differences in nearly the entire 0 to 2π region are achievable, with small deviations as approaches 2π, as seen in Fig. 1b.
Separation of oscillator ensembles into phase clusters
Uniquely attractive phase patterns are desired, where a common input synchronizes the oscillators to a pattern independently of their initial phases. The fixed point of equation (2) must be unique for each j, which is achieved when the interaction function crosses each horizontal line Δω_{j} from above only once at . This is possible when the phase offsets are monotonically ordered as for ω_{1}<ω_{2}<…<ω_{N}, as demonstrated by segregation of 20 inhomogeneous electrochemical oscillators into clusters in the experiments described in Fig. 2. An antiphase configuration is achieved for electrodes in balanced (N_{1}, N_{2})=(10, 10) and unbalanced (N_{1}, N_{2})=(1, 19) clusters in Fig. 2a,b, respectively. In these twocluster examples, the interaction function decreases in two steps, of which the top and bottom correspond to clusters with slower (blue) and faster (red) natural frequencies. Figure 2c shows the formation of four balanced clusters of (N_{1}, N_{2}, N_{3}, N_{4})=(5, 5, 5, 5) oscillators with the phase structure radians. The phase offsets are increasing as −Δω_{j} decreases (and ω_{j} increases), and the designed interaction function decreases monotonically as it crosses the required intersections. Observe that the assumption of a common PRC is reasonable, because the functions in Fig. 1a are very similar, yet the oscillator frequencies are sufficiently heterogeneous for our technique to create the patterns in Fig. 2.
Control of pattern transitions in an ensemble
In addition, we establish patterns without monotone phase ordering by designing an interaction function of the form at the bottom of Fig. 3d, which crosses a horizontal line at −Δω_{j} from above twice, yielding multiple possible entrained phases and dependence on initial conditions. We apply precursor waveforms to steer subsets of the ensemble into attractive regions for the desired phase offsets , then finalize and hold the pattern with an ultimate input. This procedure is applied to steer an ensemble between spatially associated clusters by alternating selective inputs. Fig. 3 illustrates input design for itinerant formation of letters in the word ‘OK’ in the array of 20 electrochemical oscillators used in the experiments in Fig. 2. We produce antiphase assignment between clusters to display the letter ‘O’, then switch the input to produce the letter ‘K’. Rhythmic elements are assigned desired phase offsets of or , which correspond to ‘on’ (pattern) or ‘off’ (background) states, respectively, that are visualized using a colour scale where 0 (‘on’) is blue and π (‘off’) is yellow. Switching between two patterns is accomplished using four numbered clusters, where 1 is ‘on’ for both, 2 switches from ‘on’ to ‘off', 3 switches from ‘off’ to ‘on’, and 4 is always ‘off’. Electrodes in clusters of (N_{1}, N_{2}, N_{3}, N_{4})=(7, 7, 3, 3) elements with mean natural frequencies (ω_{1}, ω_{2}, ω_{3}, ω_{4})=(0.390, 0.406, 0.427, 0.442) Hz are positioned in the spatial arrangement in Fig. 3a. Figure 3b–d each have panels that show, from top to bottom, the spatial distribution of phase offsets, the structure on the unit circle, and a sketch of the ideal interaction function.
The pattern ‘O’, shown in Fig. 3b, is realized using the phases , which are achieved by an interaction function as in Fig. 2a. The phases used for ‘K’ are not monotonically ordered, so a precursor waveform is applied to generate globally attractive phase offsets and for clusters 2 and 3. This antiphase pattern establishes initial conditions for the final input waveform, while clusters 1 and 4 lose their entrainment, as shown in Fig. 3c. Figure 3d illustrates the design of the finalizing control for pattern ‘K’, where the phase assignments for clusters 1 and 4 are globally attractive, as seen in the bottom panel of Fig. 3d, while clusters 2 and 3 stay at phase offsets established in the precursor stage. The transition from pattern ‘K’ back to ‘O’ is accomplished by applying the initial control for the pattern in Fig. 3b. We provide additional descriptions of interaction function construction (Supplementary Note 4; Supplementary Figs 1 and 2), pattern realizability (Supplementary Note 5), control design for monotonically ordered patterns (Supplementary Note 6; Supplementary Fig. 3), and precursor waveform engineering (Supplementary Note 7; Supplementary Figs 2 and 4).
Measurements of the ‘O→K→O’ pattern switching experiment appear in Fig. 4. Figure 4a shows current oscillations for the reaction units given zero input, when no pattern forms. When the controls (shown above the current) are applied, the ensemble is entrained within several cycles, selectively forming the patterns for ‘O’, the precursor, and ‘K’. These results are shown in Fig. 4b–d, and correspond to Fig. 3b–d, respectively. The ensemble is returned to pattern ‘O’, as shown in Fig. 4e, to demonstrate switching. The automatically generated interaction functions and control waveforms are presented in section ‘Automatic control waveforms generated in experiments’. An animation produced using the experimental current traces and oscillation phases is included as Supplementary Movie 1.
Discussion
Phaseselective entrainment enables the use of a single global signal to robustly assign elements of a noisy nonlinear oscillator ensemble to specific phases without using coupling or feedback information. Control design using interaction functions simplifies the creation of complex synchronization patterns to drawing or automatically generating curves through sets of crossing points and computing the resulting controls with a simple formula, which is an accessible technique for experimentalists. Greater relaxation in the oscillation and ensemble heterogeneity increases pattern controllability, and performance is improved as the oscillations move farther away from the Hopf bifurcation. The asymptotic nature of entrainment yields robustness to noise, disturbances and model parameter variability while preserving the intrinsic nature of the ensemble.
Such resilience is required for nondestructive control of underactuated, noisy and uncertain biological and chemical ensembles that cannot be readily observed. For example, an effective technology for neurological treatment of Parkinson’s disease^{17} is provided by deep brain stimulation, which alleviates pathological synchronization in the brain. Selective entrainment could be extended to ensembles with weak coupling to design robust desynchronization inputs, which would potentially benefit noninvasive neurostimulation technology^{41}. The goal could be a target distribution that is found to be optimal for leveraging neuroplasticity to prevent resynchronization after the stimulus is ended. The technique could also improve phase regulation to treat cardiac arrhythmias^{16} and sleep irregularities^{45}. The formalism could also represent the entrainment that occurs in circadian timekeeping^{14}.
We note that a simple sinusoidal forcing of the form v(Ωt)=sin (Ωt) results in a sinusoidal interaction function, because of orthogonality of the trigonometric Fourier basis. Sinusoidal forcing can thus be used to create monotone ordered phase patterns, and could also be used for desynchronization. However, because such an interaction function is decreasing on an interval of length π, the maximum achievable distance between extremal phase offsets and is . Thus, a sinusoidal input cannot produce antiphase synchronization. Our approach enables more versatile manipulation of phase relationships beyond this limitation. We describe the application to desynchronization in Supplementary Note 8 and Supplementary Fig. 5, and quantify how our approach increases the achievable relative phase desynchronization difference over sinusoidal forcing. More rigorous mathematical characterization of ensemble desynchronization by periodic inputs is a compelling direction for further investigation.
In our methodology, we take advantage of approximations that are possible in the specific experimental setting. In our experiments, the distribution of natural frequencies of ensemble oscillations varies by ±20% from the (nonzero) mean, the oscillators are very weakly coupled, the amplitude of the required external forcing signal is relatively small, and the entrainment process is approximated well by averaged phase models. In addition, although the ensemble subsystems are slightly heterogeneous and noisy, with variation in natural frequencies and dynamic properties, the phase response curves are very similar. We expect the methodology to function well in other experimental settings in which these conditions are satisfied. Moreover, the method holds promise for extension to other scenarios, e.g., in sub and superhamornic entrainments (the oscillations are locked to different frequency ratios), where the interaction functionbased phase description is possible.
The arrangement of frequency clusters in an oscillator ensemble can also be viewed as encoding information within the spatial pattern produced when selective entrainment is applied. One of several encoded messages can then be retrieved using the phaseselective entrainment process, for which the passkey for retrieving the message is the temporal information contained in a periodic input signal. The passkey is constructed using the PRC and natural frequencies of the dynamical subsystems, and after that input signal is applied, the spatial phase pattern produced in the ensemble reveals the message. This approach may be incorporated in neurocomputing architectures^{46} that mimic neural systems in nature. Future investigation is required to understand how network coupling could be suppressed or taken advantage of to improve pattern resilience and information capacity, or effectively encrypt the message by preventing estimation of PRCs and natural frequencies of oscillators in the spatial array.
Methods
Experimental apparatus
Our methodology for controlling the phase structure of an ensemble of heterogeneous oscillators is experimentally verified in the electrochemical dissolution of nickel in 3 mol l^{−1} sulfuric acid solution using potential actuation. A schematic description of the experimental setup is depicted in Fig. 5. The apparatus consists of 20 nickel wires that function as working electrodes (WE), with diameters of 0.69 mm, spaced by 2.0 mm and embedded in epoxy resin. Prior to the electrochemical measurements, the WE were polished with sandpaper in six levels of roughness, ranging from 180 to 4,000 grit. The current of all electrodes is monitored independently. Under such conditions, the spontaneous formation of selfsustained electrochemical oscillations driven by the presence of a negative differential resistance^{47} is observed in the anodic dissolution of nickel.
Once the WE are placed in the electrochemical cell, a slow positive sweep of 0.01 V s^{−1} from 0 V to a constantly applied potential V_{0} was conducted to form a thin passive layer on the electrode surface. This baseline is set in reference to an Hg/Hg_{2}SO_{4}/sat. K_{2}SO_{4} reference electrode (RE) in an electrochemical cell, containing a 1.6 mm diameter Pt coated Ti wire counter electrode at constant temperature of 10 °C. The potential V_{0} was initially set using a potentiostat (Bank Instruments) at a value for which the oscillation is close to the Hopf bifurcation, which is ∼1.15 V. Inputs to the oscillating system consist of an additional potential u superimposed onto the baseline potential V_{0} using the potentiostat.
Soon after, the PRCs were measured simultaneously for the WE by the automatic procedure of applying a pseudorandomly timed potential pulse sequence of −0.20 V with pulsewidth of 0.05 s and postprocessing the observed current using the pulse perturbation procedure as described in section ‘Entrainment of an ensemble’. Measurements of the current oscillation in the electrochemical reactions were carried out by a realtime data acquisition by a highspeed multifunction M Series DAQ PXI6255 (National Instruments) interface with a sample rate of 200 Hz. Simultaneously, the periodic potential perturbations u were applied in the electrochemical oscillator ensemble using the potentiostat to superimpose the waveform onto the applied baseline voltage V_{0}. Each control waveform was generated based on the PRC obtained preliminary to the experiment and the targeted interaction function generated from the desired set of phase offsets for the oscillators.
Automatic control waveforms generated in experiments
Figure 6 displays the interaction function (left column) and the respective control waveform (right column) used to entrain the electrochemical oscillators for the pattern switching procedure ‘O→K→O’ described in Fig. 4 (see main text for additional details). An electrochemical oscillator ensembles with clusters of size (N_{1}, N_{2}, N_{3}, N_{4})=(7, 7, 3, 3) were selected by tuning the mean natural frequency of the oscillations in the following order: cluster 1 with ω_{1}=0.390 Hz (electrodes j=1 to 7 in blue), cluster 2 with ω_{2}=0.406 Hz (electrodes j=8 to 14 in cyan), cluster 3 with ω_{3}=0.427 Hz (electrodes j=15 to 17 in yellow) and cluster 4 with ω_{4}=0.442 Hz (electrodes j=18 to 20 in red). The phase assignments for pattern switching ‘O→K→O’ are listed in Fig. 6.
Additional information
How to cite this article: Zlotnik, A. et al. Phaseselective entrainment of nonlinear oscillator ensembles. Nat. Commun. 7:10788 doi: 10.1038/ncomms10788 (2016).
References
 1
Strogatz., S. Nonlinear Dynamics And Chaos: With Applications To Physics, Biology, Chemistry, And Engineering Westview Press (2001).
 2
Pikovsky, A., Rosenblum, M. & Kurths., J. Synchronization: A Universal Concept in Nonlinear Science Cambridge Univ. Press (2001).
 3
Kuramoto., Y. Chemical Oscillations, Waves, and Turbulence Springer (1984).
 4
Epstein, I. R. & Pojman., J. A. An Introduction to Nonlinear Chemical Dynamics: Oscillations, Waves, Patterns, and Chaos (Topics in Physical Chemistry) Oxford Univ. Press (1998).
 5
Glass, L. Synchronization and rhythmic processes in physiology. Nature 410, 277–284 (2001).
 6
Glass, L. & Mackey., M. C. From clocks to chaos: The rhythms of life Princeton Univ. Press (1988).
 7
Zlotnik, A. & Li, J.S. Optimal subharmonic entrainment of weakly forced nonlinear oscillators. SIAM J. Appl. Dynam. Syst. 13, 1654–1693 (2014).
 8
Arenas, A., DazGuilera, A., Kurths, J., Moreno, Y. & Zhou, C. Synchronization in complex networks. Phys. Rep. 469, 93–153 (2008).
 9
Dörfler, F. & Bullo, F. Synchronization in complex networks of phase oscillators: A survey. Automatica 50, 1539–1564 (2014).
 10
Dörfler, F., Chertkov, M. & Bullo, F. Synchronization in complex oscillator networks and smart grids. Proc. Natl Acad. Sci USA 110, 2005–2010 (2013).
 11
Izhikevich., E. Dynamical Systems in Neuroscience MIT Press (2007).
 12
Uhlhaas, P. J. & Singer, W. Neural synchrony in brain disorders: relevance for cognitive dysfunctions and pathophysiology. Neuron 52, 155–168 (2006).
 13
Buzsáki, G. & Draguhn, A. Neuronal oscillations in cortical networks. Science 304, 1926–1929 (2004).
 14
Leloup, J.C. & Goldbeter, A. Toward a detailed computational model for the mammalian circadian clock. Proc. Natl Acad. Sci. USA 100, 7051–7056 (2003).
 15
Hunter, J. & Milton, J. Amplitude and Frequency Dependence of Spike Timing: Implications for Dynamic Regulation. J. Neurophysiol. 90, 387–394 (2003).
 16
Guevara, M. & Glass., L. Phase locking, period doubling bifurcations and chaos in a mathematical model of a periodically driven oscillator. J. Math. Biol. 14, 1–23 (1982).
 17
Hammond, C., Bergman, H. & Brown, P. Pathological synchronization in parkinson's disease: networks, models and treatments. Trends Neurosci. 30, 357–364 (2007).
 18
Good, L. Control of synchronization of brain dynamics leads to control of epileptic seizures in rodents. Int. J. Neural Syst. 19, 173–196 (2009).
 19
Reinkensmeyer, D. J. et al. Tools for understanding and optimizing robotic gait training. J. Rehab. Res. Dev 43, 657–670 (2014).
 20
Strogatz., S. H. From Kuramoto to Crawford: exploring the onset of synchronization in populations of coupled oscillators. Phys. D 143, 1–20 (2000).
 21
Izhikevich, E. M. Synchronization of elliptic bursters. SIAM Rev. 43, 315–344 (2001).
 22
Moehlis, H., Brown, E. & Rabitz, H. Optimal inputs for phase models of spiking neurons. J. Comput. Nonlin. Dyn. 1, 358–367 (2006).
 23
Dasanayake, I. & Li, J.S. Optimal design of minimumpower stimuli for phase models of neuron oscillators. Phys. Rev. E 83, 061916 (2011).
 24
Efimov, D., Sacré, P. & Sepulchre, R. Proceedings of the 48th IEEE Conference on Decision and Control, 7692–7697 (Shanghai, China, 2009).
 25
Sepulchre, R., Paley, D. A. & Leonard, N. E. Stabilization of planar collective motion: Alltoall communication. IEEE Trans. Autom. Control 52, 811–824 (2007).
 26
Kiss, I. Z., Rusin, C. G., Kori, H. & Hudson, J. L. Engineering complex dynamical structures: sequential patterns and desynchronization. Science 316, 1886–1889 (2007).
 27
Kori, H., Rusin, C. G., Kiss, I. Z. & Hudson, J. L. Synchronization engineering: Theoretical framework and application to dynamical clustering. Chaos 18, 026111 (2008).
 28
Rosenblum, M. G. & Pikovsky, A. S. Controlling synchronization in an ensemble of globally coupled oscillators. Phys. Rev. Lett. 92, 114102 (2004).
 29
Orosz, G. Decomposition of nonlinear delayed networks around cluster states with applications to neurodynamics. SIAM J. Appl. Dyn. Syst. 13, 1353–1386 (2014).
 30
Li, J.S. & Khaneja, N. Control of inhomogeneous quantum ensembles. Phys. Rev. A 73, 030302 (2006).
 31
Kobzar, K., Luy, B., Khaneja, N. & Glaser, S. J. Pattern pulses: design of arbitary excitation profiles as a function of pulse amplitude and offset. J. Magn. Reson. 173, 229–235 (2005).
 32
Kiss, I. Z., Zhai, Y. & Hudson, J. Emerging coherence in a population of chemical oscillators. Science 296, 1676–1678 (2002).
 33
Winfree., A. T. The Geometry of Biological Time SpringerVerlag (1980).
 34
Ermentrout, B. Type I membranes, phase resetting curves, and synchrony. Neural Comput. 8, 979–1001 (1996).
 35
Galan, R. F., Ermentrout, G. B. & Urban, N. N. Efficient estimation of phaseresetting curves in real neurons and its significance for neuralnetwork modeling. Phys. Rev. Lett. 94, 158101 (2005).
 36
Strogatz, S. H. Human sleep and circadian rhythms: a simple model based on two coupled oscillators. J. Math. Biol. 25, 327–347 (1987).
 37
Nakao, H., Yanagita, T. & Kawamura., Y. Phasereduction approach to synchronization of spatiotemporal rhythms in reactiondiffusion systems. Phys. Rev. X 4, 021032 (2014).
 38
Li, J.S., Dasanayake, I. & Ruths, J. Control and synchronization of neuron ensembles. IEEE Trans. Autom. Control 58, 1919–1930 (2013).
 39
Dasanayake, I. S. & Li, J.S. Constrained chargebalanced minimumpower controls for spiking neuron oscillators. Syst. Control Lett. 75, 124–130 (2015).
 40
Zhai, Y., Kiss, I. Z., Tass, P. A. & Hudson, J. L. Desynchronization of coupled electrochemical oscillators with pulse stimulations. Phys. Rev. E 71, 065202 (2005).
 41
Tass., P. A. A model of desynchronizing deep brain stimulation with a demandcontrolled coordinated reset of neural subpopulations. Biol. Cybern. 89, 81–88 (2003).
 42
Zlotnik, A. & Li, J.S. Optimal entrainment of neural oscillator ensembles. J. Neural Eng. 9, 046015 (2012).
 43
Zlotnik, A., Chen, Y., Kiss, I. Z., Tanaka, H.A. & Li, J.S. Optimal waveform for fast entrainment of weakly forced nonlinear oscillators. Phys. Rev. Lett. 111, 024102 (2013).
 44
Hoppensteadt, F. & Izhikevich., E. Weakly connected neural networks SpringerVerlag (1997).
 45
Gooley, J. J. et al. Spectral responses of the human circadian system depend on the irradiance and duration of exposure to light. Sci. Trans. Med. 2, 31ra33 (2010).
 46
Hoppensteadt, F. & Izhikevich, E. Oscillatory Neurocomputers with Dynamic Connectivity. Phys. Rev. Lett. 82, 2983 (1999).
 47
Koper., M. T. M. The theory of electrochemical instabilities. Electrochim. Acta. 37, 1771–1778 (1992).
Acknowledgements
This work was supported by the U.S. National Science Foundation under the awards CMMI1301148, CMMI1462796, ECCS1509342, and CHE1465013, the CNPq Award 229171/20133, and the U.S. Department of Energy under contract DEAC5206NA25396.
Author information
Affiliations
Contributions
The authors contributed equally to developing the methodology presented here. A.Z. developed computational algorithms and did primary writing for the manuscript and R.N. performed the experiments and created figures.
Corresponding author
Ethics declarations
Competing interests
The authors declare no competing financial interests.
Supplementary information
Supplementary Information
Supplementary Figures 15, Supplementary Notes 18 and Supplementary References. (PDF 186 kb)
Supplementary Movie 1
This animation demonstrates pattern formation in an oscillating chemical reaction on a multielectrode array. The image in the circle at Left illustrates the phase of each oscillator using color as well as a dot on the unit circle. The plots at Right show the control applied to the entire array (top) and the measured current traces (bottom). Initially, when no control is applied, there is no pattern formed. A control is then applied to form the pattern "O", in which phases of oscillators forming the letter are synchronized antiphase with respect to the others. The control is subsequently modified to form the pattern "K", and finally the initial control is repeated to revert to the pattern "O". (MOV 3792 kb)
Rights and permissions
This work is licensed under a Creative Commons Attribution 4.0 International License. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license, users will need to obtain permission from the license holder to reproduce the material. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/
About this article
Cite this article
Zlotnik, A., Nagao, R., Kiss, I. et al. Phaseselective entrainment of nonlinear oscillator ensembles. Nat Commun 7, 10788 (2016). https://doi.org/10.1038/ncomms10788
Received:
Accepted:
Published:
Further reading

Leveraging deep learning to control neural oscillators
Biological Cybernetics (2021)

Optimization of periodic input waveforms for global entrainment of weakly forced limitcycle oscillators
Nonlinear Dynamics (2021)

Phase modelbased neuron stabilization into arbitrary clusters
Journal of Computational Neuroscience (2018)

Towards Topological Mechanisms Underlying Experience Acquisition and Transmission in the Human Brain
Integrative Psychological and Behavioral Science (2017)
Comments
By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.