Phase-selective entrainment of nonlinear oscillator ensembles

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.

C omplex 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 man-made 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 ensemble-level properties while withstanding experimental uncertainties.
Control-theoretic 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 open-loop 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, open-loop 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 o j ¼ 2p/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 c j (refs 2,33). A value of c j ¼ 0 (equivalently 2p) 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 phase-selective 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 When the forcing frequency is near the natural frequencies in the ensemble, averaging theory 44 states that the mean difference j j between each phase c j and the forcing phase y ¼ Ot follows the time-independent dynamic equation where Do j ¼ o j À O is called the frequency detuning, and is a 2p-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 j Ã j that satisfies , then the phase of the jth oscillator becomes entrained to the forcing phase y with an average offset of j Ã j . 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(Ot) so that each system in equation (1) is entrained to a forcing frequency O, e.g., the mean of natural frequencies o 1 oo 2 oyoo N , and such that the jth oscillator cycles its orbit with a phase offset of j Ã j relative to the forcing phase y. The set of pairs ðo j ; j Ã j Þ constitutes a pattern for selective entrainment. We require that _ j j ¼ Do j þ L v ðj Ã j Þ¼ 0 eventually holds for each oscillator, that is, equation (2) has an attractive fixed point at j Ã j for all j at which the slope of the interaction function L v is negative 43 . The function that best satisfies these ideal conditions has steep decreases at phase values j Ã j where entrainment must occur, and crosses (from above) horizontal lines at frequency detuning values À Do j . This creates the desired attractive fixed points for equation (2). Because the interaction function is periodic, it must then increase so that L v (2p) ¼ L v (0) holds. Crossings of À Do 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 Dj Ã ¼j Ã 2 À j Ã 1 is assigned between two entrained oscillators. In Fig. 1c we desire in-phase synchronization at phase offsets of ðj Ã 1 ; j Ã 2 Þ¼ðp; pÞ, so the ideal interaction function has one steep decrease at j Ã 1 ¼j Ã 2 ¼p that intersects horizontal lines through À Do 1 (blue) and À Do 2 (red) at p radians. Figure 1d illustrates anti-phase synchronization with ðj Ã 1 ; j Ã 2 Þ¼ð0; pÞ, where the interaction function has two steep decreases at j Ã 1 ¼0 and j Ã 2 ¼p that intersect horizontal lines at À Do 1 (blue) and À Do 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 in-phase and anti-phase arrangements by the input waveform (shown above). These configurations are achieved regardless of initial oscillator phases, because the interaction function crosses the line À Do 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 2p region are achievable, with small deviations as Dj Ã approaches 2p, 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 j Ã j of equation (2) must be unique for each j, which is achieved when the interaction function crosses each horizontal line Do j from above only once at j Ã j . This is possible when the phase offsets are monotonically ordered as as demonstrated by segregation of 20 inhomogeneous electrochemical oscillators into clusters in the experiments described in Fig. 2. An anti-phase configuration ðj 1 ; j 2 Þ¼ð0; pÞ 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 two-cluster 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 , 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 À Do 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 j Ã j , 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 anti-phase 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 j Ã j ¼ 0 or j Ã j ¼ p, which correspond to 'on' (pattern) or 'off' (background) states, respectively, that are visualized using a colour scale where 0 ('on') is blue and p ('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 ) ¼  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 ðj Ã 1 ; j Ã 2 ; j Ã 3 ; j Ã 4 Þ¼ð0; 0; p; pÞ, which are achieved by an interaction function as in Fig. 2a. The phases ðj Ã 1 ; j Ã 2 ; j Ã 3 ; j Ã 4 Þ¼ð0; p; 0; pÞ used for 'K' are not monotonically ordered, so a precursor waveform is applied to generate globally attractive phase offsets j Ã 2 ¼p and j Ã 3 ¼0 for clusters 2 and 3. This anti-phase 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    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
Phase-selective 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(Ot) ¼ sin (Ot) 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 p, the maximum achievable distance between extremal phase offsets j Ã 1 and j Ã N is j Ã N À j Ã 1 op. Thus, a sinusoidal input cannot produce anti-phase 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 (non-zero) 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 phase-selective 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 set-up 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 self-sustained 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 B1.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 pulse-width of 0.05 s and post-processing 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 real-time data acquisition by a high-speed multifunction M Series DAQ PXI-6255 (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. The control waveform is a potential signal generated by the multifunction M Series DAQ PXI-6255 (National Instruments) and implemented in the WE by a potentiostat (Bank Instruments) as a superimposed signal on the applied baseline potential V 0 . Each nickel wire is connected in series to an individual resistance R ind ¼ 2,500 Ohm.  Fig. 4. Left: the best approximation to the desired interaction function; Right: one cycle of the resulting control waveform, obtained using the control design procedure described in section 'Control of pattern transitions in an ensemble', which is applied to entrain the electrochemical oscillators. The ensemble is grouped into four clusters of quantities (N 1 , N 2 , N 3 , N 4 ) ¼ (7, 7, 3, 3). Cluster 1 for oscillators j ¼ 1 to 7 (blue), cluster 2 for j ¼ 8 to 14 (cyan), cluster 3 for j ¼ 15 to 17 (yellow) and cluster 4 for j ¼ 18 to 20 (red). Phase assignment for: (a): Pattern 'O': ðj Ã 1 ; j Ã 2 ; j Ã 3 ; j Ã 4 Þ¼ð0; 0; p; pÞ. (b): Precursor of pattern 'K': ðj Ã 2 ; j Ã 3 Þ ¼ ðp; 0Þ. (c): Pattern 'K': ðj Ã 1 ; j Ã 2 ; j Ã 3 ; j Ã 4 Þ¼ð0; p; 0; pÞ.
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 o 1 ¼ 0.390 Hz (electrodes j ¼ 1 to 7 in blue), cluster 2 with o 2 ¼ 0.406 Hz (electrodes j ¼ 8 to 14 in cyan), cluster 3 with o 3 ¼ 0.427 Hz (electrodes j ¼ 15 to 17 in yellow) and cluster 4 with o 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.