Controlling collective synchrony in oscillatory ensembles by precisely timed pulses
CControlling collective synchrony in oscillatory ensembles by precisely timedpulses
Michael Rosenblum a) Institute of Physics and Astronomy, University of Potsdam, Karl-Liebknecht-Str. 24/25, 14476 Potsdam-Golm,Germany (Dated: 3 September 2020)
We present an efficient technique for control of synchrony in a globally coupled ensemble by pulsatile action.We assume that we can observe the collective oscillation and can stimulate all elements of the ensemblesimultaneously. We pay special attention to the minimization of intervention into the system. The key idea isto stimulate only at the most sensitive phase. To find this phase we implement an adaptive feedback control.Estimating the instantaneous phase of the collective mode on the fly, we achieve efficient suppression usinga few pulses per oscillatory cycle. We discuss the possible relevance of the results for neuroscience, namelyfor the development of advanced algorithms for deep brain stimulation, a medical technique used to treatParkinson’s disease.
Networks of highly-interconnected oscillatory el-ements are popular models for various systems,either manufactured or natural. It is well-known that, for sufficiently strong interaction, theunits of the network synchronize, and the sys-tem as a whole exhibits a collective rhythm. Fre-quently this rhythm is detrimental and shall besuppressed: the examples include oscillation ofpedestrian bridges and some pathological brainactivity. On the contrary, if the interaction withinthe network is too weak to induce collective os-cillation, enhancement of synchrony may be de-sirable, e.g., to ensure coherent oscillation ofmany low-power sources so that they producea high-power output. These two related prob-lems call for efficient control techniques, and var-ious schemes have been designed for this purpose.Here we elaborate on a special case when the con-trol action shall be pulsatile, which is a commonrequirement for neuroscience applications. Wedevelop a feedback-based adaptive technique thatachieves suppression of undesired collective syn-chrony with only one or two pulses per oscillationcycle. A slightly modified version of this tech-nique enhances collective synchrony if required.We discuss a possible application to a clinicaltechnique, deep brain stimulation, widely used totreat several neurological diseases.
I. INTRODUCTION
The nonlinear science community has paid a lot ofattention to research on large populations of interact-ing self-oscillatory units. Hundreds (if not thousands)of research articles followed the pioneering publications a) Electronic mail: [email protected] on this topic . Many of them exploited the analyti-cally tractable model of globally coupled phase oscilla-tors . Theoretical, numerical, and experimental studiesdescribed and analyzed many interesting phenomena. Anincomplete list includes the emergence of the collectivemode, clustering, quasiperiodic dynamics, appearance ofheteroclinic cycles, and chimera states, see reviews and references therein.The most important and most studied effect is theemergence of the collective oscillation in the populationdue to the synchronization of individual units. Collec-tive synchrony can be important for maintaining high-power output in a population of low-power generatorsand is known to play a significant role in the genera-tion of both vital and pathological biological rhythms.Therefore, control of synchrony, i.e. either suppressionor enhancement of the collective mode, is a challengingproblem. In particular, the suppression task is motivatedby a possible relevance to a widely used clinical proce-dure, deep brain stimulation (DBS). DBS implies high-frequency pulse stimulation of some brain areas and aimsat an improvement of motor symptoms in Parkinsonianpatients as well as in the case of some other patholo-gies . Though the mechanisms of DBS remain in the fo-cus of research in neuroscience , many researchers fromthe nonlinear science community have adopted a work-ing hypothesis that views DBS as a desynchronizationtask . This hypothesis has been exploited in a numberof model studies suggesting open-loop and closed-looptechniques for suppression . In this paper, we followthis line of research and consider both the suppressionand the enhancement task for a globally coupled network.We extend our previous studies on feedback-based con-trol , concentrating on the case of pulsatile stim-ulation. With the goal to minimize the intervention intothe controlled system, we employ precisely timed pulses,applied at a vulnerable phase that is determined on thefly. In this way, we efficiently desynchronize the oscilla-tory activity by a few pulses per oscillatory cycle.The paper is organized as follows. In Section II wepresent the simplest model of globally coupled Bonho- a r X i v : . [ n li n . C D ] S e p effer – van der Pol oscillators and use it to introduceand illustrate the main idea. Here we also discuss howthe phase of the collective oscillation can be obtainedin real-time. In Section III we present the algorithm foradaptive tuning of the feedback parameters and illustrateits performance with the help of the ensemble of chaoticR¨ossler oscillators. Section IV takes into account limita-tions inherent to neuroscience and presents suppressionby the so-called charge-balanced pulses. Section V is de-voted to the enhancement of collective synchrony whileSection VI summarizes and discusses the results. II. PULSES APPLIED AT A VULNERABLE PHASEA. The basic model and the main idea
We introduce the approach using as an example a sim-ple model of N globally coupled Bonhoeffer–van der Poloscillators: (cid:40) ˙ x k = x k − x k / − y k + I k + εX + cos ψ · P ( t ) , ˙ y k = 0 . x k − . y k + 0 .
7) + sin ψ · P ( t ) , (1)where k is the oscillator index, k = 1 , . . . , N , and theterm εX describes the global coupling. Here X is themean field, X = N − (cid:80) k x k , and the coupling coeffi-cient ε explicitly describes the interaction between theelements of the ensemble. The oscillators are not iden-tical: their frequencies are determined by the parameter I k that is Gaussian-distributed with the mean 0 . . P ( t ) is external pulsatile actionapplied to the ensemble; it will be specified below. Fi-nally, the parameter ψ describes how the external pulsesact on the system. This parameter is considered to beunknown, to imitate the uncertainty in stimulation of areal-world system without any knowledge of its model.Figure 1 illustrates the dynamics of the autonomousensemble, P ( t ) = 0; here we plot Y = N − (cid:80) k y k vs. X for ε = 0 .
03 and N = 1000. A symbol at X ≈ − . , Y ≈ .
55 shows the unstable fixed pointof the globally coupled system . In this representation,suppression of the collective oscillation X ( t ) means thatthe system is put into and is kept in a vicinity of the un-stable fixed point. Suppose the applied external pulsesact along a certain direction, indicated by dashed linesin Figure 1. Obviously, the pulses applied to the systemat phase angles close to θ and the pulses of an oppositepolarity applied at approximately θ + π are most effi-cient for reducing the collective oscillation, and, hence,for desynchronization. On the contrary, the oscillationamplitude is much less affected by the pulses appliedaround θ ± π/
2. This qualitative discussion presentsthe main idea of our approach: in order to achieve thecontrol goal with minimal intervention we have to stimu-late only in a small interval around the vulnerable phase θ . For the rest of this Section we assume that θ isknown, while in Section III we drop this assumption andshow how θ can be found. -2 -1 0 1 X Y θ FIG. 1. Qualitative explanation of the approach. Blue solidline is the limit cycle of the collective mode of the sys-tem (1) (the collective oscillation is periodic, up to finite-size fluctuations). Suppression of the collective mode canbe achieved if the applied pulse pushes the system towardsthe unstable fixed point, shown by a small filled circle at X ≈ − . , Y ≈ .
55. The direction of the applied pulsescannot be chosen: it is predetermined by the equations of thesystems and by the way the stimulation enters these equa-tions; this direction is shown by dotted lines. Obviously, theoscillation amplitude is mostly affected by a pulse, appliedwhen the system’s state is close to the phase angle θ or θ + π ,and less affected if the phase angle is close to θ ± π/
2. Thus,efficient suppression can be achieved by a repetitive applica-tion of pulses of certain polarity at about θ and of inversepolarity at θ + π . B. Phase estimation
For efficient stimulation, we have to monitor the in-stantaneous phase of the collective oscillation on the fly,assuming that we observe only a scalar time series. Be-low we suppose that X ( t ) is measured. To this end, wefollow and introduce a “device” consisting of a har-monic linear oscillator and an integrating unit¨ u + α ˙ u + ω u = X ( t ) , (2) µ ˙ d + d = ˙ u . (3)The role of the harmonic oscillator Eq. (2) is twofold.First, it acts as a band-pass filter and extracts the oscil-latory mode of our interest from its mixture with noise.Second, it yields signal ˙ u which phase is close to that ofthe input X ( t ), provided the frequency ω is chosen tobe close to the mean frequency of X ( t ). The integrat-ing unit Eq. (3) provides a signal, shifted by π/ u . It is convenient to introduce two auxiliaryvariables ˆ x = α ˙ u and ˆ y = αω µd ; their amplitudes areclose to that of X ( t ) while their phases are delayed by0 and π/
2, respectively, cf. . Hence, we can estimatethe instantaneous (proto)phase of X ( t ) as θ ( t ) = arctan(ˆ y/ ˆ x ) . (4) t -202 q , Q FIG. 2. Estimation of the phase of the collective mode of theautonomous ensemble Eq. (1) with the help of the “measuringdevice” described by Eqs. (2,3). Red line shows phase θ com-puted according to Eq. (4) while the black bold line showsthe angle variable Θ = arctan[( Y − . / ( X + 0 . θ is obtainedin real time. In the following, we will also need the instantaneous am-plitude a ( t ) = (cid:112) ˆ x + ˆ y . (5)Figure 2 illustrates how the algorithm for phase esti-mation works with the system (1). The parameter valuesused here are: ω = 2 π/ . α = 0 . ω , and µ = 500. C. Timing and strength of stimuli
To determine when and how to stimulate, we trace theinstantaneous phase θ ( t ) and check whether | θ ( t ) − θ | < Θ tol or | θ ( t ) − θ − π | < Θ tol . (6)If one of these conditions is fulfilled at time instant t n then a pulse of a certain strength A n is applied to allelements of the ensemble. Here Θ tol is the tolerance pa-rameter. For a fixed width of stimulation pulses Θ tol determines whether one pulse (if Θ tol is small) or severalpulses (if Θ tol is sufficiently large) are applied around θ or θ + π , respectively. The strength of each pulse, A n ,is limited by the maximal allowed value, | A n | ≤ A , andis determined by the current value of the instantaneousamplitude a ( t n ): A = ± max( ε fb a ( t n ) , − A ) , (7)where positive and negative signs correspond to stimula-tion around θ and θ + π , respectively. Here ε fb < D. A numerical example
For the first illustration of the approach, we considerthe model (1) with N = 1000 and ε = 0 .
03 and try tosuppress the collective oscillation by rectangular pulses Δ δ A n A n +1 t n −1 t n +1 t n Δ FIG. 3. Stimulation by rectangular pulses. Fixed parameters δ and ∆ determine the pulse width and the minimal inter-pulse interval, respectively. The pulse amplitude, A n , variesfrom pulse to pulse, as determined by Eq. (7) according to theinstantaneous amplitude a at time instant t n . Notice that A n can be both positive and negative. t -101 X , P FIG. 4. Suppression of the collective mode in system (1) byrectangular pulses. The stimulation is switched on at t =1000. Pulse width and minimal inter-pulse distance are δ =0 . T ≈ . of the constant width δ and minimal inter-pulse interval∆, see Fig. 3. We set ψ = 0, and stimulate with negativepulses around θ = 0 and with positive pulses around π .Other parameters are ε fb = − .
05, Θ tol = 0 . π , and A = 0 .
2. Figure 4 demonstrates efficient suppression ofthe collective oscillation. For the chosen Θ tol there arethree (sometimes four) stimuli in a bunch around θ or θ + π .Before proceeding with the further details of our ap-proach, we discuss the meaning of the a priori unknownparameter ψ . It describes the distribution of the stimula-tion between the equations and is related to phase shift,inherent to stimulation. The latter also depends on theproperty of individual oscillators and of the coupling be-tween them, see a discussion in and references therein.Thus, ψ is related to θ , though is not exactly equal toit. To illustrate this and to analyse sensitivity of ourtechnique to the choice of θ we compute the suppressioncoefficient S as a function of θ , for ψ = ± π/ S = std( X ) / std( X s ) , where std means standard deviation and X and X s arethe mean fields in the unstimulated and stimulated sys-tem, respectively. The results show that choice of θ iscrucial and therefore we need a technique for tuning θ q S a) -0.05-0.15-0.3 q b) -0.05-0.15-0.3 FIG. 5. Suppression coefficient for system (1) as a function of θ for ψ = π/ ψ = − π/ ε fb = − . , − . , − .
3. Θ tol = 0 . π , A = 0 . δ = 0 .
2, and ∆ = 1. Vertical arrows touching thehorizontal axis indicate the optimal values of θ detected byan automated algorithm described in Section III. as well as the feedback strength ε fb automatically. Thistechnique is presented in the next Section. III. AUTOMATIC TUNING OF SUPPRESSIONPARAMETERS
For a proper tuning of the feedback-based suppressionalgorithm we adapt the approach developed in our pre-vious publication . Namely, we adjust parameters θ , ε fb after each complete cycle, according to the averagedvalue ¯ a of the instantaneous amplitude a ( t ), see Eq. 5. Tobe exact, the latter is averaged over all points within onecycle, except for the interval where the system is stimu-lated (i.e. except for the points where | θ − θ | < Θ tol and | θ − θ − π | < Θ tol ). The update rules are θ → θ + k ¯ a (1 + tanh[ k (¯ a − a stop )] , (8) ε fb → ε fb − k ¯ a/ cosh( k ε fb ) , (9)where k i and a stop are parameters. The initial conditions,if not said otherwise, are θ ( t ) = 0, ε fb ( t ) = 0.An example of suppression with an automated tun-ing of parameters is illustrated in Fig. 6, for ψ = − π/ θ here is θ ≈ .
03, cf.Fig. 5a; the suppression factor is S = 52 .
6. Stimula-tion is turned on smoothly and its onset is followed bya temporal increase of synchrony, because θ is sweptthrough the interval of angles that are beneficial for en-hancement. For ψ = π/ S = 37 . θ ≈ .
56, cf. Fig. 5b. Parameters are k = 0 . k = 500, k = 0 . k = 5. The parame-ter a stop is taken as 20% of the average amplitude of theautonomous system, i.e. before the feedback is turnedon. A. An example: ensemble of chaotic R¨ossler oscillators
With this example we demonstrate that the approachcan be also applied to more complicated models and that -2-1012 X , P a)036 q b)0 2000 4000 6000 8000 t -0.4-0.20 e f b c)-2 0 2 x k -0.500.511.5 y k d) -2 0 2 x k e) FIG. 6. Suppression of the collective mode in system (1) by anadaptive technique. Panel (a) shows the mean field and stimu-lation that is smoothly switched on at t = 1000. Panels (b,c)show the time evolution of two feedback parameters, θ and ε fb that vary unless the mean field X is suppressed. Panels(d,e) show snapshots of the ensemble in the synchronous state(before the feedback is switched on) and after suppression isachieved, respectively. The snapshots demonstrate that inthe desynchronized state the individual units continue to os-cillate, though not coherently. suppression can be achieved with only two pulses peroscillatory cycle. Next, we explore the dependence of theperformance on most important parameters.We consider an ensemble of globally coupled chaoticR¨ossler oscillators: ˙ x k = − ω k y k − z k + εX + cos ψ · P ( t ) , ˙ y k = ω k x k + 0 . y k + sin ψ · P ( t ) , ˙ z k = 0 . z k ( x k − . , (10)where frequencies ω k are Gaussian distributed with themean ω = 1 and standard deviation 0 .
02. Without stim-ulation the system exhibits the Kuramoto synchroniza-tion transition at the critical coupling ε cr ≈ . . For ε > ε cr the mean-field dynamics is nearly periodic, whilefor ε < ε cr one observes small finite-size fluctuations of X .First, in Fig. 7 we illustrate suppression of synchronyin a system of N = 5000 units, with ε = 0 . ψ = π/
4. Parameters of the feedback system are: ω = 1, k = 0 . k = 0 . tol = 0 . π , A = 2, δ =0 .
2, and ∆ = 0 . tol only two pulses per cyclesare applied (as can be seen in Fig. 7) and the adaptivealgorithm converges to θ ≈ .
47 and ε fb ≈ − . Thesuppression coefficient is S = 33 . -5 0 5 X -505 Y FIG. 7. Suppression of the collective mode in the ensembleof R¨ossler oscillators (10). Black solid curve shows a pieceof trajectory of the unforced system in the mean-field coordi-nates X and Y = N − (cid:80) k y k . Dashed curve shows trajectoryof the controlled system. (The trajectory is omitted for smallamplitudes for better visibility.) Black circles (red squares)indicate the points where negative (positive) pulses are ap-plied. w S a) 0.1 0.3 0.5 d b) FIG. 8. Suppression coefficient S in dependence on the fre-quency ω of the linear oscillator Eq. (2) (a) and on the pulsewidth, δ , for k = 5 (circles) and k = 0 . ω in the interval ±
15% of the collectivemode frequency provides a good suppression with S (cid:38) Next, we check the dependence of S on most impor-tant parameters, starting with the frequency of the lin-ear oscillator, ω , see Eq. 2. Figure 8a presents the re-sults. This plot demonstrates that the technique worksfor a rather broad range of ω . This feature is importantfor treatment of real-world systems with drifting averagefrequency. The second test shows how the performancedepends on the pulse width δ and on the parameter k (Fig. 8b). The latter determines saturation level for ε fb ,so that we can expect that the smaller k the larger ε fb and, correspondingly, S . We also expect that broadeningthe pulse increases efficiency of suppression. Figure 8bindicates that this expectation is correct unless the pulsesbecome too wide and do not any more fit the interval ofvulnerable phases.Finally, we check whether the approach works forstrongly coupled R¨ossler ensemble, ε = 0 . δ = 0 . k = 5 the techniquefails, also with A = 4. For δ = 0 . k = 0 .
5, and A = 4 we achieve suppression with S ≈
11. Helpful is also ini-tial increase of the feedback, i.e. taking ε fb ( t ) = − k = 5. (We alsocompare the final values of ε fb : for k = 5 it remains ≈
1; for k = 0 . − . α determines the width ∆ f of the bandpass, ∆ f = α/ π .Thus, if, e.g., the rhythm in question contains frequenciesbetween 10 and 13 Hz, then ∆ f shall be larger than 3Hz, i.e., for this example α (cid:38) . ω . Parameter of theintegrating unit Eq. (3) shall fulfill µ (cid:29)
1; the value µ = 500 ensures correct integration.The stimulation parameter Θ tol determines the num-ber of stimuli in a bunch: the larger Θ tol , the largeris the number of stimuli, and, correspondingly, the sup-pression factor S , see Fig. 10 below. On the other hand,the more stimuli in the bunch, the large is an interven-tion into the system. Hence, the optimal choice of Θ tol depends on whether high values of S or minimal inter-vention are preferred in a particular application. Themaximal stimulation amplitude A in Eq. (7) also de-pends on the application: it shall be sufficiently small toensure non-destructive action on the system.Finally, we discuss parameters k i in Eqs. (8,9). Theproduct k ¯ a determines the adaptation step for the vari-able θ . If this step is too big, then the adaptation rulemay miss the optimal phase. If the step is too small,then the suppression is slow. A reasonable choice is totake k a ¯ a ≈ .
01. Similarly, the product k ¯ a determinesthe adaptation step for the feedback strength ε fb , andthis step also shall be neither too big nor too small. Thevalue k ¯ a ∼ .
01 works well with the tested models. Thechoice of the parameter k (cid:29) k = 500, so that the tanh-function lookslike the step function. IV. CHARGED-BALANCED PULSES
The electrical stimulation of living systems requires aspecial form of pulses. Since the accumulation of electri-cal charge in the cells can be harmful, the pulses must bebipolar and charge-balanced. Figure 9 provides the sim-plest example of such stimuli. In the rest of this Section,we explore desynchronization with charge-balanced stim-ulation. For the test system, we again take the ensembleof globally coupled Bonhoeffer – van der Pol oscillators,see Eqs. (1). If not said otherwise, the parameters arethe same as in Section IIIb.We begin with the stimuli shown in Fig. 9a. We fix δ = 0 . A n /A n, − = −
10, thenthe charge-balance condition yields ∆ = 10 δ . It turnsout that the result of stimulation essentially depends on∆ . If the negative part of the stimulus immediately fol-lows the positive one (or vice versa) then their actionscompensate each other. Indeed, for ∆ = 0 and ∆ = 2 A n +1 δ Δ Δ A n ,− A n t n t n +1 δ Δ Δ Δ A n +1 A n ,− A n t n t n +1 a)b) FIG. 9. Examples of bipolar charge-balanced stimuli. Panel(a) illustrates the simplest considered shape. Here two stim-uli initiated at t n and t n +1 are shown. Each stimulus con-sists of two rectangular pulses of opposite polarity and thecharge-balance requirement means that the blue horizontallystriped area equals the yellow vertically striped one, i.e. A n δ = A n, − ∆ . Panel (b) illustrates a generalization of theshape in (a). Now N b narrow blue (horizontally striped) rect-angular pulses are followed by one yellow (vertically striped)pulse of opposite polarity (the case N b = 2 is shown). Thecharge-balance condition becomes now N b A n δ = A n, − ∆ . there is no suppression, S ≈ for a detailed modelstudy on the suppression efficacy in dependence on thegap ∆ ). However, for ∆ = 6, to be compared withthe average oscillation period T ≈ .
5, the suppressionfactor is S ≈
40. It means that a narrow pulse comes inthe vulnerable phase while the compensating wide pulseappears close to the least sensitive phase. The efficiencyof the suppression can improve if stimuli shown in Fig. 9bare used. We tested stimulation with δ = ∆ = 0 . = 6, and N b = 2 and N b = 3. As expected, themaximal suppression with S ≈
52 is achieved for N b = 3.Moreover, stimulation with N b = 3 and only once per pe-riod also succeeds to suppress the collective oscillation.We summarize our results in the diagram shown inFig. 10. For comparison, we present here both the re-sults for simple rectangular as well as for charge-balancedpulses. In the trials 1-3 we exploit simple rectangularpulses. The tolerance parameter here is Θ tol = 0 . π ,0 . π , and 0 . π , what yields 3, 2, and 1 pulse in a burst,respectively. In cases 4-9 we use Θ tol = 0 . π and dif-ferent values of N b and ∆ , see figure. We notice thatin case 6 one observes waning and waxing patterns. Infact, these patterns can be often obtained if the system isbrought close to the border of suppression, by decreasing ε fb . This result might be interesting for neuroscience ap-plications because some observations indicate that suchregimes correspond to an improvement in the state ofParkinsonian patients . S rectangular charge-balanced
123 456 Δ = 0Δ = 6 w a n i n g a n d w a x i n g Δ = 2 Δ = 6Δ = 6Δ = 6 FIG. 10. Summary of the results for the Bonhoeffer – vander Pol model. The suppression coefficient is shown for 9different cases, presented schematically above the bar chart.Cases 1-3 and 4-9 correspond to stimulation by unipolar andcharge-balanced pulses, respectively. In cases 1-3 there are3, 2, and 1 rectangular pulse in a burst around the vulnera-ble phase. Cases 4-6 correspond to stimulation with charge-balanced pulses with N b = 1 and a different delay between thenarrow pulses and compensating wide one. In case 7 N b = 2,while in cases 8 and 9 N b = 3. Notice that in case 9 thestimulation is applied only once per period. V. ENHANCEMENT OF COLLECTIVE OSCILLATION
The simplest and most reliable way to increase ensem-ble synchrony by pulsatile stimulation is to apply thestimuli periodically, with some frequency ν . This tech-niques is known as injection locking. It works for net-works of periodic or chaotic oscillators, even if they areuncoupled. However, the frequency ν shall be chosenin a proper way and this may be not an easy task ifthe frequency of the collective oscillation is not knownbeforehand and only finite-size fluctuations of the asyn-chronous ensemble are observed, especially in the pres-ence of noise. The dependence of the standard deviation σ of the mean field on the frequency of the drive has atypical resonance-like shape (see the solid line in Fig. 11).Enhancement can be also achieved via the feedbacktechnique with slightly modified update rules: θ → θ + k ( A stop − ¯ A )(1 + tanh[ k ( A sat − ¯ A )] , (11) ε fb → ε fb + k ( A stop − ¯ A/ ) cosh( k ε fb ) , (12)where A sat is the saturation value. Certainly, this ap-proach also has a frequency parameter, namely the fre-quency of the linear oscillator, ω . If the frequency ofthe collective oscillation is not known a priori , ω shallbe guessed. However, the results are not very sensitiveto the choice of ω , as illustrated in Fig. 11 for the model(10) with the sub-threshold coupling ε = 0 .
02. Otherparameters are A = 1, A stop = 5, A sat = 2. The rect-angular pulses ( δ = 0 .
2, ∆ = 1) were used. w ,n s t d ( X ) FIG. 11. Enhancement of the collective mode in the ensem-ble of R¨ossler oscillators with the sub-threshold coupling. Redsolid line shows dependence of σ = std( X ) on the frequency ofthe periodic pulsatile forcing ν , while the symbols illustratethe enhancement via the adaptive feedback-based approachwith the rule Eq. (11,12); here σ is plotted vs. the param-eter ω . Red circles and blue squares correspond to waningand waxing patterns and to stationary chaotic oscillation, re-spectively. The dashed line indicates the level of finite-sizefluctuations in the autonomous system. VI. DISCUSSION AND CONCLUSIONS
In summary, we presented and tested a closed-loopapproach for control of collective activity in a globallycoupled ensemble . The global coupling conjecture de-scribes many natural phenomena, see, e.g., and, inparticular, models pathological brain activity in Parkin-son’s disease . Global coupling is a good approximationfor highly-interconnected networks, e.g., for randomlycoupled neuronal networks . Furthermore, numericalstudies demonstrate that the feedback schemes developedfor the globally coupled ensembles are efficient for con-trolling highly interconnected networks of excitatory andinhibitory neurons, e.g., of STN and GPe cells with arealistic coupling scheme .The main advantage of our approach is that effi-cient control is accomplished by rare precisely timedpulses. So, we have shown that desynchronization canbe achieved and maintained by only one stimulus per os-cillatory cycle. The control parameters – the feedbackcoefficient ε fb and the value of the phase θ when thesystem is most sensitive to stimulation – are adjustedautomatically. An essential feature of the approach isthat the collective oscillation phase is estimated on thefly, using only previous values of the measured signal.The ”device” for phase estimation is simple: it consistsof a linear oscillator and integrator and, therefore, can beeasily implemented either via an electronic circuit or dig-itally. Another essential property of the feedback schemeis that it ensures the vanishing stimulation and main-tains the desynchronized state by small-amplitude stim-uli. Next, the developed technique aims at the desyn-chronization on the collective level while preserving theoscillation of individual units, i.e., the goal of the controlis to avoid oscillation death of ensemble elements, but de-stroy the coherence of their activity. In this respect, ourstudy follows the line of the research reported in . All these studies assume that no access to individual units orconnections is possible, neither for measurement nor forstimulation: it is supposed that only the collective modecan be observed and that control acts on the whole popu-lation, in contradistinction to the techniques like pinningor push-pull control that rely on access to a subset ofthe units.The idea of applying pulses at a vulnerable phase goesback to P.A. Tass publications about twenty years ago.The significant improvement brought by our approachis due to the feedback loop. Closed-loop control auto-matically detects the vulnerable phase and reduces thestimulus amplitude when the suppression is achieved. Asa result, the desired state is maintained by only one ortwo weak stimuli per oscillation period. Last but notleast, our approach does not rely on phase approxima-tion. Hence, it can be applied not only to ensembles ofcoupled limit-cycle oscillators but also to ensembles ofchaotic units, e.g., bursting neuronal models. For theformer case, when the phase sensitivity curve can be in-troduced, one can associate the vulnerable phase withthe phase interval where the slope of the phase sensitiv-ity curve is most steep. In this case, a common stimula-tion acts differently on units with close phases and shiftsthem apart, resulting in desynchronization.Finally, we discuss the possible application of the pro-posed approach to neuroscience. We rely on a quite gen-eral assumption that the rhythms to be controlled emergedue to interaction in a large neuronal population. In thisrespect, we remain in the framework of the working hy-pothesis frequently exploited by the nonlinear commu-nity. Though we tested the approach on rather simplemodels, we believe that it works for more sophisticatedones, as long as the rhythms appear due to synchroniza-tion in a highly-interconnected network. Certainly, thedynamics of the human brain is much more complex, andthe synchronization hypothesis may turn out too simplis-tic. However, our model-based approach can be useful asit is or in combination with ad hoc model-free closed-loop techniques for DBS that are nowadays under devel-opment in the neuroscience community . In this con-text we especially mention the phase-specific stimulationsuggested and implemented in . We believe that incor-porating the stimulus amplitude adaptation mechanismalong with the automated tuning of phase relation shallbe an essential improvement.We emphasize several properties of our technique thatrender it suitable for DBS application. (i) It works withrealistic charge-balanced stimuli. Though we have notsearched for the optimal shape of stimuli, we have shownthat stimulation is efficient if pulses of opposite polarityappear at the most and least sensitive phases, respec-tively, cf. . It means that the condition of charge bal-ance is fulfilled on a time scale of about one-fourth ofthe oscillatory cycle. (ii) Stimulation and measurementare separated in time. Indeed, the adaptive algorithmrelies only on the values of the instantaneous amplitudebetween the epochs where stimulation is applied. (iii)Since the optimal phase for stimulation is determinedautomatically, it does not matter whether the mean field X or its phase-shifted version is measured. This propertyis useful if the stimulation and measurement of brain ac-tivity are performed at different cites. (iv) The linearoscillator used for phase estimation also acts as a band-pass filter and, therefore, extracts the rhythm of interestfrom the raw signal. However, the bandwidth of the fil-ter is quite large – the property required to deal with thesignals with drifting frequency.As a direction for further improvement, we mentionthe modification of the adaptation rule to allow for bothincrease and decrease of θ . This modification will reducethe transient time for desynchronization and will help toavoid a temporal increase of synchrony in the process ofadaptation. ACKNOWLEDGMENTS
The author acknowledges fruitful discussions withA. Pikovsky, S. Yanchuk, L. Feldmann, A. K¨uhn, andW.-J. Neumann.
DATA AVAILABILITY
Data sharing is not applicable to this article as no newdata were created or analyzed in this study. A. T. Winfree, J. Theor. Biol. , 15 (1967); The Geometry ofBiological Time (Springer, Berlin, 1980). Y. Kuramoto, in
International Symposium on MathematicalProblems in Theoretical Physics , edited by H. Araki (SpringerLecture Notes Phys., v. 39, New York, 1975) p. 420;
ChemicalOscillations, Waves and Turbulence (Springer, Berlin, 1984). S. H. Strogatz, Physica D , 1 (2000);
Sync: The EmergingScience of Spontaneous Order (Hyperion, NY, 2003). A. Pikovsky, M. Rosenblum, and J. Kurths,
Synchronization. AUniversal Concept in Nonlinear Sciences. (Cambridge UniversityPress, Cambridge, 2001); A. Pikovsky and M. Rosenblum, Chaos , 097616 (2015). J. A. Acebron, L. L. Bonilla, C. J. P. Vicente, F. Ritort, andR. Spigler, Reviews of Modern Physics , 137 (2005); G. Os-ipov, J. Kurths, and C. Zhou, Synchronization in OscillatoryNetworks (Springer-Verlag, Berlin Heidelberg, 2007); M. Break-spear, S. Heitmann, and A. Daffertshofer, Frontiers in HumanNeuroscience , 190 (2010). A. Benabid, P. Pollak, C. Gervason, D. Hoffmann, D. Gao,M. Hommel, J. Perret, and J. De Rougemont, Lancet , 403(1991); A. Benabid, S. Chabardes, J. Mitrofanis, and P. Pollak,Lancet Neurol. , 67 (2009); A. K¨uhn and J. Volkmann, Mov.Disorders. , 11 (2017). M. D. Johnson, S. Miocinovic, C. C. McIntyre, and J. L. Vitek,Neurotherapeutics , 294 (Apr 2008); V. Gradinaru, M. Mogri,K. R. Thompson, J. M. Henderson, and K. Deisseroth, Science , 354 (2009); J.-M. Deniau, B. Degos, C. Bosch, and N. Mau-rice, European Journal of Neuroscience , 1080 (2010). P. A. Tass,
Phase Resetting in Medicine and Biology. StochasticModelling and Data Analysis. (Springer-Verlag, Berlin, 1999);Prog. Theor. Phys. Suppl. , 301313 (2000); Europhys Lett. , 15 (2001); P. Tass, Europhys. Lett. (2001); Phys. Rev. E , 036226 (2002). M. G. Rosenblum and A. S. Pikovsky, Phys. Rev. Lett. ,114102 (2004); Phys. Rev. E. , 041904 (2004). O. Popovych, C. Hauptmann, and P. A. Tass, Phys. Rev. Lett. , 164102 (2005). N. Tukhlina, M. Rosenblum, A. Pikovsky, and J. Kurths, Phys.Rev. E. , 011019 (2007). C. Hauptmann and P. A. Tass, J Neural Eng. , 016004 (2009);O. V.Popovych and P. A. Tass, Front Hum Neurosci. , 58 (2012). G. Montaseri, M. Javad Yazdanpanah, A. Pikovsky, andM. Rosenblum, Chaos , 033122 (2013). W. Lin, Y. Pu, Y. Guo, and J. Kurths, EPL (Europhysics Let-ters) , 20003 (apr 2013); S. Zhou, P. Ji, Q. Zhou, J. Feng,J. Kurths, and W. Lin, New Journal of Physics , 083004 (aug2017); D. Wilson and J. Moehlis, PLOS Computational Biology , 1 (12 2016); A. Holt, D. Wilson, M. Shinn, J. Moehlis, andT. Netoff, PLoS Comput Biol. , e1005011 (2016). O. Popovych, B. Lysyansky, M. Rosenblum, A. Pikovsky, andP. Tass, PLOS One , e0173363 (2017). D. Krylov, D. Dylov, and M. Rosenblum, Chaos: An Inter-disciplinary Journal of Nonlinear Science , 033126 (2020), https://doi.org/10.1063/1.5128909 . The fixed point can be found by simulating the ensemble for ε = 0. We emphasize that phase angle θ is not the true phase of theself-sustained oscillatory system but only a protophase, see ,but this distinction is not important for our problem. We alsostress, that θ is related to the parameter ψ but is not equal to it,as will be discussed below. A. S. Pikovsky, M. G. Rosenblum, and J. Kurths, EurophysicsLetters , 165 (1996). Here, in order to shorten the transient we took ε fb ( t ) = − . O. Popovych, B. Lysyansky, and P. Tass, Sci. Rep., 1033(2017). G. Tinkhauser, A. Pogosyan, S. Little, M. Beudel, D. Herz,H. Tan, and P. Brown, Brain , 1053 (2017); G. Tinkhauser,F. Torrecillos, Y. Duclos, H. Tan, A. Pogosyan, P. Fischer,R. Carron, M.-L. Welter, C. Karachi, W. Vandenberghe, B. Nut-tin, T. Witjas, J. R´egis, J.-P. Azulay, A. Eusebio, and P. Brown,Neurobiol. Dis. , 217 (2018). The approach certainly applies to a single oscillator as well. Inthis case, the stimulation quenches the self-sustained oscillationby keeping the system in a small vicinity of the unstable fixedpoint. M. G. Rosenblum, N. Tukhlina, A. S. Pikovsky, and L. Cim-poneriu, Int. J. of Bifurcation and Chaos , 1989 (2006). H. Su and X. Wang,
Pinning Control of Complex Networked Sys-tems: Synchronization, Consensus and Flocking of NetworkedSystems Via Pinning (Springer-Verlag, Berlin Heidelberg, 2013);Z. He, X. Wang, G.-Y. Zhang, and M. Zhan, Phys. Rev. E ,012909 (2014). B. Rosin, M. Slovik, R. Mitelman, M. Rivlin-Etzion, S. N. Haber,Z. Israel, E. Vaadia, and H. Bergman, Neuron , 370 (2011);S. Little, A. Pogosyan, S. Neal, Z. B., L. Zrinzo, M. Hariz,T. Foltynie, P. Limousin, K. Ashkan, J. FitzGerald, A. Green,T. Aziz, and P. Brown, Ann Neurol. , 449 (2013). H. Cagnan, J.-S. Brittain, S. Little, T. Foltynie, P. Limousin,L. Zrinzo, M. Hariz, C. Joint, J. Fitzgerald, A. Green, T. Aziz,and P. Brown, Brain , 30623075 (2013); H. Cagnan,D. Pedrosa, S. Little, A. Pogosyan, B. Cheeran, T. Aziz,A. Green, J. Fitzgerald, T. Foltynie, P. Limousin, L. Zrinzo,M. Hariz, K. Friston, T. Denison, and P. Brown, BRAIN , 132145 (2017); A. B. Holt, E. Kormann, A. Gulberti,M. P¨otter-Nerger, C. G. McNamara, H. Cagnan, M. K.Baaske, S. Little, J. A. K¨oppen, C. Buhmann, M. Westphal,C. Gerloff, A. K. Engel, P. Brown, W. Hamel, C. K. Moll, andA. Sharott, Journal of Neuroscience ; C. G. McNa-mara, M. Rothwell, and A. Sharott, bioRxiv(2020), doi:“bibinfodoi 10.1101/2020.05.21.102335, ; B. Duchet, G. Weerasinghe, H. Cagnan, P. Brown, C. Bick, and R. Bogacz,J. Math. Neurosci. , 4 (2020). B. Kralemann, L. Cimponeriu, M. Rosenblum, A. Pikovsky, andR. Mrowka, Phys. Rev. E , 055201 (2007);77