PHASE-SELECTIVE ENTRAINMENT OF NONLINEAR OSCILLATOR ENSEMBLES

Abstract

A system for entraining an oscillator ensemble is disclosed that includes a plurality of oscillators in an entrained phase pattern. The system includes an entrainment device operatively coupled to each non-linear oscillator of the oscillator ensemble, and the entrainment control device is configured to deliver a 2π-periodic control signal v(θ) to all oscillators of the plurality of oscillators to induce the entrained phase pattern.

Claims

1-10. (canceled)

11. A system for entraining a non-linear oscillator ensemble in an entrained phase pattern, the system comprising: an entrainment control device operatively coupled to each of N non-linear oscillators of the non-linear oscillator ensemble, the entrainment control device is configured to deliver a 2π-periodic control signal v(θ) to each of the N non-linear oscillators, wherein each of the N non-linear oscillators endogenously oscillates at a natural frequency ω.sub.d with a corresponding natural phase φ.sub.j wherein j ranges from 1 to N; an entrainment pattern generator configured to determine a 2π-periodic interaction function Λ.sub.v comprising a relationship between an entrained frequency shift Δω and an entrained phase offset φ for the N non-linear oscillators, the 2π-periodic interaction function Λ.sub.v defining the entrained phase pattern comprising an entrained phase offset φ.sub.j for each of the N non-linear oscillators, each entrained phase offset φ.sub.j consisting of a difference between the natural phase ψ.sub.j and an entrained phase θ; an entrainment control signal generator configured to determine the 2π-periodic control signal v(θ) using the 2π-periodic interaction function Δ.sub.v and a phase response curve Z(θ), the phase response curve Z(θ) characterizing a phase shift of the N non-linear oscillators in response to a weak perturbation, wherein the 2π-periodic control signal v(θ) is configured to: detune the natural frequency ω.sub.j by the entrained frequency shift Δω.sub.j for each of the N non-linear oscillators, the entrained frequency shift Δω.sub.j consisting of a difference between the natural frequency ω.sub.j and an entrained frequency Ω; and shift each of the natural phases ψ.sub.j by the entrained phase offset φ.sub.j for each of the N non-linear oscillators.

12. The system of claim 11, wherein the N non-linear oscillators are detuned to the same entrained frequency Ω consisting of a mean of the natural frequencies ω.sub.j of the N non-linear oscillators.

13. The system of claim 11, wherein each entrained phase offset φ.sub.j is independently selected from a range of entrained phase offsets ranging between 0 and 2π.

14. The system of claim 11, wherein the 2π-periodic interaction function Λ.sub.v comprises a plurality of phase/detuning pairs (φ.sub.j*, Δω.sub.j), each phase/detuning pair (φ.sub.j*, Δω.sub.j) defining an average entrained phase offset φ.sub.j* and corresponding entrained frequency shift Δω.sub.j for each of the N non-linear oscillators in accordance with the entrained phase pattern, and wherein the 2π-periodic interaction function Λ.sub.v* further comprises a negative derivative $\frac{d}{d.Math..Math.\varphi }.Math.{\Lambda }_v<0$ at each phase/detuning pair (φ.sub.j*, Δω.sub.j).

15. The system of claim 14, wherein the 2π-periodic interaction function Λ.sub.v comprises an interaction function Λ.sub.v*, wherein: ${\Lambda }_v^{*}\left(\varphi \right)=r_1+\left(r_1-r_{N+1}\right).Math.\sigma \left(\frac{\varphi -s_1}{h_{N+1}}\right)+\underset{j=2}{\overset{N}{.Math.}}.Math.\left(r_{j+1}-r_j\right).Math.\sigma \left(\frac{\varphi -{\varphi }_j^{*}}{h_j}\right)+\underset{j=2}{\overset{N}{.Math.}}.Math.\left(r_{j+1}+r_j\right).Math.\sigma \left(\frac{\varphi +2.Math.\pi -{\varphi }_j^{*}}{h_j}\right)+\left(r_1-r_{N+1}\right).Math.\sigma \left(\frac{\varphi -s_{N+1}}{h_{N+1}}\right),.Math.r_1=-{\mathrm{\Delta \omega }}_1+\frac{1}{2}.Math.\left(-{\mathrm{\Delta \omega }}_{N-1}+{\mathrm{\Delta \omega }}_2\right),r_j=\frac{1}{2}.Math.\left(-{\mathrm{\Delta \omega }}_1+{\mathrm{\Delta \omega }}_{j-1}\right).Math..Math.\mathrm{for}.Math..Math.j=2,\mathrm{.Math.}.Math.,N,.Math..Math.r_{N+1}=-{\mathrm{\Delta \omega }}_N+\frac{1}{2}.Math.\left(-{\mathrm{\Delta \omega }}_{N-1}+{\mathrm{\Delta \omega }}_N\right),s_1=\frac{1}{2}.Math.\left({\varphi }_1^{*}+{\varphi }_N^{*}-2.Math.\pi \right),.Math..Math.s_j=\frac{1}{2}.Math.\left({\varphi }_j^{*}+{\varphi }_{j-1}^{*}\right).Math..Math.\mathrm{for}.Math..Math.j=2,\mathrm{.Math.}.Math.,N,s_{N+1}=\frac{1}{2}.Math.\left({\varphi }_1^{*}+{\varphi }_N^{*}+2.Math.\pi \right),.Math.h_1=\min .Math.\left\{s_2-{\varphi }_1^{*},\left({\varphi }_1^{*}-s_1\right)/2\right\},h_j=\min .Math.\left\{s_{j+1}-{\varphi }_j^{*},{\varphi }_j^{*}-s_j\right\},\mathrm{for}.Math..Math.j=2,\mathrm{.Math.}.Math.,N,.Math..Math.h_{N+1}=\min .Math.\left\{\frac{1}{2}.Math.\left(s_{N+1}-{\varphi }_N^{*}\right),\frac{1}{2}.Math.\left({\varphi }_1^{*}-s_1\right)\right\},.Math..Math.\sigma \left(x\right)=\frac{1}{2}.Math.\left(\mathrm{erf}\left(2.Math.x\right)+1\right),\mathrm{erf}\left(x\right)=\frac{2}{\sqrt{\pi }}.Math.{\int }_0^x.Math.e^{-t^2}.Math.\mathrm{dt}..Math.$

16. The system of claim 15, wherein the 2π-periodic control signal v(θ) is given by: $v\left(\theta \right)\approx v^r\left(\theta \right)=\frac{c_0}{2}+\underset{n=1}{\overset{r}{.Math.}}.Math.\left[c_n.Math.\cos \left(n.Math..Math.\theta \right)+d_n.Math.\sin \left(n.Math..Math.\theta \right)\right],$ wherein $c_0=2.Math.\frac{f_0}{a_0}.Math.X\left[a_0\ne 0\right],c_n=2.Math.\frac{f_n.Math.a_n+b_n.Math.g_n}{a_n^2+b_n^2}.Math.X\left[a_n^2+b_n^2\ne 0\right],d_n=2.Math.\frac{f_n.Math.b_n+a_n.Math.g_n}{a_n^2+b_n^2}.Math.X\left[a_n^2+b_n^2\ne 0\right],$ χ[A]=1 if A is true, and χ[A]=0 otherwise, a.sub.0, a.sub.n, and b.sub.n are obtained from a truncated Fourier-series expression of a phase response curve given by: $Z\left(\theta \right)\approx Z^r\left(\theta \right)=\frac{a_0}{2}+\underset{n=1}{\overset{r}{.Math.}}.Math.\left[a_n.Math.\cos \left(n.Math..Math.\theta \right)+b_n.Math.\sin \left(n.Math..Math.\theta \right)\right],$ and f.sub.0, f.sub.n, and g.sub.n are obtained from a truncated Fourier-series expression of an interaction function Λ.sub.v.sup.r given by: ${\Lambda }_v^r\left(\varphi \right)=\frac{f_0}{2}+\frac{1}{2}.Math.\underset{n=1}{\overset{r}{.Math.}}.Math.\left[f_n.Math.\cos \left(n.Math..Math.\varphi \right)+g_n.Math.\sin \left(n.Math..Math.\varphi \right)\right],$ wherein $f_0=\frac{a_0.Math.c_0}{2},f_n=a_n.Math.c_n+b_n.Math.d_n,\mathrm{and}.Math..Math.g_n=b_n.Math.c_n-a_n.Math.d_n.$

17. The system of claim 11, wherein the phase response curve Z(θ) consists of a mean of individual phase response curves Z.sub.j(θ), each individual phase response curve Z.sub.j(θ) corresponding to the j.sup.th non-linear oscillator of the N non-linear oscillators.

18. The system of claim 17, further comprising an oscillator characterization device operatively coupled to each of the N non-linear oscillators, the oscillator characterization device configured to determine, for each of the N non-linear oscillators, at least one of the natural frequency ω.sub.j, the natural phase ψ.sub.j, the phase response curve Z.sub.j(θ), and any combination thereof.

19. The system of claim 18, wherein the oscillator characterization device is configured to determine a phase response curve Z(ψ(t.sub.o)) for each of the N oscillators by: measuring a baseline phase ψ.sub.o(t.sub.o+NT) at N cycles after the absence of an applied pulse and a perturbed phase ψ.sub.1(t.sub.o+NT) at N cycles after the applied pulse, the applied pulse comprising an impulse with a duration Δt and a magnitude M; and determining Z(ψ(t.sub.o) according to: $Z\left(\psi \left(t_0\right)\right)=\frac{{\psi }_1\left(t_0+\mathrm{NT}\right)-{\psi }_0\left(t_0+\mathrm{NT}\right)}{M.Math..Math.\Delta .Math..Math.t}.$

20. The system of claim 11, wherein the system is configured to entrain the oscillator ensemble in a series of entrained phase patterns by sequentially determining a series of 2π-periodic interaction functions Λ.sub.v defining the series of entrained phase patterns, sequentially determining the series of 2π-periodic control signals v(θ) using the series of 2π-periodic interaction functions Λ.sub.v and the phase response curve Z(θ), and sequentially delivering the series of 2π-periodic control signals v(θ) to each of the N non-linear oscillators.

21. A method of entraining a non-linear oscillator ensemble in an entrained phase pattern, the method comprising: determining a 2π-periodic interaction function Λ.sub.v defining the entrained phase pattern for N non-linear oscillators of the non-linear oscillator ensemble, the entrained phase pattern comprising an entrained phase offset φ.sub.j relative to a natural phase ψ.sub.j for j=1, . . . , N, wherein the 2π-periodic interaction function Λ.sub.v comprises a relationship between an entrained frequency shift Δω and the entrained phase offset φ for the non-linear oscillator ensemble; determining the 2π-periodic control signal v(θ) using the 2π-periodic interaction function Λ.sub.v and a phase response curve Z(θ), the phase response curve Z(θ) characterizing a phase shift of the N non-linear oscillators in response to a weak perturbation; delivering the 2π-periodic control signal v(θ) to each of the N non-linear oscillators, wherein the 2π-periodic control signal v(θ) is configured to: detune the natural frequency ω.sub.j by the entrained frequency shift Δω.sub.j for each of the N non-linear oscillators, the entrained frequency shift Δω.sub.j consisting of a difference between the natural frequency ω.sub.j and an entrained frequency Ω; and shift each of the natural phases ψ.sub.j by the entrained phase offset φ.sub.j for each of the N non-linear oscillators, the entrained phase shift φ.sub.j consisting of a difference between the natural phase ψ.sub.j and an entrained phase θ;

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0009] The patent or application file contains at least one drawing executed in color. Copies of this patent or patent application publication with color drawing(s) will be provided by the Office upon request and payment of the necessary fee.

[0010] The figures described herein below illustrate various aspects of the disclosure.

[0011] FIG. 1 is a schematic illustration showing an electrochemical cell and associated circuitry used to demonstrate control of a phase structure of an ensemble of heterogeneous oscillators.

[0012] FIG. 2A is a graph illustrating an interaction function used to generate a control waveform to produce an “O” pattern in an ensemble of heterogeneous oscillators according to one aspect.

[0013] FIG. 2B is a graph illustrating a control waveform produced using the interaction function illustrated in FIG. 2A.

[0014] FIG. 2C is a graph illustrating an interaction function used to generate a control waveform to produce a precursor pattern in an ensemble of heterogeneous oscillators according to one aspect.

[0015] FIG. 2D is a graph illustrating a control waveform produced using the interaction function illustrated in FIG. 2C.

[0016] FIG. 2E is a graph illustrating an interaction function used to generate a control waveform to produce a “K” pattern in an ensemble of heterogeneous oscillators according to one aspect.

[0017] FIG. 2F is a graph illustrating a control waveform produced using the interaction function illustrated in FIG. 2E.

[0018] FIG. 3A is a graph showing phase response curves (PRCs) measured simultaneously for twenty working electrodes, and an inset graph showing observed current oscillations for the twenty working electrodes measured from the electrochemical cell and associated circuitry illustrated in FIG. 1.

[0019] FIG. 3B is graph showing the experimentally measured phase differences for the twenty working electrodes measured in FIG. 3A as a function of the corresponding desired phase differences.

[0020] FIG. 3C is a graph showing an ideal interaction function (dotted-line) and a respective best achievable approximation (solid-line) used to construct an in-phase assignment for two nonidentical electrode oscillators from the electrochemical cell and associated circuitry illustrated in FIG. 1.

[0021] FIG. 3D contains an upper graph showing a periodic control signal used to entrain the two nonidentical electrode oscillators in an in-phase pattern corresponding to the interaction function shown in FIG. 3C, and a lower graph showing a time-series of the output current measured from the two nonidentical electrode oscillators entrained in-phase: (φ.sub.1*, φ.sub.2*)=(π,π), with natural frequencies (ω.sub.1,ω.sub.2)=(0.330,0.348) Hz, and Ω=0.339 Hz with each oscillator response shown in blue and red; an inset graph shows the entrained in-phase oscillator phases on a unit circle.

[0022] FIG. 3E is a graph showing an ideal interaction function (dotted-line) and a respective best achievable approximation (solid-line) used to construct an anti-phase assignment for two nonidentical electrode oscillators from the electrochemical cell and associated circuitry illustrated in FIG. 1.

[0023] FIG. 3F contains an upper graph showing a periodic control signal used to entrain the two nonidentical electrode oscillators in an anti-phase pattern corresponding to the interaction function shown in FIG. 3E, and a lower graph showing a time-series of the output current measured from the two nonidentical electrode oscillators entrained anti-phase: (φ.sub.1*, φ.sub.2*)=(0,π), with natural frequencies (ω.sub.1,ω.sub.2)=(0.443,0.480) Hz, and Ω=0.462 Hz with each oscillator response shown in blue and red; an inset graph shows the entrained in-phase oscillator phases on a unit circle.

[0024] FIG. 4A is a graph showing an interaction function and the intersections of each phase offset φ.sub.j* with the interaction function for a phase assignment of an oscillator ensemble (N.sub.1.N.sub.2)=(10. 10) including cluster 1 (blue) with (φ.sub.1=0 and ω.sub.1=0.450 Hz as well as cluster 2 (red) with φ.sub.2=π and φ.sub.2=0.471 Hz. both clusters were forced at Ω=0.463 Hz.

[0025] FIG. 4B includes a lower graph showing an entrained time-series corresponding to the phase assignment illustrated in FIG. 4A, an upper graph showing a periodic waveform of a control signal used to produce the phase assignment illustrated in FIG. 4A, and an inset graph showing the phase offsets projected onto a unit circle.

[0026] FIG. 4C is a graph showing an interaction function and the intersections of each phase offset φ.sub.j* with the interaction function for a phase assignment of an oscillator ensemble (N.sub.1.N.sub.2)=(1. 19) including electrode 1 (blue) with φ.sub.1=0 and φ.sub.1=0.419 Hz as well as cluster 2 (red) with ω.sub.2=π and ω.sub.2=0.454 Hz. both clusters were forced at Ω=0.450 Hz.

[0027] FIG. 4D includes a lower graph showing an entrained time-series corresponding to the phase assignment illustrated in FIG. 4C, an upper graph showing a periodic waveform of a control signal used to produce the phase assignment illustrated in FIG. 4C, and an inset graph showing the phase offsets projected onto a unit circle.

[0028] FIG. 4E is a graph showing an interaction function and the intersections of each phase offset φ.sub.j* with the interaction function for a phase assignment of an oscillator ensemble (N.sub.1.N.sub.2.N.sub.3.N.sub.4)=(4.4.4.4) including cluster 1 (blue) with φ.sub.1=0 and ω.sub.1=0.386 Hz, cluster 2 (red) with φ.sub.2=1.1 rad and ω.sub.2=0.404 Hz, cluster 3 (green) with φ.sub.3=2.1 rad and ω.sub.3=0.421 Hz, and cluster 4 (yellow) with φ.sub.4=π and ω.sub.4=0.440 Hz. all clusters were forced at Ω=0.413 Hz.

[0029] FIG. 4F includes a lower graph showing an entrained time-series corresponding to the phase assignment illustrated in FIG. 4E, an upper graph showing a periodic waveform of a control signal used to produce the phase assignment illustrated in FIG. 4C, and an inset graph showing the phase offsets projected onto a unit circle.

[0030] FIG. 5A is a color map showing the natural frequencies for four spatially distributed oscillator clusters associated with the electrochemical cell and associated circuitry shown in FIG. 1, in which the clusters have a size distribution of (N.sub.1.N.sub.2.N.sub.3.N.sub.4)=(7.7.3.3) as arranged by natural frequency in ascending order (ω.sub.1<ω.sub.2<ω.sub.3<ω.sub.4). Cluster 1 contains electrodes 1-7 (blue), cluster 2 contains electrodes 8-14 (cyan), cluster 3 contains electrodes 15-17 (yellow), and cluster 4 contains electrodes 18-20 (red).

[0031] FIG. 5B is a color map showing a snapshot phase assignment of the four spatially distributed oscillator clusters shown in FIG. 5A at a forcing phase θ=0 and corresponding to an “O” pattern, in which (φ.sub.1*,φ.sub.2*,φ.sub.3*,φ.sub.4*)=(0,0, π, π).

[0032] FIG. 5C is a color map showing a snapshot phase assignment of the four spatially distributed oscillator clusters shown in FIG. 5A at a forcing phase θ=0 and corresponding to a precursor of a “K” pattern, in which (φ.sub.2*,φ.sub.3*)=(π,0).

[0033] FIG. 5D is a color map showing a snapshot phase assignment of the four spatially distributed oscillator clusters shown in FIG. 5A at a forcing phase θ=0 and corresponding to a “K” pattern, in which (φ.sub.1*,φ.sub.2*,φ.sub.3*,φ.sub.4*)=(0,π,0,π).

[0034] FIG. 5E is a graph showing the “O” phase assignment pattern shown in FIG. 5B projected onto a unit circle.

[0035] FIG. 5F is a graph showing the precursor of the “K” pattern phase assignment pattern shown in FIG. 5C projected onto a unit circle.

[0036] FIG. 5G is a graph showing the “K” pattern phase assignment pattern shown in FIG. 5D projected onto a unit circle.

[0037] FIG. 5H is a graph showing an interaction function corresponding to the “O” phase assignment pattern shown in FIG. 5B.

[0038] FIG. 5I is a graph showing an interaction function corresponding to the precursor of the “K” pattern phase assignment pattern shown in FIG. 5C.

[0039] FIG. 5J is a graph showing an interaction function corresponding to the “K” pattern phase assignment pattern shown in FIG. 5D.

[0040] FIG. 6A contains a lower graph showing the current oscillations measured from the electrochemical cell and associated circuitry shown in FIG. 1, and an upper graph showing the control signal used to obtain the current oscillations; no control signal was applied to the four electrode clusters (N.sub.1;N.sub.2;N.sub.3;N.sub.4)=(7;7;3;3) with natural frequencies of (ω.sub.1,ω.sub.2,ω.sub.3,ω.sub.4)=(0.390,0.406,0.427,0:442) Hz.

[0041] FIG. 6B contains a graph showing the phase pattern corresponding to the measured current oscillations shown in FIG. 6A at time t.sub.s=92 sec, as projected on to a unit circle.

[0042] FIG. 6C contains a graph showing a map of the phase pattern corresponding to the measured current oscillations shown in FIG. 6A at time t.sub.s=92 sec.

[0043] FIG. 6D contains a lower graph showing the current oscillations measured from the electrochemical cell and associated circuitry shown in FIG. 1, and an upper graph showing the forcing control signal at Ω=0.408 Hz used to obtain the current oscillations shown in the lower graph corresponding to an “O” phase distribution pattern characterized by (φ.sub.1*,φ.sub.2*,φ.sub.3*,φ.sub.4*)=(0,0, π, π).

[0044] FIG. 6E contains a graph showing the phase assignment corresponding to the “O” phase distribution pattern shown in FIG. 6D at time t.sub.s=351 sec, as projected on to a unit circle.

[0045] FIG. 6F contains a graph showing a map of the “O” phase distribution pattern corresponding to the measured current oscillations shown in FIG. 6D at time t.sub.s=351 sec.

[0046] FIG. 6G contains a lower graph showing the current oscillations measured from the electrochemical cell and associated circuitry shown in FIG. 1, and an upper graph showing the forcing control signal at Ω=0.408 Hz used to obtain the current oscillations shown in the lower graph corresponding to a precursor to a “K” phase distribution pattern characterized by (φ.sub.2*,φ.sub.3*)=(π,0).

[0047] FIG. 6H contains a graph showing the phase assignment corresponding to the precursor to a “K” phase distribution pattern shown in FIG. 6G at time t.sub.s=570 sec, as projected on to a unit circle.

[0048] FIG. 6I contains a graph showing a map of the precursor to the “K” phase distribution pattern shown in FIG. 6G at time t.sub.s=570 sec.

[0049] FIG. 6J contains a lower graph showing the current oscillations measured from the electrochemical cell and associated circuitry shown in FIG. 1, and an upper graph showing the forcing control signal at Ω=0.408 Hz used to obtain the current oscillations shown in the lower graph corresponding to a “K” phase distribution pattern characterized by (φ.sub.1*,φ.sub.2*,φ.sub.3*,φ.sub.4*)=(0, π, 0, π).

[0050] FIG. 6K contains a graph showing the phase assignment corresponding to the “K” phase distribution pattern shown in FIG. 6J at time t.sub.s=721 sec, as projected on to a unit circle.

[0051] FIG. 6L contains a graph showing a map of the “K” phase distribution pattern shown in FIG. 6J at time t.sub.s=721 sec.

[0052] FIG. 6M contains a lower graph showing the current oscillations measured from the electrochemical cell and associated circuitry shown in FIG. 1, and an upper graph showing the forcing control signal at Ω=0.408 Hz used to obtain the current oscillations shown in the lower graph corresponding to a second “O” phase distribution pattern characterized by (φ.sub.1*,φ.sub.2*,φ.sub.3*,φ.sub.4*)=(0,0, π, π).

[0053] FIG. 6N contains a graph showing the phase assignment corresponding to the second “O” phase distribution pattern shown in FIG. 6M at time t.sub.s=1001 sec, as projected on to a unit circle.

[0054] FIG. 6O contains a graph showing a map of the second “O” phase distribution pattern corresponding to the measured current oscillations shown in FIG. 6D at time t.sub.s=1001 sec.

[0055] FIG. 7 is a graph showing several properties of an interaction function Λ.sub.v (φ); the maximum and minimum values Λ.sub.v (φ.sup.+) and Λ.sub.v (φ.sup.−), which occur at the phases φ.sup.+ and φ.sup.−, respectively, determine the range of frequency detuning for which the oscillator can be entrained using weak forcing. The solutions φ.sub.1.sup.∞ and φ.sub.2.sup.∞ of the equation Λ.sub.v(φ)=−Δω that satisfy Λ′(φ)<0 determine the average phase shift, relative to Ωt, at which the oscillation stabilizes from a given initial phase. Initial phases in the regions (φ.sub.1.sup.∞) (grey) and (φ.sub.2.sup.∞) (green) result in asymptotic phase shifts of φ.sub.1.sup.∞ and φ.sub.2.sup.∞, respectively. When Λ.sub.v(φ)>−Δω then {dot over (φ)}>0 so the phase increases, and when Λ.sub.v(φ)<−Δω the {dot over (φ)}<0 so the phase decreases, as indicated by the arrows to illustrate equation (10);

[0056] FIG. 8 is a graph showing an interaction function Λ.sub.v (φ) that satisfies the conditions (13)-(14) with φ.sub.11.sup.∞ for j=1, φ.sub.12.sup.∞ or φ.sub.22.sup.∞ for j=2, φ.sub.13.sup.∞ or φ.sub.23.sup.∞ for j=3, and yon for j=4. The attractive regions /(φ.sub.ij.sup.∞) are indicated in grey for i=1 and green for i=2.

[0057] FIG. 9 is a graph showing an interaction function Λ.sub.v (φ) designed using the monotonicity method that satisfies the conditions in equations (13x)-(14x) with a globally attractive phase pattern. The unique fixed points for each oscillator are (φ11.sup.∞, φ.sub.2.sup.∞, φ.sub.33.sup.∞, and φ.sub.41.sup.∞ for j=1; 2; 3; 4, respectively. The attractive regions are A.sub.j(φ.sub.1j.sup.∞)=[0,2π) for all j, as indicated.

[0058] FIG. 10 is a graph showing an interaction function Λ.sub.v (φ) corresponding to a precursor input that entrains oscillators j=2 and 3 to the desired phase offsets φ.sub.2*=φ.sub.22.sup.∞ and φ.sub.3*=φ.sub.33.sup.∞. The attractive regions A.sub.2 (φ.sub.22.sup.∞)=[0,2π) and A.sub.3(φ.sub.12.sup.∞)=[0,2π) are global. The waveform generated from the interaction function in FIG. 8 repeated here as an orange dashed line, can be subsequently applied to finalize the pattern.

[0059] FIG. 11 is a graph showing an interaction function Λ.sub.v (φ) corresponding to desynchronization of a large ensemble of oscillators with frequencies distributed uniformly on the interval [ω.sub.min,ω.sub.max]. The frequency detuning values are distributed uniformly on [−ω.sub.max,ω.sub.min], so that the control generated from Λ.sub.v (φ) will disperse the phase offsets uniformly on the interval [φ.sub.min.sup.∞, φ.sub.max.sup.∞] The asymptotic offset φ.sub.j.sup.∞ for an oscillator with an intermediate frequency ω.sub.jε[ω.sub.min,ω.sub.max] is indicated.

[0060] FIG. 12A is a graph showing an interaction function Λ.sub.v.sup.r(φ) designed using an r=5 order Fourier series fit to an ideal Λ.sub.v(φ) designed using the monotone step method. Inset: PRC of Hodgkin-Huxley (HH) System.

[0061] FIG. 12B is a graph showing an input waveform v(θ) resulting from Fourier inversion of Λ.sub.v.sup.r(φ).

[0062] FIG. 12C is a graph showing an ideal, designed and actual phase pattern resulting from applying v(θ) to a three uncoupled HH phase models with frequencies at 1.01, 1, and 0.99 times the natural frequency ω.

[0063] FIG. 13A is a graph showing an interaction function Λ.sub.v.sup.r(φ) designed using a r=5 order Fourier series fit to an ideal Λ.sub.v(φ) designed using the monotone step method. Inset: PRC of Hodgkin-Huxley (HH) System.

[0064] FIG. 13B is a graph showing an input waveform v(θ) resulting from Fourier inversion of Λ.sub.v.sup.r(φ).

[0065] FIG. 13C is a graph showing an ideal, designed and actual phase pattern resulting from applying v(θ) to a network of three HH neurons with frequencies at 1.01, 1, and 0.99 times the natural frequency ω.

[0066] FIG. 14A is a graph showing an interaction function Λ.sub.v.sup.r(φ) designed using a Fourier series fit to an ideal Λ.sub.v(φ) designed using the monotone step method. Inset: PRC of Morris-Lecar (ML) System.

[0067] FIG. 14B is a graph showing an input waveform v(θ) resulting from Fourier inversion of Λ.sub.v.sup.r(φ).

[0068] FIG. 14C is a graph showing an ideal and an actual phase pattern resulting from applying v(θ) to a network of five weakly coupled ML neurons with frequencies at 0.98, 0.99, 1, 1.01, 1.02 times the natural frequency ω.

[0069] FIG. 15A is a schematic diagram showing a configuration of a network of 15 chemical oscillators.

[0070] FIG. 15B is a graph showing a distribution of frequencies of the oscillators.

[0071] FIG. 15C is a graph showing a designed input waveform v(θ).

[0072] FIG. 15D is a graph showing an ideal and an actual phase pattern resulting from applying v(θ).

[0073] FIG. 16 illustrates an example configuration of a remote device system.

[0074] FIG. 17 illustrates an example configuration of a server system.

DETAILED DESCRIPTION OF THE INVENTION

[0075] A methodology is presented 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. A 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. The theoretical methodology is demonstrated in practice with experiments to control complex electrochemical ensembles whose dynamics are variable and unknown, and for which state information is unavailable.

[0076] External stimuli that alter neuronal spiking activities are designed. Practically applicable techniques have been derived to control dynamic structures, such as synchrony or clustering, in a population of neurons in the brain.

[0077] The developed method can, among other things, be (1) used for the advancement of the neurostimulation technology for the treatment of neurological disorders, such as Parkinson's disease and epilepsy; (2) adopted to design neurocomputing architectures.

A. Phase Model Approximation

[0078] The effect of inputs on periodic dynamical systems may be approximated using phase models which can be computed for systems with known dynamics or experimentally inferred in practice when the dynamics are unknown. These models are used to characterize circadian cycles, cardiac rhythms, and phenomena in neural and chemical systems, and their simplicity has enabled control design for neuron models 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. 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 neighbors. However, for the manipulation and desynchronization of electrochemical and neural systems, 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, i.e., asymptotic phase structure, are manipulated rather than steering the system states directly.

B. Entrainment of an Ensemble

[0079] 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.sub.j, and is represented by the Winfree phase model

{dot over (ψ)}.sub.j=ω.sub.j+Z(ψ.sub.j)u  (1)

where φ.sub.j=2π/T.sub.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 ψ.sub.j. A value of ψ.sub.j=0 (equivalently 2π) corresponds to a measurement of the j.sup.th system reaching its maximum.

[0080] i) Phase Coordinate Transformation

[0081] The technique for phase coordinate transformation stems from the early work of Andronov and Vitt, and can be derived directly from the method of Malkin. An oscillator is formally defined as a system of smooth ordinary differential equations

{dot over (x)}=f(x,u),x(0)=x.sub.0,tε[0,∞)  (1x) [0082] where x(t)ε.sup.n is the state and u(t)εis a control function. The oscillatory property of equation (1x) implies the existence of an attractive, non-constant limit cycle γ(t)=γ(t+T) which is a solution to the unperturbed system {dot over (γ)}=f(γ,0), on the periodic orbit defined by Γ={γε.sup.n: y=γ(t) for 0≦t<T}⊂.sup.n. The key idea is to define a mapping between γ(t) to ψ(t) in a way such that the phase of the unperturbed system advances proportionally to time along the periodic orbit, i.e., ψ(t)=ωt, where ω=2π/T is the natural frequency of oscillation for the unforced oscillator. The dynamics equation (1) can then be approximated in phase coordinates by a scalar equation

{dot over (ψ)}=ω+Z(ψ)u,  (2x)

which is called a phase model, where Z is the phase response curve (PRC), also called the infinitesimal PRC or iPRC, which quantifies the asymptotic phase shift due to an infinitesimal perturbative input applied at a given phase on the limit cycle. The phase variable ψ(t) is associated to the isochronous manifold corresponding to γ(t)εΓ, from which all undisturbed trajectories asymptotically synchronize. ψ(t)=0 (mod2π) when an observed variable in the state vector y=γ(t) attains its maximum on Γ. The phase model equation (2x) accurately represents the dynamics equation (1x) for sufficiently weak inputs u(t) such that the solution x(t, x.sub.0, u) to equation (1x) remains within a small enough neighborhood of Γ. The PRC for a system with known dynamics of the form in equation (1x) can be computed by various procedures, including an approach similar to that of Malkin that is described in Appendix C of Malkin.

[0083] Phase selective entrainment techniques are developed using experiments in which nickel electrodes undergo anodic dissolution in sulfuric acid and exhibit electrochemical oscillations, for which the experimental apparatus is described in Section IV A. Prior to these experiments, the PRCs of the ensemble elements, which are shown FIG. 2A, 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, described below.

[0084] 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 states that the mean difference φ.sub.j between each phase ψ.sub.j and the forcing phase θ=Ωt follows the time-independent dynamic equation

{dot over (φ)}.sub.j=Δω.sub.j+Λ.sub.v(φ.sub.j)  (2)

where Δω.sub.j=ω.sub.j−Ω is called the frequency detuning, and

[00001]$\begin{array}{cc}{\Lambda }_v\left(\varphi \right)=\frac{1}{2.Math.\pi }.Math.{\int }_0^{2.Math.\pi }.Math.Z\left(\theta +\varphi \right).Math.v\left(\theta \right).Math.d.Math.\theta & \left(3\right)\end{array}$

is a 2π-periodic interaction function that characterizes the average effect of the periodic input on the oscillation phase. Such ergodic averaging is discussed in more detail herein below and illustrated in FIG. 7. If equation (2) has a unique attractive fixed point φ.sub.j* that satisfies Δω.sub.j+Λ.sub.v(φ.sub.j*)=0, then the phase of the j.sup.th oscillator becomes entrained to the forcing phase θ with an average offset of φ.sub.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. This approach is reversed by choosing the desired asymptotic behavior, constructing a suitable interaction function, and using the PRC to obtain the input by circular de-convolution of equation (3).

C. Interaction Functions for Phase Selection

[0085] The input v(Ωt) is designed so each system in equation (1) is entrained to a forcing frequency Ω, e.g., the mean of natural frequencies ω.sub.1<ω.sub.2< . . . <ω.sub.N, and such that the j.sup.th oscillator cycles its orbit with a phase offset of φ.sub.j*; relative to the forcing phase θ. The set of pairs (ω.sub.j,φ.sub.j*) constitutes a pattern for selective entrainment. It is required that {dot over (φ)}.sub.j=Δω.sub.j+Λ.sub.v(φ.sub.j*)=0 eventually holds for each oscillator, i.e., equation (2) has an attractive fixed point at φ.sub.j* for all j at which the slope of the interaction function Λ.sub.v is negative. The function that best satisfies these ideal conditions has steep decreases at phase values φ.sub.j* where entrainment must occur, and crosses (from above) horizontal lines at frequency detuning values−Δω.sub.j. This creates the desired attractive fixed points for equation (2). Because the interaction function is periodic, it must then increase so that Λ.sub.v(2π)=Λ.sub.v(0) holds. Crossings of −Δω.sub.j from below are unstable fixed points, and do not affect convergence.

[0086] The concept is illustrated in FIG. 6A, FIG. 6B, FIG. 6C, FIG. 6D, FIG. 6E, FIG. 6F, FIG. 6G, FIG. 6H, FIG. 6I, FIG. 6J, FIG. 6K, FIG. 6L, FIG. 6M, FIG. 6N, and FIG. 6O, which describe an experiment where a phase difference Δφ*=φ.sub.2*−φ.sub.1* is assigned between two entrained oscillators. In-phase synchronization at phase offsets of (φ.sub.1*,φ.sub.2*)=(π,π) were desired, so the ideal interaction function has one steep decrease at φ.sub.1*=φ.sub.2*=π that intersects horizontal lines through −Δω.sub.1 (blue) and −Δω.sub.2 (red) at it radians. FIG. 6D illustrates anti-phase synchronization with (φ.sub.1*,φ.sub.2*)=(0,π), where the interaction function has two steep decreases at (φ.sub.1*=0 and φ.sub.2*=π that intersect horizontal lines at −Δω.sub.1 (blue) and −Δω.sub.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. 3D and FIG. 3F 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 −Δω.sub.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. 3C.

D. Separation of Oscillator Ensembles into Phase Clusters

[0087] Uniquely attractive phase patterns are desired, where a common input synchronizes the oscillators to a pattern independently of their initial phases. The fixed point φ.sub.j* of equation (2) must be unique for each j, which is achieved when the interaction function crosses each horizontal line Δω.sub.j from above only once at φ.sub.j*. This is possible when the phase offsets are monotonically ordered as φ.sub.1*<φ.sub.2*< . . . <φ.sub.N* , for ω.sub.1<ω.sub.2< . . . <ω.sub.N, as demonstrated by segregation of twenty inhomogeneous electrochemical oscillators into clusters in the experiments described in FIG. 4A, FIG. 4B, FIG. 4C, FIG. 4D, FIG. 4E, and FIG. 4F. An anti-phase configuration (φ.sub.1, φ.sub.1)=(0,π) is achieved for electrodes in balanced (N.sub.1, N.sub.2)=(10,10) and unbalanced (N.sub.1, N.sub.2)=(1,19) clusters in FIG. 4A and FIG. 4C, 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. FIG. 4E shows the formation of four balanced clusters of (N.sub.1, N.sub.2, N.sub.3, N.sub.4)=(5, 5, 5, 5) oscillators with the phase structure (φ.sub.1*, φ.sub.2*, φ.sub.3*, φ.sub.4*)=(0, 1.1, 2.1, π) radians. The phase offsets φ.sub.j* are increasing as −Δω.sub.j decreases (and φ.sub.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. 2A are very similar, yet the oscillator frequencies are sufficiently heterogeneous for the technique to create the patterns in FIG. 4A, FIG. 4B, FIG. 4C, FIG. 4D, FIG. 4E, and FIG. 4F.

E. Control of Pattern Transitions in an Ensemble

[0088] Additionally, patterns were established without monotone phase ordering by designing an interaction function of the form at the bottom of FIG. 5I, which crosses a horizontal line at −Δω.sub.j from above twice, yielding multiple possible entrained phases and dependence on initial conditions. Precursor waveforms are applied to steer subsets of the ensemble into attractive regions for the desired phase offsets φ.sub.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. 5A, FIG. 5B, FIG. 5C, FIG. 5D, FIG. 5E, FIG. 5F, FIG. 5G, FIG. 5H, FIG. 5I, and FIG. 5J 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. 4A, FIG. 4B, FIG. 4C, FIG. 4D, FIG. 4E, and FIG. 4F. Anti-phase assignment between clusters are produced to display the letter “O”, then switch the input to produce the letter “K”. Rhythmic elements are assigned desired phase offsets of φ.sub.j*=0 or φ.sub.j*=π, which correspond to “on” (pattern) or “off” (background) states, respectively, that are visualized using a color scale where 0 (“on”) is blue and it (“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.sub.1, N.sub.2, N.sub.3, N.sub.4)=(7, 7, 3, 3) elements with mean natural frequencies (ω.sub.1, φ.sub.2, ω.sub.3, ω.sub.4)=(0.390, 0.406, 0.427, 0.442) Hz are positioned in the spatial arrangement in FIG. 5A. FIG. 5B, FIG. 5C, and FIG. 5D show the spatial distribution of phase offsets, FIG. 5E, FIG. 5F, and FIG. 5G show the corresponding structures on the unit circle, and FIG. 5H, FIG. 5I, and FIG. 5J show corresponding sketches of ideal interaction functions.

[0089] The pattern “O”, shown in FIG. 5B, is realized using the phases (φ.sub.1*,φ.sub.2*,φ.sub.3*,φ.sub.4*)=(0,0,π,π), which are achieved by an interaction function as in FIG. 5H. The phases (φ.sub.1*,φ.sub.2*,φ.sub.3*,φ.sub.4*)=(0,π,0,π) used for “K” are not monotonically ordered, so a precursor waveform is applied to generate globally attractive phase offsets φ.sub.2*=π and φ.sub.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. 5F. FIG. 5J 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 FIG. 5G, 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. 5H. Additional descriptions are provided of interaction function construction, pattern realizability, design for monotonically ordered patterns, and precursor waveform engineering, described herein below.

[0090] Measurements of the “O.fwdarw.K.fwdarw.O” pattern switching experiment appear in FIG. 6A, FIG. 6B, FIG. 6C, FIG. 6D, FIG. 6E, FIG. 6F, FIG. 6G, FIG. 6H, FIG. 6I, FIG. 6J, FIG. 6K, FIG. 6L, FIG. 6M, FIG. 6N, and FIG. 6O. FIG. 6A 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. 6G, FIG. 6J and FIG. 6M, and correspond to FIG. 5B, FIG. 5C and FIG. 5D, respectively. The ensemble is returned to pattern “O”, as shown in FIG. 6I, to demonstrate switching. The automatically generated interaction functions and control waveforms are herein below.

III. Discussion

[0091] 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.

[0092] Such resilience is required for nondestructive control of under-actuated, noisy, and uncertain biological and chemical ensembles that cannot be readily observed. For example, an effective technology for neurological treatment of Parkinson's disease 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. 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 and sleep irregularities. The formalism could also represent the entrainment that occurs in circadian timekeeping.

[0093] 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 φ.sub.1* and φ.sub.N*, is |φ.sub.N*−φ.sub.1*|<π. Thus, a sinusoidal input cannot produce anti-phase synchronization. This approach enables more versatile manipulation of phase relationships beyond this limitation. The application to desynchronization is described herein below and FIG. 11, and quantify how this 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.

[0094] The disclosed methodology takes advantage of approximations that are possible in the specific experimental setting. In the experiments described herein, 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. The methodology is expected 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 function based phase description is possible.

[0095] 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 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.

ii. Model Identification for Oscillating Systems

[0096] Several approaches to model identification have been used to estimate the PRC in real oscillating systems subject to noise and disturbances. A technique is described for phase model identification, which has been applied previously in experiments on electrochemical systems. This model identification procedure does not require online state observations or feedback, but rather relies on post-processing of a pseudo-random input sequence and observations of the effect on a state variable. Suppose that a brief impulse of duration Δt and magnitude M is applied to an oscillator at time t.sub.0, when the phase is ψ(t.sub.0). Specifically, suppose u(t)=M for tε[t.sub.0,t.sub.0+Δt], and u(t)=0 elsewhere. Let ψ.sub.1(t.sub.0+NT) and ψ.sub.0(t.sub.0+NT) represent the phase value when the oscillator has completed N cycles after the pulse is applied, and in the absence of a pulse, respectively. Because Γ is strongly attractive, the system will relax back to the periodic orbit several cycles after the pulse is applied. Assuming that the pulse duration Δt is brief, Z(ψ(t)) is approximately constant on tε[t.sub.0,t.sub.0+Δt], so that integrating equation (2x) results in

ψ.sub.0(t.sub.0+NT)=ωNT+ψ(t.sub.0),  (3x)

ψ.sub.1(t.sub.0+NT)=ωNT+ψ(t.sub.0)+∫.sub.t.sub.0.sup.t.sup.0.sup.+ΔtZ(ψ(t)Mdt≈ωNT+ψ(t.sub.0)+Z(ψ(t.sub.0)MΔt  (4x)

[0097] Subtracting equation (3) from equation (4) and solving for Z(ψ(t.sub.0)) results in

[00002]$\begin{array}{cc}Z(\psi \left(t_0\right)=\frac{{\psi }_1\left(t_0+\mathrm{NT}\right)-{\psi }_0\left(t_0+\mathrm{NT}\right)}{M.Math..Math.\Delta .Math..Math.t}& \left(5.Math.x\right)\end{array}$

[0098] The phases ψ.sub.0 and ψ.sub.1 are valued at 0 (equivalently 2π) at peaks of the observed state variable, and their values between peaks are reasonably approximated by linear interpolation. If y(t) is the observed oscillator output, let T.sub.−1 and T.sub.N denote the latest time at which y(t) reaches a peak before time t.sub.0 and the time when y(t) reaches the N.sup.th peak after the pulse is applied. Because phase on the limit cycle is linearly proportional to time, i.e., ψ(t)=ωt, the phases at two distinct time points t.sub.1, t.sub.2ε[0, T) will satisfy ψ(t.sub.1)−ψ(t.sub.2)=ω.Math.(t.sub.1−t.sub.2). As a result, the linear interpolation yields:

[00003]$\begin{array}{cc}\begin{array}{c}{\psi }_1\left(t_0+\mathrm{NT}\right)-{\psi }_0\left(t_0+\mathrm{NT}\right)\approx .Math.\omega .Math.\left(T_N-T_{-1}\right)-\omega .Math.\mathrm{NT}\\ =.Math.2.Math.\pi .Math.\frac{\left(T_N-T_{-1}\right)-\mathrm{NT}}{T},\end{array}& \left(6.Math.x\right)\end{array}$

where N is greater than the sufficient number of limit cycles for the perturbed system ψ.sub.1 to have relaxed to the periodic orbit Γ. In one implementation, a sequence of pulses is created at intervals of T(N+r.sub.j), where r.sub.jε(0,1) are pseudo-random, so that post-processing results in a time series Z(ψ.sub.j)ψ.sub.j=ψ.sub.j−1+ωT(N+r.sub.j)(mod2π). The pulse sequence is applied, and the measured peak times are used to estimate the phase differences in equation (6x) in a postprocessing step. The data are fitted using a Fourier series to estimate the PRC Z(ψ). One important advantage of this approach is to enable the concurrent identification of the PRCs for a large ensemble of rhythmic systems with simultaneously measurable output by using a single globally applied pulse sequence applied to the entire ensemble. The measurement state variable signal for each oscillator can then subsequently be processed off-line to produce a collection of phase models. This procedure is used to estimate the PRCs of each electrochemical oscillator in an experimental array, which are shown in FIG. 2A of the main text.
iii) Ergodic Averaging for Oscillating Systems

[0099] The central mechanism of the disclosed control methodology is entrainment, which refers to the dynamic synchronization of an oscillating system to a periodic input. The properties that characterize entrainment of the system in equation (2x) to a periodic forcing control u(t)=v(Ωt), where v is 2π-periodic and Ω is the forcing frequency, are closely approximated by ergodic averaging when the forcing signal is weak, i.e., v=εv.sub.1 where v.sub.1 has unit energy and ε<<1. Given such an input, the system in equation (1x) is guaranteed to remain in a neighborhood of Γ in which the phase model in equation (2x) remains valid. A slow phase variable is defined by φ(t)=ψ(t)−Ωt, and the difference Δω=ω−Ω is called the frequency detuning. The dynamic equation for the slow phase is then

{dot over (φ)}={dot over (ψ)}−Ω=Δω+Z(Ωt+φ)v(Ωt),  (7x)

where {dot over (φ)} is called the phase drift. The asymptotic behavior of equation (7x) is obtained by ergodic averaging, which eliminates the explicit dependence on time on the right hand side. An averaging operator :.fwdarw. is defined where is the set of 2π-periodic functions on , by

[00004]$\begin{array}{cc}〈x〉=\frac{1}{2.Math.\pi }.Math.{\int }_0^{2.Math.\pi }.Math.x\left(\theta \right).Math.d.Math..Math.\theta .& \left(8.Math.x\right)\end{array}$

[0100] Then the weak ergodic theorem for measure-preserving dynamical systems on the torus implies that for any periodic function v, the interaction function

[00005]$\begin{array}{cc}\begin{array}{c}\Lambda \left(\phi \right).Math.\stackrel{△}{=}.Math..Math.〈Z\left(\theta +\phi \right).Math.v\left(\theta \right)〉\\ =.Math.\frac{1}{2.Math.\pi }.Math.{\int }_0^{2.Math.\pi }.Math.Z\left(\theta +\phi \right).Math.v\left(\theta \right).Math.d.Math..Math.\theta \\ =.Math.{\lim }_{T.fwdarw.\infty }.Math.\frac{1}{T}.Math.{\int }_0^T.Math.Z\left(\Omega .Math..Math.t+\phi \right).Math.v\left(\Omega .Math..Math.t\right).Math.\mathrm{dt}\end{array}& \left(9.Math.x\right)\end{array}$

exists as a continuous, 2π-periodic function in where θ=Ωt is called the forcing phase.

[0101] The formal averaging theorem permits us to approximate equation (7x) by the averaged system

{dot over (φ)}=Δω+Λ.sub.v(φ),  (10x)

in the sense that there exists a change of variables φ=φ+εh(φ,φ) that maps solutions of equation (7x) to those of equation (10x). Without being limited to any particular theory, it has been proved previously that this approximation is of order (ε.sup.2), and hence it is appropriate in the case of weak forcing (v=εv.sub.1 with ε<<1). The averaged equation (10x) is autonomous, and approximately characterizes the asymptotic behavior of the system in equation (2x) under periodic forcing. Specifically, the system is entrained by a control u=v(Ωt) when the phase drift equation (10x) satisfies {dot over (φ)}=0, which will occur as t.fwdarw.∞ if there exists a phase φ.sup.∞ that satisfies Δω=Λ.sub.v(φ.sup.∞)=0. Conversely, when both the control waveform v and PRC Z are non-zero, the function Λ.sub.v(φ) is not identically zero, so when the system is entrained there exists at least one φ.sup.∞ε[0,2π) that is an attractive fixed point of equation (10x). Such stable fixed points {φ.sub.i.sup.∞} of equation (10x) are the roots of the equation Δω+Λ.sub.v(φ)=0 that also satisfy Λ.sub.v′(φ)<0, where

[00006]${\Lambda }_v^{\prime }\left(\varphi \right)=\frac{d}{d.Math..Math.\varphi }.Math.{\Lambda }_v.$

The values {φ.sub.i.sup.∞} determine the average phase shift, relative to the forcing phase θ=Ωt, at which the oscillation stabilizes from an initial phase difference φ(0). It is denoted by (φ.sub.i.sup.∞)⊂[0,2π) that the set of initial phases φ(0) that result in convergence of the oscillator to φ.sub.i.sup.∞. In addition, the phases φ.sup.+=arg max.sub.φΛ.sub.v(φ) and φ.sup.−=arg min.sub.φΛ.sub.v(φ) at which the interaction function achieves its maximum and minimum values, respectively, are defined. In order for entrainment to occur, −Λ.sub.v(φ.sup.+)≦Δω≦−Λ.sub.v(φ.sup.−) must hold, so that at least one stable fixed point of Λ.sub.v exists. Thus the range of the interaction function determines which values of the frequency detuning Δω yield phase locking, and the shape of Λ.sub.v determines the possible asymptotic phase shifts φ.sub.i.sup.∞. These properties are illustrated in FIG. 7.

iv) Ensemble Interaction Function and Control Input Construction

[0102] The interaction function of a structurally similar ensemble is examined, and methods of constructing and processing the interaction function to create a control input with a desired phase selection effect is disclosed herein. The interaction function formalism is invoked to derive the long-run behavior of a collection of N nonlinear oscillators with phase-reduced dynamics of the form in equation (2x) given by

={{dot over (ψ)}.sub.j=ω.sub.j+Z(ψ.sub.j)u,j=1, . . . ,N}  (11x)

which all share an identical PRC. It is supposed, without loss of generality, that the natural frequencies of the ensemble elements are ordered according to ω.sub.1<ω.sub.2< . . . ω.sub.N. If an input of the form u=v(Ωt), where v is 2π-periodic, is applied to all the oscillators in , the collection of average phase drift dynamics of the form in equation (10x) is obtained as given by

={{dot over (φ)}.sub.1=Δω.sub.j+Λ.sub.v(φ.sub.j), j=1, . . . ,N}.  (12x)

[0103] If each average drift equation in has a fixed point, then the input waveform v(Ωt) entrains all of the oscillators in to the frequency Ω, and the pattern that emerges in the phase offsets φ.sub.j.sup.∞ of the oscillators relative to the forcing phase θ=Ωt can be inferred from the interaction function, the frequencies ω.sub.j, and the initial conditions φ.sub.j(0). Specifically, the average phase drift {dot over (φ)}.sub.j of each synchronized oscillator will be zero, and the average phase offset from the forcing phase θ=Ωt will be φ.sub.ij.sup.∞, depending on which set (φ.sub.ij.sup.∞) the initial phase φ.sub.j(0) belongs to. Here .sub.j(φ.sub.ij.sup.∞) denotes the set of initial phases for which the j.sup.th oscillator of is attracted to the asymptotic phase offset φ.sub.ij.sup.∞. When synchronization of the collection to the fixed points φ.sub.ij.sup.∞ occurs, then {dot over (φ)}.sub.j=Δω.sub.i+Λ.sub.v(φ.sub.j)=0 holds, with

Λ.sub.v(φ.sub.ij.sup.∞)=Δω.sub.j j=1, . . . ,N  (13x)

Λ.sub.v′(φ.sub.ij.sup.∞)<0 j=1, . . . ,N,  (14x)

where φ.sub.ij.sup.∞ is one of the possible asymptotic phases for the jth oscillator.

[0104] The condition in equation (14x) guarantees that the dynamical configuration is locally attractive and stable, as illustrated for an individual oscillator in FIG. 7. Conversely, if Λ.sub.v satisfies equations (13x)-(14x), and the initial phases of the oscillators in satisfy ψ.sub.j(0)=φ.sub.ij.sup.∞ for j=1, . . . , N, then the synchronization pattern will be maintained. The pattern is also established and maintained when the initial conditions are relaxed to φ.sub.j(0)ε(φ.sub.ij.sup.∞) for j=1, . . . , N. These regions are illustrated for an example interaction function in FIG. 8. Note that for a given collection of entrained oscillators in equation (12x) the stable synchronization pattern may not be unique, but depends on the initial conditions ψ.sub.j(0) of the collection.

[0105] The above analysis describes how the interaction function between a forcing waveform v and the PRC Z common to the collection characterizes the asymptotic phase structure of an oscillator ensemble with heterogeneous natural frequencies. The construction of ideal interaction functions for desired phase patterns is then examined, and the properties of the ensemble that determine their realizability. Because the Z(θ), v(θ), and θ.sub.v(φ) are all 2π-periodic, they may be expressed using Fourier series, and the Fourier series for Λ.sub.v(φ) is readily computed using equation (9x), trigonometric identities, and the orthogonality of the Fourier basis. The functions Z and v may be represented using the truncated series expansions

[00007]$\begin{array}{cc}Z\left(\theta \right)\approx Z^r\left(\theta \right)=\frac{a_0}{2}+{.Math.}_{n=1}^r.Math.\left[a_n.Math..Math.\cos .Math..Math.\left(n.Math..Math.\theta \right)+b_n.Math..Math.\sin \left(n.Math..Math.\theta \right)\right],& \left(15.Math.x\right)\\ v\left(\theta \right)\approx v^r\left(\theta \right)=\frac{c_0}{2}+{.Math.}_{n=1}^r.Math.\left[c_n.Math..Math.\cos .Math..Math.\left(n.Math..Math.\theta \right)+d_n.Math..Math.\sin \left(n.Math..Math.\theta \right)\right],& \left(16.Math.x\right)\end{array}$

where the appropriate order r will be discussed below. Applying trigonometric angle sum identities to equation (9) and the orthogonality of the Fourier basis to eliminate terms yields

[00008]$\begin{array}{cc}{\Lambda }_v^r\left(\varphi \right)=\frac{f_0}{2}+\frac{1}{2}.Math.{.Math.}_{n=1}^r.Math.\left[f_n.Math..Math.\cos .Math..Math.\left(n.Math..Math.\varphi \right)+g_n.Math..Math.\sin \left(n.Math..Math.\varphi \right)\right],& \left(17.Math.x\right)\\ f_0=\frac{a_0.Math.c_0}{2},.Math.f_0=a_n.Math.c_n+b_n.Math.d_n,.Math.g_n=b_n.Math.c_n-a_n.Math.d_n.& \left(18.Math.x\right)\end{array}$

[0106] Therefore given truncated Fourier series expansions Λ.sub.v.sup.r(θ) and Z.sup.r(θ), the Fourier coefficients of the corresponding truncated series for the control waveform v.sup.r(θ) are given by

[00009]$\begin{array}{cc}{c}_0=2.Math.\frac{f_0}{a_0}.Math.{\chi }_{\left[a_0\ne 0\right]},.Math.c_n=2.Math.\frac{f_n.Math.a_n+b_n.Math.g_n}{a_n^{2+}.Math.b_n^2}.Math.{\chi }_{\left[a_n^2+b_n^2\ne 0\right]},.Math.d_n=2.Math.\frac{f_n.Math.b_n-a_n.Math.g_n}{a_n^2+b_n^2}.Math.{\chi }_{\left[a_n^2+b_n^2\ne 0\right]},& \left(19.Math.x\right)\end{array}$

where χ.sub.A=1 if A is true, and χ.sub.A=0 otherwise.

[0107] An ideal interaction function Λ.sub.v(φ) that corresponds to the desired phase assignment task in equations (13x)-(14x) can be designed using a sum of scaled and shifted sigmoid functions. The error function is used, which is given by

[00010]$\begin{array}{cc}\mathrm{erf}\left(x\right)=\frac{2}{\sqrt{\pi }}.Math.{\int }_0^x.Math.e^{-t^2}.Math.\mathrm{dt},.Math.& \left(20.Math.x\right)\end{array}$

as the basic element of the construction of the designed ideal interaction function that passes precisely through the coordinates (φ.sub.j*, Δω.sub.j). A smooth approximation to the unit step function is first defined by

σ(x)=½(erf(2x)+1),  (21x)

which satisfies σ(−1)≈0, σ(1)≈1, and σ(0)=½. The range points are then created,

r.sub.1=−Δω.sub.1+½(−Δω.sub.1+−Δω.sub.2)  (22x)

r.sub.j=½(−Δω.sub.j−Δω.sub.j−1), j=2, . . . ,N  (23x)

r.sub.N+1=−Δω.sub.N+½(−Δω.sub.N−1+Δω.sub.N),  (24x)

as well as the domain points,

s.sub.j=½(φ.sub.1*+φ.sub.N*−2π),  (25x)

s.sub.j=½(φ.sub.j*+φ.sub.j−1*), j=2, . . . ,N  (26x)

s.sub.N+1=½(φ.sub.1*+φ.sub.N*+2π),  (27x)

and the domain scaling factors

h.sub.1=min{s.sub.2−φ.sub.1*,(φ.sub.1*−s.sub.1)/2},  (28x)

h.sub.j=min{s.sub.j+1−φ.sub.1*,(φ.sub.1*−s.sub.1)/2}, j=2, . . . N  (29x)

h.sub.N+1=min{½(s.sub.N+1−φ.sub.N*,}.  (30x)

A suitable interaction function Λ.sub.v* is then constructed by

[00011]$\begin{array}{cc}{\Lambda }_v^{*}\left(\varphi \right)=r_1+\left(r_1-r_{N+1}\right).Math.\sigma \left(\left(\varphi -s_1\right)/h_{N+1}\right)+\underset{j=2}{\overset{N}{.Math.}}.Math.\left(r_{j+1}-r_j\right).Math.\sigma \left(\left(\varphi -{\varphi }_j^{*}\right)/h_j\right)+\underset{j=2}{\overset{N}{.Math.}}.Math.\left(r_{j+1}+r_j\right).Math.\sigma \left(\left(\varphi +2.Math.\pi -{\varphi }_j^{*}\right)/h_j\right)+\left(r_1-r_{N+1}\right).Math.\sigma \left(\left(\varphi -s_{N+1}\right)/h_{N+1}\right).Math.& \left(31.Math.x\right)\end{array}$

[0108] The function in equation (31x) will thus pass close to the phase and detuning pairs (φ.sub.j,−Δω.sub.j). The sums in (31x) are repeated with the added offset factors 2π and −2π because the sigmoid function (21x) is not a perfect Heaviside unit step, but a periodic interaction function Λ.sub.v* must nevertheless be created. The particular methodology used to automatically create a suitable interaction function Λ.sub.v* that approximately satisfies (13x)-(14x) is not crucial in itself.

[0109] As long as the designed interaction function satisfies these properties, the Fourier coefficients (19x) obtained by circular de-convolution of this interaction function Λ.sub.v* and the nominal ensemble PRC Z yields the control that produces the desired phase pattern. Equations (20x)-(31x) are provided as a suggested example algorithm, which has been used for automatic control design for the experiments disclosed herein. A satisfactory interaction function curve Λ.sub.v* may not be unique, and could be constructed using many different algorithms, and the need to satisfy (13x)-(14x) is the only common requirement. Control inputs may be designed indirectly by constructing an interaction function.

[0110] The resulting function Λ.sub.v* is approximated by a truncated r-order Fourier series Λ.sub.v.sup.r as in equation (17x), which is then used to compute the input waveform in equation (19x). The maximal order r is determined by the number of square amplitudes .sub.n.sup.2 of the Fourier series terms in the PRC that exceed a tolerance δ.sub.Z. In order for Λ.sub.v.sup.r(φ) to satisfy the design conditions in equations (13x)-(14x), it must sufficiently well approximate Λ.sub.v(φ), and therefore the realizability of the design in equations (13x)-(14x) ultimately depends on the complexity of the PRC, in addition to the frequency distribution of the ensemble and the desired phase offsets. When these conditions are satisfied, the control waveform v.sup.r(θ) that results in the desired phase pattern can be easily synthesized using the Fourier coefficient formula in equation (19x).

iv) Interaction Function Realizability

[0111] The nonlinear complexity of the common PRC of oscillating subsystems in an ensemble determines the realizability of a desired interaction function for the ensemble. Specifically, a greater number of significant terms in the Fourier series approximation of the PRC are obtained as the oscillation exhibits higher relaxational behavior. Typically, this occurs as the oscillation moves farther past a Hopf bifurcation.

[0112] By inspecting equation (19x), it is evident that v.sup.r(θ) exists for a pair Λ.sub.v.sup.r(θ) and Z.sup.r(θ) only if a.sub.n.sup.2+b.sub.n.sup.2>0 for all n for which f.sub.n>0 or g.sub.n>0, and when a.sub.0≠0 if f.sub.0≠0. Therefore, any terms of the Fourier series for Λ.sub.v for which these conditions are not satisfied must be removed from the construction. For actual biological and chemical oscillating systems, the magnitudes .sub.n.sup.2=a.sub.n.sup.b+b.sub.n.sup.2 decay approximately exponentially. It is thus crucial for the expansion order r to be appropriately chosen such that .sub.n.sup.2>δ.sub.Z for n≦r, where δ.sub.Z is an appropriate tolerance, and any order n for which .sub.n=0 occurs must be omitted from the truncated series. This guarantees that numerical conditioning errors do not arise, and also that the designed v.sup.r(θ) satisfies the weak forcing assumption. An appropriate empirical value for the tolerance is δ.sub.Z=Z.sup.2.Math.10.sup.−4. Hence given the input u(t)=v.sup.r(Ωt), the asymptotic configuration of the entrained oscillators of equation (12x) will approximately satisfy the conditions in equations (13x)-(14x).

[0113] For the experimental electrochemical oscillators described in the main text, the appropriate series truncation is r=5. This limits the types of phase patterns that can be achieved for ensembles of such oscillators, and in particular, explains why the observed and desired values of Δφ in the parity diagram in FIG. 3C of the main text diverge near Δφ=⅘.Math.2π. The interaction function must be decreasing at Λ.sub.v(φ.sub.2*)=Δω.sub.2 then return and again decrease at Λ.sub.v(φ.sub.1*)=−Δω.sub.1. Because only 5 Fourier series terms can be used to synthesize Λ.sub.v, the derivative

[00012]$\frac{d}{d.Math..Math.\varphi }.Math.{\Lambda }_v\left(\varphi \right)$

is limited when v has limited amplitude, so there is a minimum phase gap required for this return. This capacity for return of Λ.sub.v disappears as Δφ approaches 2π, as seen in the left panel of FIG. 3E of the main text.

[0114] The realizability of the interaction function in the above sense corresponds directly to realizability of a particular phase pattern in the ensemble. Thus, the existence of a periodic waveform that entrains an ensemble to such a phase pattern depends on the PRC of the ensemble sub-elements. In particular, a greater number of significant terms in the Fourier series expansion of Z corresponds to greater flexibility in phase pattern construction. Thus, more highly nonlinear and relaxational oscillation dynamics reflect greater controllability of the ensemble. Further work to more rigorously characterize this property of oscillatory ensembles is compelling.

v) Interaction Functions for Globally Attractive Patterns

[0115] A control input waveform is globally attractive for a phase pattern if frequency detuning values for the oscillators are monotone decreasing as desired phase offsets increase, and the corresponding interaction function can be realized. The natural frequencies of the ensemble may be ordered according to ω.sub.1<ω.sub.2< . . . ω.sub.N. Then the frequency detuning values in satisfy −Δω.sub.1>−Δω.sub.2> . . . Δω.sub.N. Suppose also that the desired phase assignment satisfies φ.sub.1*<φ.sub.2*< . . . φ.sub.N*. Then an ideal interaction function that satisfies the appropriate conditions in equations (13x)-(14x) can be designed to have a unique set of asymptotic average phase offsets φ.sub.j.sup.∞=φ.sub.j*; that are globally attractive, with .sub.j(φ.sub.j*)=[0,2π). The desired interaction function is monotone decreasing for phases on [δ.sub.1,φ.sub.N*+δ.sub.2)(mod2π) and monotone increasing on (φ.sub.N*+δ.sub.2,2π) to satisfy periodicity, where δ.sub.1 and δ.sub.2 are sufficient gaps for

[00013]$\frac{d}{d.Math..Math.\varphi }.Math.{\Lambda }_v\left(\varphi \right)$

to change sign. An example of such a Λ.sub.v, which is of the same form as those used in experiments described herein previously is illustrated in FIG. 9.

[0116] Application of the waveform v that produces the interaction function Λ.sub.v in FIG. 9 through the circular convolution in equation (9x) will asymptotically guide a rhythmic ensemble of four oscillators into a coherent phase configuration with phase offsets φ.sub.11.sup.∞, φ.sub.12.sup.∞, φ.sub.13.sup.∞, and φ.sub.14.sup.∞ for oscillators j=1, 2, 3, 4. In essence, the heterogeneity in natural frequencies ω.sub.j and the nonlinearity in Z are leveraged to design a weak, low-amplitude waveform that achieves and maintains such a synchronization structure without altering the fundamental dynamics of individual rhythmic units. Monotone ordering of the assigned phases in the same order as natural frequencies allows the synthesis of an ideal interaction function that achieves the desired stable fixed-point average phases φ.sub.1j.sup.∞=lim.sub.t.fwdarw.∞φ.sub.j(t)=φ.sub.j*. Moreover, such a Λ.sub.v satisfies .sub.j(φ.sub.j*)=[0,2π) for all j=1, . . . , N so that the pattern will be established for any set of initial conditions Imo), so that the phase structure is globally attractive. Finally, suppose that the coefficients of the Fourier series of the PRC Z of the oscillators in satisfy .sub.n.sup.2>δ.sub.Z=Z.sup.2.Math.10.sup.−4 for all n<N, so that the deconvolution of Λ.sub.v as given in equations (15x)-(19x) may be done using truncated series with r=N terms.

Precursor Waveforms for Non-Globally Attractive Patterns

[0117] If the desired phase pattern is not monotonically arranged in the sense described herein below, then a control waveform does not exist such that the pattern is globally attractive. In that case, a precursor waveform is required to bring the phases of ensemble elements into an arrangement after which an ultimate waveform is applied to finalize and hold the phase pattern. When the desired phase offsets φ.sub.j* corresponding to natural frequencies ω.sub.j are not uniformly increasing on [0,2π), a globally attractive phase structure cannot be achieved. In such cases, as the one illustrated in FIG. 8, one or more precursor inputs must be applied to steer subsets of the ensemble elements into the desired attractive regions .sub.j(φ.sub.j*), after which a final input is applied to finalize and subsequently hold the pattern. For the example in FIG. 8, a precursor waveform can be applied to entrain oscillators j=2 and 3 to designed phase offsets φ.sub.2*=φ.sub.22.sup.∞ and φ.sub.3*=φ.sub.13.sup.∞ as shown in FIG. 10. When this preliminary entrainment step is accomplished with respect to the chosen forcing phase, the control input is changed to the one corresponding to the Λ.sub.v for the full pattern in FIG. 8, which is outlined using a dashed line in FIG. 10 as well. This two-step procedure guarantees that the final outcome is the desired phase assignment, rather than one of several possible patterns. This is the situation that characterizes the experiment “O.fwdarw.K.fwdarw.O” as illustrated in FIG. 5A, FIG. 5B, FIG. 5C, FIG. 5D, FIG. 5E, FIG. 5F, FIG. 5G, FIG. 5H, FIG. 5I, and FIG. 5J and FIG. 6A, FIG. 6B, FIG. 6C, FIG. 6D, FIG. 6E, FIG. 6F, FIG. 6G, FIG. 6H, FIG. 6I, FIG. 6J, FIG. 6K, FIG. 6L, FIG. 6M, FIG. 6N, and FIG. 6O.

vii) Interaction Function for Desynchronization

[0118] Suppose it is wished to accomplish phase desynchronization of a very large ensemble of oscillators , where N is large but uncertain with the frequencies ω.sub.j distributed uniformly on an interval [φ.sub.min,φ.sub.max] The ideal interaction function Λ.sub.v corresponding to a waveform that will spread the phase offsets of the ensemble uniformly across a large domain is shown in FIG. 11. Note that because the interaction function satisfies the monotonicity property, the desynchronizing effect is global for any initial configuration of the ensemble. The Fourier series coefficients of the desynchronizing control waveform v are then obtained from the formula in equation (19x) for circular deconvolution from Λ.sub.v. A natural continuation of the work on generation of coherent phase structures using global controls without feedback disclosed herein pertains to the effect of such controls on networks of weakly coupled oscillators. Specifically, in the presence of weak coupling, the averaged dynamic equation in (12x) for the j.sup.th oscillator becomes {dot over (φ)}.sub.j=Δω.sub.j+Λ.sub.v(φ.sub.j)+εΣ.sub.iH.sub.ij(φ.sub.j−φ.sub.i), where H.sub.ij characterizes the interaction between the i.sup.th and j.sup.th oscillators. If the coupling strength is weak relative to the interaction function Λ.sub.v, i.e., 0<ε<<1, then the coupling terms will have little influence on the ability to construct phase patterns in the ensemble by external inputs. A mathematical problem is posed by characterizing the desynchronization action of this type of control and how it is affected as coupling strength c increases. A phase transition may occur after a critical coupling strength is reached after which the desynchronizing action is overpowered by synchronization due to coupling. As a result, a continued investigation in the direction of effectively desynchronizing entrainment controls that account for couplings and topology is warranted. In the context of affecting pathologically synchronized neuronal ensembles, progress in this direction will have particularly important implications for research and clinical neuroscience applications.

[0119] Consider that a simple sinusoidal forcing of the form v(Ωt)=sin(Ωt) can be used to desynchronize the ensemble with phase offsets on [0,π). This is because it results in an interaction function that decreasing on an interval of length π so the maximum phase difference achievable is |φ.sub.max.sup.∞−φ.sub.min.sup.∞|<π. Hence, the disclosed approach enables more versatile manipulation of phase relationships beyond this limitation. In fact, it is possible to quantify how the disclosed approach increases the achievable relative phase desynchronization difference relative to sinusoidal forcing. If the PRC has non-trivial terms in its Fourier series up to order r, then it is possible to create an input waveform that results in an interaction function that is monotone on the interval

[00014]$\left[0,\left(2-\frac{1}{r}.Math.\pi \right)\right].$

Thus, an approach can be used to achieve maximal desynchronization of an oscillator ensemble by increasing the phase offset distribution range by a factor of

[00015]$2-\frac{1}{r}$

period.

NEW1. Extension of the Phase-Selective Entrainment Technique

[0120] The disclosed method may be used to selectively entrain an oscillator network in which the oscillators are weakly coupled. Specifically, consider the network

{dot over (φ)}.sub.j=Δω.sub.j+Λ.sub.v(φ.sub.j)+Σ.sub.i=1.sup.Nα.sub.ijK.sub.ij(φ.sub.j−φ.sub.i), j=1, . . . ,N,  (1z)

where α.sub.ij and K.sub.ij are the coupling influence of oscillator i on oscillator j and the coupling function between them, respectively, and

[00016]${\Lambda }_v\left(\varphi \right)={\lim }_{T.fwdarw.\infty }.Math.\frac{1}{T}.Math.{\int }_0^T.Math.Z\left(\Omega .Math..Math.t+\varphi \right).Math.v\left(\Omega .Math..Math.t\right).Math.\mathrm{dt}$

is the interaction function in which Z is the phase response curve (PRC). The network of oscillators is defined as synchronized if the conditions

(i) Λ.sub.v(φ.sub.j.sup.∞)=Δω.sub.j−Σ.sub.i=1.sup.Nα.sub.ijK.sub.ij(φ.sub.j−φ.sub.i),  (2z)

(ii) Λ.sub.v′(φ.sub.j.sup.∞)<−α.sub.ijK.sub.ij′(φ.sub.j−φ.sub.i), j=1, . . . ,N,  (3z)

hold, where condition (ii) guarantees that the dynamical configuration is locally attractive and stable. Note that if the oscillators are not connected, then the coupling functions K.sub.ij do not exist so that the equation (1z) and conditions in (2z) and (3z) are identical to the uncoupled case, namely, Λ.sub.v(φ.sub.j.sup.∞)=−Δω.sub.j and Λ.sub.v′(φ.sub.j.sup.∞)<0.

[0121] In various aspects, the methods described herein may be implemented using a remote and/or local computing device as described herein below. FIG. 16 illustrates an example configuration of a remote device system 300 and depicts an exemplary configuration of a remote or user computing device 302, such as requestor computing device. Computing device 302 includes a processor 304 for executing instructions. In some embodiments, executable instructions are stored in a memory area 306. Processor 304 may include one or more processing units (e.g., in a multi-core configuration). Memory area 306 is any device allowing information such as executable instructions and/or other data to be stored and retrieved. Memory area 306 may include one or more computer-readable media.

[0122] Client computing device 302 also includes at least one media output component 308 for presenting information to a user 310. Media output component 308 is any component capable of conveying information to user 310. In some embodiments, media output component 308 includes an output adapter such as a video adapter and/or an audio adapter. An output adapter is operatively coupled to processor 304 and operatively coupleable to an output device such as a display device (e.g., a liquid crystal display (LCD), organic light emitting diode (OLED) display, cathode ray tube (CRT), or “electronic ink” display) or an audio output device (e.g., a speaker or headphones). In some embodiments, media output component 308 is configured to present an interactive user interface (e.g., a web browser or client application) to user 310.

[0123] In some embodiments, client computing device 302 includes an input device 312 for receiving input from user 310. Input device 312 may include, for example, a keyboard, a pointing device, a mouse, a stylus, a touch sensitive panel (e.g., a touch pad or a touch screen), a camera, a gyroscope, an accelerometer, a position detector, and/or an audio input device. A single component such as a touch screen may function as both an output device of media output component 308 and input device 312.

[0124] Computing device 302 may also include a communication interface 314, which is communicatively coupleable to a remote device such as SE computing device. Communication interface 314 may include, for example, a wired or wireless network adapter or a wireless data transceiver for use with a mobile phone network (e.g., Global System for Mobile communications (GSM), 3G, 4G, or Bluetooth) or other mobile data network (e.g., Worldwide Interoperability for Microwave Access (WIMAX)).

[0125] Stored in memory area 306 are, for example, computer-readable instructions for providing a user interface to user 310 via media output component 308 and, optionally, receiving and processing input from input device 312. A user interface may include, among other possibilities, a web browser and client application. Web browsers enable users 310 to display and interact with media and other information typically embedded on a web page or a website from a web server associated with a merchant. A client application allows users 310 to interact with a server application.

[0126] FIG. 17 illustrates an example configuration of a server system 400, and depicts an example configuration of a server computing device 402. Server computing device 402 includes a processor 404 for executing instructions. Instructions may be stored in a memory area 406, for example. Processor 404 may include one or more processing units (e.g., in a multi-core configuration).

[0127] Processor 404 is operatively coupled to a communication interface 408 such that server computing device 402 is capable of communicating with a remote device such as computing device 302 shown in FIG. 16 or another server computing device 402. For example, communication interface 408 may receive requests from requestor computing device 106 via the Internet.

[0128] Processor 404 may also be operatively coupled to a storage device 410. Storage device 410 is any computer-operated hardware suitable for storing and/or retrieving data. In some embodiments, storage device 410 is integrated in server computing device 402. For example, server computing device 402 may include one or more hard disk drives as storage device 410. In other embodiments, storage device 410 is external to server computing device 402 and may be accessed by a plurality of server computing devices 402. For example, storage device 410 may include multiple storage units such as hard disks or solid state disks in a redundant array of inexpensive disks (RAID) configuration. Storage device 410 may include a storage area network (SAN) and/or a network attached storage (NAS) system.

[0129] In some embodiments, processor 404 is operatively coupled to storage device 410 via a storage interface 412. Storage interface 412 is any component capable of providing processor 404 with access to storage device 410. Storage interface 412 may include, for example, an Advanced Technology Attachment (ATA) adapter, a Serial ATA (SATA) adapter, a Small Computer System Interface (SCSI) adapter, a RAID controller, a SAN adapter, a network adapter, and/or any component providing processor 404 with access to storage device 410.

[0130] Memory areas 306 (shown in FIG. 3) and 406 may include, but are not limited to, random access memory (RAM) such as dynamic RAM (DRAM) or static RAM (SRAM), read-only memory (ROM), erasable programmable read-only memory (EPROM), electrically erasable programmable read-only memory (EEPROM), and non-volatile RAM (NVRAM). The above memory types are for example only, and are thus not limiting as to the types of memory usable for storage of a computer program.

EXAMPLES

[0131] The following examples demonstrate various aspects of the disclosure.

Example 1: Electrochemical Dissolution of Nickel Using Potential Actuation

[0132] To demonstrate the method for controlling the phase structure of an ensemble of heterogeneous oscillators, the following experiments were conducted. Electrochemical oscillations associated with the electrochemical dissolution of nickel in concentrated sulfuric acid solution using potential actuation was controlled using the method described above.

[0133] A schematic description of the experimental setup is depicted in FIG. 1. As illustrated in FIG. 1, the working electrode (WE) was composed of an array of nickel wires, the counter electrode (CE) was a Pt coated Ti wire, and the reference electrode (RE) was Hg/Hg2SO4/sat. K2SO4. The control waveform was a potential signal generated by a high-speed multifunction data acquisition (DAQ) device (M Series DAQ PXI-6255, National Instruments). The control waveform was implemented in the WE by a potentiostat (Bank Instruments) as a superimposed signal on the applied baseline potential V0. Each nickel wire was connected in series to an individual resistance Rind.

[0134] The apparatus consisted of twenty nickel wires that functioned 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, each nickel wire of the WE was polished with sandpaper with one of six levels of roughness, ranging from 180 to 4000 grit. The current of each working electrode was monitored independently. During the anodic dissolution of nickel within the WE, the spontaneous formation of self-sustained electrochemical oscillations driven by the presence of a negative differential resistance was observed.

[0135] Once the WE were placed in the electrochemical cell, a slow positive sweep of 0.01 V s−1 from 0 V to a constantly applied potential V0 was conducted in order to form a thin passive layer on the electrode surface. This baseline was set in reference to an Hg/Hg2SO4/sat. K2SO4 reference electrode (RE) in the electrochemical cell, containing a 1.6 mm diameter Pt-coated Ti wire counter electrode (CE), and 3.0 mol L−1 sulfuric acid solution at constant temperature of 10 deg C. The potential V0 was initially set using the potentiostat (Bank Instruments) at a value for which the oscillation was close to the Hopf bifurcation, which was approximately 1.15 V. Inputs to the oscillating system consisted of an additional potential V superimposed onto the baseline potential V0 using the potentiostat.

[0136] Soon after, the PRCs were measured simultaneously for the WE by the automatic procedure of applying a pseudo-randomly timed potential pulse sequence of −0.20 V with pulse-width of 0.05 seconds and post-processing the observed current using the pulse perturbation procedure disclosed herein above. Measurements of the current oscillations in the electrochemical reactions were obtained by an interface to the high-speed multifunction data acquisition (DAQ) device with a sample rate of 200 Hz. Simultaneously, the periodic potential perturbations v with determined amplitude scaling A and frequency Ω were applied in the electrochemical oscillator ensemble using the controlled potentiostat to superimpose the waveform onto the applied baseline voltage V0. Both interfaces were written in LabView. Each control waveform was generated based on the PRC obtained prior to conducting the experiment and the targeted interaction function was generated from the desired set of phase offsets for the oscillators.

[0137] FIG. 2A, FIG. 2B, FIG. 2C, FIG. 2D, FIG. 2E, and FIG. 2F are graphs showing a series of interaction functions and control waveforms used to enable the switching pattern “O.fwdarw.K.fwdarw.O” as described in additional detail below. FIG. 2A, FIG. 2C, and FIG. 2E illustrate the best approximations to the desired interaction functions, and FIG. 2B, FIG. 2D, and FIG. 2F; Right: one cycle of the resulting control waveform, obtained using the control design procedure described herein, which is applied to entrain the electrochemical oscillators. The ensemble is grouped into four clusters of quantities (N1;N2;N3;N4)=(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 pattern “O”: (φ.sub.1*,φ.sub.2*,φ.sub.3*,φ.sub.4*)=(0,0,π,π); (b): Precursor of pattern “K”: (φ.sub.2*,φ.sub.3*)=(π, 0); (c): Pattern “K”: (φ.sub.1*,φ.sub.2*,φ.sub.3*,φ.sub.4*)=(0, π, 0, π). FIG. 2A, FIG. 2C, and FIG. 2E display the interaction functions and FIG. 2B, FIG. 2D, and FIG. 2F display the respective control waveforms (used to entrain the electrochemical oscillators for the pattern switching procedure OKO described in FIG. 5A, FIG. 5B, FIG. 5C, FIG. 5D, FIG. 5E, FIG. 5F, FIG. 5G, FIG. 5H, FIG. 5I, and FIG. 5J. An electrochemical oscillator ensembles with clusters of size (N1;N2;N3;N4)=(7;7;3;3) were selected by tuning the mean natural frequency of the oscillations in the following order: cluster 1 with ω.sub.1=0:390 Hz (electrodes j=1 to 7 in blue), cluster 2 with φ.sub.2=0:406 Hz (electrodes j=8 to 14 in cyan), cluster 3 with ω.sub.3=0:427 Hz (electrodes j=15 to 17 in yellow) and cluster 4 with ω.sub.4=0:442 Hz (electrodes j=18 to 20 in red). The phase assignments for pattern switching OKO are listed in FIG. 5A-5J.

Example 2: Phase Assignment Applications to Different Systems

[0138] The phase assignment method disclosed herein was applied to entrain oscillatory systems arising from different domains. In particular, robust global entrainment waveforms were designed to form phase patterns (synchronization, desynchronization, and clustering) in a collection of oscillatory systems, including Hodgkin-Huxley neurons, Morris-Lecar neurons, and electro-chemical oscillators.

Entrainment of Hodgkin-Huxley Neurons

[0139] The Hodgkin-Huxley (HH) model (see Appendix A) describes the propagation of action potentials in neurons, specifically the squid giant axon, and is used as a canonical example of neural oscillator dynamics. The phase-selective entrainment technique was applied to entrain HH neurons.

Hodgkin-Huxley Model

[0140] The Hodgkin-Huxley model describes the propagation of action potentials in neurons, specifically the squid giant axon, and is used as a canonical example of neural oscillator dynamics. The equations are

{dot over (V)}=I.sub.b+I(t)−g.sub.NaH(V−V.sub.Na)m.sup.3−g.sub.K(V−V.sub.k)n.sup.4−g.sub.L(V−V.sub.L),

{dot over (m)}=a.sub.m(V)(1−m)−b.sub.m(V)(m),

{dot over (h)}=a.sub.h(V)(1−h)−b.sub.h(V)(h),

{dot over (n)}=a.sub.n(V)(1−n)−b.sub.n(V)(n),

a.sub.m(V)=0.1(V+40)/(1−exp(−(V−40)/10)),

b.sub.m(V)=4exp(−(V+65)/18),

a.sub.h(V)=0.07exp(−(V+65)/20),

b.sub.h(V)=1/(1+exp(−(V+35)/10)),

a.sub.n(V)=0.01(V+55)/(1−exp(−(V+55)/10)),

b.sub.n(V)=0.125exp(−(V+65)/80).

[0141] The variable V is the voltage across the axon membrane, m, h, and n are the ion gating variables, I.sub.b is a baseline current that induces the oscillation, and I(t) is the control input. The units of V are millivolts and the units of time are milliseconds. Nominal parameters are V.sub.Na=50 mV, V.sub.K=−77 mV, V.sub.L=−54.4 mV, g.sub.Na=120 mS/cm.sup.2, g.sub.K=36 mS/cm.sup.2, g.sub.L=0.3 mS/cm.sup.2, I.sub.b=10 μA/cm.sup.2, and c=1 μF/cm.sup.2, for which the period of oscillation is T=14.63842±10.sup.−5 ms.

Uncoupled HH Neurons

[0142] An entrainment phase pattern {φ.sub.1.sup.∞,φ.sub.2.sup.∞,φ.sub.3.sup.∞}={π/3, π, 5π/3} for three oscillators (N=3) with the frequencies {ω.sub.1, ω.sub.2, ω.sub.3}={0.99ω, ω, 1.01ω} was created, where ω is the natural frequency. The designed interaction function N.sub.v.sup.r(φ) and the optimal waveform v.sup.r(θ) are shown in FIG. 12A and FIG. 12B, respectively. The desired phase structure stabilizes after about 50 input cycles, as shown in FIG. 12C.

HH Neuronal Networks

[0143] The same entrainment phase pattern {φ.sub.1.sup.∞,φ.sub.2.sup.∞,φ.sub.3.sup.∞}={π/3, π, 5π/3} was created for the same oscillators as described above with their frequencies {ω.sub.1, ω.sub.2, ω.sub.3}={0.99ω, ω, 1.01ω} and consider the case when they are weakly coupled through the coupling functions K.sub.ij(φ)=cos(φ) with the coupling strength α.sub.ij=0.05 for i,j=1, 2, 3. The designed interaction function and the optimal waveform Λ.sub.v.sup.r(v) are shown in FIG. 13A and FIG. 13B, respectively. The desired phase structure stabilizes after about 50 input cycles, as shown in FIG. 13C.

Entrainment of Morris-Lecar Neurons

[0144] The Morris-Lecar model was originally proposed to capture the oscillating voltage behavior of giant barnacle muscle fibers. Over the past years this model has been extensively studied and used as a standard model for representing many different real neurons that are experimentally observable. For example, it has been found that Morris-Lecar PRC is extremely similar to the experimentally observed PRCs of Aplysia motoneuron.

Morris-Lecar Model

[0145] The dynamics of the Morris-Lecar neuron is described by

C{dot over (V)}+(I.sup.b+I)=g.sub.Cam.sub.∞(V.sub.Ca−V)+g.sub.Kw(V.sub.K−V)+g.sub.Lw(V.sub.L−V),

{dot over (w)}=φ(ω.sub.∞−ω)/τ.sub.ω(V),

where m.sub.∞=0.5[1+tan h((V−V.sub.1)/V.sub.2)], ω.sub.∞=0.5[1+tan h((V−V.sub.3)/V.sub.4)], τ.sub.ω=1/cos h[(V−V.sub.3)/(2V.sub.4)], and parameter values for the Morris-Lecar neuron are given as:

[0146] φ=0.5, I.sup.b=0.09 ρA/cm.sup.2, V.sub.L=−0.01 mV, v.sub.2=0.15 mV, V.sub.3=0.1 mV, v.sub.4=0.145 mV, g.sub.Ca=1 mS/cm.sup.2, V.sub.K=−0.7 mV, V.sub.L=−0.5 mV, g.sub.K=2 mS/cm.sup.2, g.sub.L=0.5 mS/cm.sup.2, C=1 μF/cm.sup.2.

[0147] A uniform entrainment pattern was created {φ.sub.1.sup.∞,φ.sub.2.sup.∞,φ.sub.3.sup.∞,φ.sub.4.sup.∞,φ.sub.5.sup.∞}={π/5, 3π/5, π, 3π/5, 7π/5, 9π/5} for five weakly coupled oscillators (N=5) with the frequencies {ω.sub.1,ω.sub.2,ω.sub.3,ω.sub.4,ω.sub.5}={0.98ω, 0.99ω, ω, 1.01ω, 1.02ω}, where ω is the natural frequency. The coupling functions are K.sub.ij(φ)=sin(φ)+cos(φ) with the coupling strength α.sub.ij=0.02 for i, j=1, . . . , 5. The designed interaction function Λ.sub.v.sup.r(φ) and the optimal waveform v.sup.r(θ) are shown in FIG. 14A and FIG. 14B, respectively. The desired desynchronization phase structure stabilizes after 150 input cycles, as shown in FIG. 13C.

Entrainment of a Network of Chemical Oscillators

[0148] A two-cluster entrainment pattern in a network of 15 chemical oscillators was created as shown in FIG. 15A. The frequency distribution of these 15 oscillators is shown in FIG. 15B. The designed optimal waveform v.sup.r(θ) and the desired phase structure are shown in FIG. 15C and FIG. 15D, respectively.

[0149] This written description uses examples to disclose the invention, including the best mode, and also to enable any person skilled in the art to practice the invention, including making and using any devices or systems and performing any incorporated methods. The patentable scope of the invention is defined by the claims, and may include other examples that occur to those skilled in the art. Such other examples are intended to be within the scope of the claims if they have structural elements that do not differ from the literal language of the claims, or if they include equivalent structural elements with insubstantial differences from the literal languages of the claims.