Adaptive elimination of synchronization in coupled oscillator

We present here an adaptive control scheme with a feedback delay to achieve elimination of synchronization in a large population of coupled and synchronized oscillators. We validate the feasibility of this scheme not only in the coupled Kuramoto’s oscillators with a unimodal or bimodal distribution of natural frequency, but also in two representative models of neuronal networks, namely, the FitzHugh–Nagumo spiking oscillators and the Hindmarsh–Rose bursting oscillators. More significantly, we analytically illustrate the feasibility of the proposed scheme with a feedback delay and reveal how the exact topological form of the bimodal natural frequency distribution influences the scheme performance. We anticipate that our developed scheme will deepen the understanding and refinement of those controllers, e.g. techniques of deep brain stimulation, which have been implemented in remedying some synchronization-induced mental disorders including Parkinson disease and epilepsy.


Introduction
Synchronization phenomena are omnipresently observed in nature and man-made systems. The initial observation of synchronization phenomenon was attributed to Huygens [1,2], and analytical investigations on synchronization, starting from the seminal work [3], have been performed systematically for deterministically coupled, stochastically coupled, time-altered coupled, and spatiotemporally configured oscillators [4][5][6][7][8]. Some of the synchronization phenomena are believed to be constructive in realization of specific physical, biological or/and ecological functions; nevertheless, some of them are regarded as fatal factors that lead to disasters or diseases in human beings. Examples abound, including the fiercely wobbly millennium Bridge in the presence of a large amount of consensus pedestrians [9] and mental disorders like Parkinson disease and epilepsy in the presence of massively synchronized bursting neurons in specific areas of the brain [10][11][12][13][14][15].
Naturally, the question arises, 'How to control and even eliminate the noxious synchronization in a large population of coupled oscillators?' In fact, several techniques and control methods for synchronization elimination have been successively developed [10,[16][17][18][19][20][21]. With an appropriate selection of control parameters, synchronization can be eliminated in models of coupled oscillators through feedback controllers. However, due to unavoidable feedback delay in controllers, the feasible parameter regions for realizing synchronization elimination are neither uniformly nor globally distributed but circumscribed in some narrow areas [20]. This may bring difficulties to practical use of these control techniques. Once the controller is switched on, a precise selection of parameters, which requires accurate information on the coupled oscillators, must be processed. However, it is not easy to instantly obtain the information with minimally invasive techniques. Therefore, from a viewpoint of practical use, an adaptive and minimally invasive technique that can automatically adjust the parameters into a feasible region for synchronization elimination is extremely desirable.
As a matter of fact, a large number of adaptive techniques have been developed and validated to be theoretically and practically useful for chaos control, oscillation death, synchronization, and parameter identification in nonlinear dynamics and even network dynamics [22][23][24][25][26][27][28][29][30][31][32]. Some particular adaptive techniques have been successfully applied for synchronization suppression in coupled oscillators [33][34][35]. These developed techniques, though being very numerically effectual, compulsorily require either only instantaneous signal input [34] or both instantaneous and time-delayed signal inputs [35], and even complicated adaptive configurations [34,35]. Only with unavoidable time-delayed signal input, an adaptive feedback control scheme of much simpler configuration and easier implementation for practical use is highly challenging. Note that in terms of dynamical order parameter of coupled oscillators, synchronization and desynchronization indicate two different steady states with different basins of attraction. This therein invites a delicate design of the adaptive scheme with a feedback delay for eliminating synchronization, different from those conventional techniques solely for archiving synchronization.
In addition, the distribution of natural frequencies, which determines the intrinsic property and classification of uncoupled oscillators, can influence the dynamics of coupled oscillators. There have been a series of works reporting on how the topological properties of the symmetric unimodal and bimodal distribution with or without time delays produce synchronization and even new dynamics [36][37][38][39][40][41][42][43][44][45]. However, in the literature, there have been only a few investigations on how the asymmetry, as well as the shift distance of the bimodal distribution of the natural frequencies, influences the elimination of synchronization in terms of the adaptive scheme with a feedback delay.
In this paper, we present a delicately designed adaptive scheme with a feedback delay for achieving synchronization elimination not only in the coupled Kuramoto oscillators, but also in the representative analogous models of neuronal networks including the coupled FitzHugh-Nagumo oscillators with spiking dynamics and the coupled Hindmarsh-Rose oscillators with bursting dynamics. Moreover, we illustrate why the designed scheme with a feedback delay is efficient, and find out how the exact topological form of the bimodal natural frequency distribution influences the scheme performance. Our findings are useful to make a clear description of the directions for further improvements and real applications of the designed scheme.

Model description: from the coupled normal forms to kuramoto oscillators
To begin with, we consider the coupled Hopf-bifurcation normal form: where z j is a complex-valued variable describing the dynamics of the jth oscillator, w j represents its natural frequency, K is the internal and global coupling strength, and = å = Z z N j N j 1 1 is the mean-field coupling term. All the natural frequencies are not identical but are supposed to obey some probability distribution. Without loss of generality, the probability distribution g(w) for the natural frequency is assumed bimodal as: Here, h 1,2 represent the widths of the two peaks of the distribution g(w), and a stands for the half distance between the two peaks. Particularly, the distribution becomes unimodal when a=0, which makes the two peaks coincide with each other. For simplicity, we assume that + = h h 2 1 2 in (2), and that = + D h 1 , the distribution is symmetric; otherwise, it becomes asymmetric, as shown in figure 1.
Without coupling (i.e. K=0), each individual oscillator rotates along with its own natural frequency in the complex plane, showing unsynchronized dynamics. With a sufficiently strong coupling K, the coupled Hopfbifurcation normal form shows some phenomena of phase synchronization, which will be illustrated in detail later. Here, the mission is to use an appropriate controller to eliminate phase synchronization. Thus, a feedback controller of the mean field is added to the coupled normal form (1), which yields a controlled system: and τ is a feedback delay, which is unavoidable in practical instruments of control. The time-altered variable L(t) is a complex-valued adaptive coupling strength, which needs to be designed to be as simple as possible for practical use. However, those control schemes proposed in the literature for updating the coupling strength are complicated and even require an instantaneous signal in addition to the time-delayed signal [34,35]. Thus, an adaptive feedback scheme, different from and simpler than the exiting methods [34,35], for updating L(t) is pending for delicate design.
To this end, inspired by [16,20], we investigate the phase dynamics of the controlled system (3) by using the i j is the order parameter and 'Im' denotes the imaginary part of a complex number.
which, by using the trigonometric formula, can be further rewritten as Here, both the adaptive coupling strengths g g g pp are time-altered, real-valued variables pending for design. Clearly, without the control term (i.e. º º C S 0), system(5) becomes the classic Kuramoto model with a bimodal natural-frequency distribution g(w). For the symmetric distribution, there exists a critical value K c such that phase synchronization could be always achieved if > K K c , as reported in [46] and references therein. For the asymmetric distribution, through analogous arguments, we still can verify phase synchronization to be achieved with a sufficiently large K. Therefore, in the following discussion, we always suppose that K is large enough to ensure the existence of phase synchronization in the uncontrolled system.

Adaptive techniques for coupled oscillators
With the preparations in section 2, we are able to develop an adaptive technique for the controlled system (5).
, which can be viewed as an energy function for the dynamics of the order parameter. In fact, the synchronization state and the desynchronization state correspond, respectively, to two steady states º ( ) r t 1 and º ( ) r t 0 for the dynamics r(t). Intuitively, the coexistent steady states are located, respectively, at the minimal and maximal energies of t ( ) Q t . In order to realize the synchronization elimination, we use the steepest descent method (SDM) [23,28,47,48], where the desynchronized steady state º ( ) r t 0, the minimal energy, is set as the unique control target. As such, starting at some position, the controlled trajectory, which moves along the downhill gradient of the energy potential, acquires the speediest convergence to the minimal energy. Accordingly, we design an adaptive scheme for variable S and C along the downhill gradient as follows: where ∂ represents a derivative operator on a given function with respect to the subscript argument, h 1,2 and ò are adjustable control parameters with h > 0 1,2 and  < < 0 1, and (·) H is the standard Heaviside step function. Combining the order parameter r and system (5) in the scheme(6) yields a more explicit adaptive scheme: where 'Re' and ' * ' denote the real part and the complex conjugate of a given complex number, respectively, and = å q = ( ) s t e N j N 1 1 2i j . As numerically shown in figure 2(a), we use the adaptive scheme (7) to achieve the synchronization elimination for the coupled dynamics (5). Also shown in figure 2(b) are not only the dynamics of the order parameter r for the entire coupled oscillators, but also the dynamics of the order parameters for the partial Figure 1. From symmetric probability distributions to asymmetric distributions for natural frequency. Here, the parameters a and Δ are set, respectively, in the legend for different distributions.
groups of coupled oscillators corresponding to both sides of the axis of the natural frequency distribution. This manifests that the synchronization elimination is realized not only in the subgroups of oscillators but also among all oscillators in the entire network. As for the Heaviside step function used in the schemes (6) or (7), it is not compulsorily required according to the classical SDM. The necessity of an inclusion of it into (6) or (7) will be illustrated in the final section.
Since we use r  1 j for obtaining the coupled phase dynamics above, keeping the polar distance variable in the mean field Z and leaving out the quantity s(t) of high order (as illustrated in figure 2(c)) yields the adaptive scheme for the original controlled system (3) as follows: i . In fact, we are able to obtain this adaptive scheme (8) directly by introducing an energy function = | | Q Z Z 2 instead of Q as defined above and using the SDM. The detailed computational procedures for getting (7) and (8) are given in appendix A. Meanwhile, we use (8) to successfully eliminate synchronization in the coupled normal form in(3), as numerically shown in figure 2(d). This manifests that the designed adaptive scheme is not only suitable for synchronization elimination in the coupled Kuramoto's oscillators but also for the original oscillator models. This also makes it feasible to use the current scheme to manipulate the synchronization dynamics in the analogous models of neuronal networks.
Furthermore, we are interested in the transition time T from synchronization to desynchronization with respect to the parameters a and Δ. Here, the transition time T is defined as on with a small constant δ, a sufficiently large constant T , and a time instant t on beyond which the control is switched on. As shown in figure 3, T shows a more stable continuum when the distance between the two peaks of the distribution is closer and the asymmetry is not very strong (i.e. a and D | | are relatively small). However, when the distance between the two peaks becomes relatively large (such as > a 1.5), T reveals a rapid fluctuation when increasing the asymmetry degree D | |. This also reveals (b) dynamics of the order parameter for all the above coupled oscillators as well as for the two groups of these oscillators corresponding to the both sides of the axis of the natural frequency distribution (a=2 and D = 0) depicted in figure 1; (c) dynamics s rapidly becomes a minor quantity of high order when the adaptive scheme (7) is switched on; (d) dynamics of the 100 controlled oscillators described by the normal forms (3) before and after the adaptive control (8) is switched on, where the parameters are set as: . The color in the color bar represents the dynamical state of each oscillating neuron. that the asymmetric degree of the natural frequency distribution slows down the transition from synchronization to desynchronization. This slowing-down procedure becomes more outstanding, when the two peaks of the distribution are at a farther distance.

Synchronization elimination in coupled analogous models of neuronal networks
We have shown the efficiency of our developed adaptive scheme for realizing synchronization elimination in the coupled oscillators. The analogous technique enables us to adaptively eliminate the synchronization of spiking or bursting dynamics in models of neuronal networks. Here, we consider two concrete coupled neuronal models, i.e., the FitzHugh-Nagumo model [49,50] and the Hindmarsh-Rose model [51]. First, we analyze the coupled FitzHugh-Nagumo models with an external controller as follows: 100, the internal noise I j is supposed to obey the Gaussian normal distribution ( ) 1, 0.01 , and the external noise All the neurons are divided into two groups: e j in (9) for the first half independently obey the uniform distribution ( ) 0.01, 0.03 , while the others independently obey ( ) 0.04, 0.06 . The settings imply that the two groups of neurons behave within two different bands of natural frequencies as without control each e j determines the natural frequency of the neurons.
In addition, the mean field for the membrane potentials of the oscillators is taken as , and analogous to (8), the adaptive scheme for L(t) in (9) is designed as As shown in figure 4, in the absence of adaptive control , the coupled FitzHugh-Nagumo models in (9) display spiking synchronization dynamics; however, when the control with h = 1 and  = 0.3 is switched on (  t t on ), the synchronization is adaptively eliminated although each individual oscillator still exhibits spiking dynamics.
Secondly, we consider a small-world network of coupled Hindmarsh-Rose neuronal models with an adaptive controller as follows: 100. The parameters are taken as: 3, 0.01 , j j satisfies the uniform distribution ( ) 0.006, 0.007 , and the external noise x ( ) t and its intensity are defined in the same manner as those in the preceding example. The internal coupling is set as , where s = 0.6 and all g kj correspond to a small-world network with N=100 (the node number), = K 2 4 (the number of nearest neighboring nodes before Here, each color represents the argument value of a complex number obtained from the membrane potential signal after the standard Hilbert transformation [52]. (b) Dynamics of the order parameter for the entire coupled oscillators as well as for the two groups of oscillators corresponding to two peaks of natural frequency distribution. (c) Selected are the membrane potential dynamics for the first and the second FitzHugh-Nagumo oscillators from the entire network, which shows the synchronized dynamics and desynchronized dynamics, respectively, before and after the control switched-on time t on . reconnection), and p=0.1 (the reconnection probability). In addition, the adaptive controller t 0, on ], the coupled Hindmarsh-Rose models in (11) display synchronization of bursting dynamics, as shown in figure 5. Once the adaptive control (12) is switched on at = t t on , the synchronized dynamics become eliminated adaptively; however, the bursting dynamics for each individual oscillator are not completely destroyed and can be still observed. This is important since preserving these dynamics is essential to sustainment of the information processing function for each neuron.
Remarks. Actually, the analytical results established in the preceding sections deal with the fully connected networks; however, the above example shows that our adaptive control scheme is also practically feasible for other network patterns, such as consumerresource, mutualism, and scale-free interactions. The analytical and numerical results under the configurations of other network patterns will be discussed elsewhere. Moreover, since we design the adaptive control scheme for the general oscillating models around their limit cycles, the control scheme is useful for the oscillations induced by limit cycles of the Hopf bifurcation. Although a large number of synchronized oscillators are induced by limit cycles in real-world systems, there are still many other types of oscillations [53,54]. Eliminating synchronization of all these oscillators could be one of our future research topics.

Feasible parameter regions
In order to analytically illustrate the feasibility of the adaptive technique proposed above, we set as a constant complex number in (4), so that we need to locate the parameter region of L in the complex plane for eliminating synchronization in (4). To this end, letting  ¥ N and using the Ott-Antonsen ansatz [46] yield:    is stable so that elimination of synchronization is achieved. First, we set D = 0, which corresponds to the symmetric probability distribution defined in (2) for the natural frequencies. In this situation, the characteristic equation obtained above is conjugately invariant, so it is sufficient to assume g p Î [ ] 0, . Without time delay (i.e. t = 0) and with a real-valued L (i.e. g p = 0 or ), the condition < - , 2 1 2 ensures the stability. With a time delay, to find the stability region of the complex-valued L, we investigate the marginal stability condition through substituting l b = i in (14), Through a calculation by using the trigonometric formulas, we further get a polynomial equation with respect to the variable β as which yields four possible roots: . Hereafter, for the parameter a, the half distance between the two peaks of the natural frequency distribution, we separate our discussions into two cases: 1': when the internal coupling is sufficiently strong (i.e., > + ( ) K a 2 1 2 ) , synchronization occurrs in the Kuramoto model without external control. Then, as  (15). This implies that there is no further stability region other than the regions obtained above. When  g p 2, there is no stability region in the | | L -τ plane for any τ.
The above discussions also reveal the stability regions in the complex plane with respect to = + L C S i for given a and K. When t = 0, the stability region is unbounded and located on the left side of two intersecting and endless curves. When t > 0, the region becomes a bounded and leaf-like area. In particular, it shrinks as τ increases, and disappears as τ increases above t = -- . This changing process of the stability region with increasing τ is depicted in figure 7(a) for particularly given a and K.
Case (ii) ' > a 1': when > + ( ) K a 2 1 2 , the changing process of the stability region is the same as the result concluded for case(i) shown in figure 7(a). When  < + ( ) K a 4 2 1 2 , the situation for stability regions becomes more complicated. In fact, from equation (15), it follows that Analogously, for simplicity, denote this piecewise-defined τ by b ( ) h L K n , , , . Through further calculations, we obtain the marginal curves for the stability regions in the | | L -τ plane: . In particular, when g p > 2, the stability region is located below t n 3, and above t + n 2, 1 for some smaller n (see the middle panel in figure 6); when  g p 2, it is located below t n 1, and above t n 4, for some smaller n (see the right panel in figure 6). The above discussions further yield a changeable stability region in the complex plane of L with increasing τ but for given a and K. Concretely, there exist two sequences t 1 , the stability region shrinks from an unbounded area (for t = 0) to a leaf-like area (for t t < < 0 1 ),   For a and k given in (b), the area of the stability region with respect to the increasing τ (c). Here, the area takes positive sign when the stability region is located on the left half of the complex plane, but takes negative sign when it appears on the right half of the plane. changing in the same manner as the stability region obtained for case (i). For t t t Î [ ] , 1 2 , the stability region vanishes. For t t t Î ( ) , 2 3 , the stability region appears in a leaf-like shape again but on the right half of the complex plane. Particularly, there exists a number t M 23 such that for t t t Î ( ] , M 2 23 , the region area is augmented; but for t t t Î [ ) , M 23 3 , it becomes reduced monotonically. For t t t Î ( ) , 3 4 the stability region vanishes again. For t t t Î ( ) , 4 5 , the stability region appears on the left hand of the complex plane again. Similarly, there exists a number t M 45 such that the stability region grows for t t t 45 5 . The left-like region appears back and forth on both half planes until τ increasingly approaches some finite critical value t c2 . As t t > c2 , the stability region no longer appears. With increasing τ, the above-described changing process of the stability region is depicted numerically in the complex plane of L (see figure 7(b) with particularly given parameters a and K ). Still with these parameters, the area of the stability region non-monotonically attenuates with increasing τ, and it eventually vanishes for t t > > 7.2 c2 , as shown in figure 7(c). Analogously, for simplicity, we denote this piecewise-defined τ by b ( ) h L K n , , , . Through further calculations, we obtain the marginal curves for the stability regions in the | | L -τ plane: . In particular, when  g p 2, the stability region is located below t n 1, and above t n 4, for some smaller n; when g p > 2, it is located below t n 3, and above t + n 2, 1 for some smaller n. These discussions further yield a changeable stability region in the complex plane of L with increasing τ but for given a and K. Concretely, there exist two sequences t { } n and t + { } ( ) n n M 1 for some smaller n. For t t Î [ ) 0, 1 , the stability region shrinks from an unbounded area (for t = 0) to a leaf-like area (for t t < < 0 1 ), changing in the same manner as the stability region obtained for Case (i). For t t t Î [ ] , 1 2 , the stability region vanishes. For t t t Î ( ) , 2 3 , the stability region appears in a leaf-like shape again but on the right half of the complex plane. Particularly, there exists a number t M 23 such that for t t t Î ( ] , M 2 23 , the region area is augmented; but for t t t Î [ ) , M 23 3 , it becomes reduced monotonically. For t t t Î ( ) , 3 4 the stability region vanishes again. For t t t Î ( ) , 4 5 , the stability region appears on the left hand of the complex plane again. Similarly, there exists a number t M 45 such that the stability region gets growing for t t t 45 5 . The left-like region appears back and forth on both half planes until τ increasingly approaches some finite critical value t c2 . As t t > c2 , the stability region no longer appears. The analytical validation of the existence t c2 is provided in appendix C. With increasing τ, the above-described changing process of the stability region is depicted numerically in the complex plane of L (see figure 7(b) with particularly given parameters a and K ). Still with these parameters, the area of the stability region non-monotonically attenuates with increasing τ, and it eventually becomes zero for t t > > 7.2 c2 (figure 7(c)). The above discussions focus on the case where the symmetric axis for the natural frequency distribution (2) is centered at w=0. However, when it is not centered at zero but at = ¹ w w 0 0 , the locations of these leaf-like stability regions found above change. As a matter of fact, a utilization of the transformation f q =w t i i 0 to system (5) yields: 0 whose natural frequency distribution indeed has a symmetric axis centered at zero. Therefore, we conclude that for realizing the elimination of synchronization, the complex-valued = , , and τ. Note that complex-valued L and L w 0 have the same norms but different arguments. The stability region of L remains in a leaf-like shape but rotates counterclockwise with the increase of w 0 . Moreover, this rotation period can be directly computed as p t 2 . Figure 8 shows these counterclockwise rotating and leaf-like regions, respectively, for a uniform group of a and K but different τ. These parameters correspond to the case shown in figure 7(b). Clearly, to realize elimination of synchronization, the external coupling strength L has to be selected from some specific stability region, and thus it has some specific forms, such as the forms of only negative real part, only positive real part, and non-real value. For given a K , , and τ, an appropriate form of L depends crucially on the choice of the symmetry axis w 0 , which is an interesting phenomenon akin to the results reported in [55].
Finally, we consider the situation of D ¹ 0, which corresponds to an asymmetric distribution g(w) as defined in (2) for the natural frequency. The absolute value of Δ measures the asymmetric degree of the distribution. Figure 1 shows that the larger the asymmetric degree is, the more asymmetric the distribution becomes. Analytical analysis of the (14) is mathematically feasible but computationally tedious. Here, we numerically show how the stability region of L is influenced by the asymmetric parameter Δ as well as by the time delay τ for a particularly given a and K. On one hand, as shown in figure 9(a), for each given small t < 0.4, the stability region area approaches the maximum at D = 0 because of the convexity of those contours. This implies that the asymmetric degree can reduce the size of the stability region for the coupling strength L. As more explicitly shown in figure 9(b), a larger asymmetric degree also leads to not only a deformation but also spin of the stability region. The real axis no longer coincides with the line crossing both vertexes on the border of the leaf-like region. On the other hand, for some fixed asymmetric degree Δ, the area of the stability region monotonically decreases with the increase of τ (figure 9(a)). A larger asymmetric degree results in a more rapid attenuation of the stability region. Moreover, the change of the area is quite smooth for a fixed t < 0.4 and with increasing D | |, which is different from the rapid change of the transition time for a fixed a with increasing D | | (see figure 3).
Additionally, as shown in figure 9(c), for larger τ, the stability region area shows a symmetry property with respect to Δ, which is analogous to the case in figure 9(a). However, one difference is that with continuously increasing τ above 0.4, the stability region disappears suddenly and then appears rapidly. The differences also include the non-monotonic area of the stability region with increasing τ, which contains a particular case for the stability area shown in figure 7(c) where D = 0. Indeed, for some fixed τ (e.g. t < < 1 1.9 in figure 9(c)), the stability region area drops down dramatically as the asymmetry parameter Δ passes through some critical value depending on τ. This is the reason why the tongue-like region can be observed in figure 9(c).
The above discussions on the stability regions are completely based on the assumption that L is a constant complex number. However, the adaptive scheme (7) we developed and validated in the preceding sections requires the complex-valued L to be always changing and updated as time passes. Although the initial value for L is not necessarily located inside the stability region, the designed adaptive scheme adjusts L at every time instant and directs it to the final destination, the leaf-like stability region. As shown by some panels in figures 8 and 9, the orbits of L, starting from the locations highlighted by the pentagram signs, automatically find their routes to enter the stability regions. Once they approach the regions, elimination of synchronization is achieved definitely and the updating rate for L tends towards zero according to the adaptive scheme (7). Therefore, all the above arguments and discussions theoretically illustrate the reasonability and feasiblity of the adaptive scheme we propose in section 3. 6. Discussions and concluding remarks 6.1. Necessity of Heaviside step function As designed in the scheme (6) or (7), the Heaviside step function is included. Although this term is not compulsorily required from a theoretical viewpoint of the SDM, the feasible stability region of the constant coupling strength º = + ( ) L t L C S i for elimination of synchronization is very limited and even distributes separately in the complex plane, as has been expatiated in the preceding section 5. Once the coupling strength = + ( ) ( ) ( ) L t C t S t i updates its value into this isolated stability region, the update rate, theoretically, should become tremendously small and approaches zero eventually; nevertheless, due to computational errors, they accumulate errors to some extent, which makes it likely to push them out of the stability region. Since the SDM provides only one direction for updating the strengths, departure of the stability region always results in an everlasting missing of this region, which eventually results in a failure of elimination of synchronization. As shown in figure 10, once the control is switched on, but in the absence of the Heaviside step function in the adaptive scheme (7), the failure of elimination of synchronization happens eventually, although at the beginning of the control being switched on, the elimination of synchronization is achieved rapidly and lasts for a certain

(c)
A tongue-like stability region with particularly alterable Δ and τ. Here, a and K are still the same as those in (a). The color values in (b) are defined in the same manner as those in figure 8, and each color value in (c) represents the stability region area after the logarithm, so that the black color indicates a disappearance of the stability region. course of time. Therefore, it is reasonable and of high practical significance to include the Heaviside step function as the manner designed in scheme (6) or (7) for obtaining a persistent and robust elimination of synchronization.

Partly elimination of synchronization
We additionally consider coupled FitzHugh-Nagumo models in (9) as well as their adaptive control (10) with different parameter configurations: e e = º = g 0.1, 0.1 j , and the internal noise I j for the neurons from j=1-50 independently obey the uniform Gaussian normal distribution ( ) 0.7, 0.01 , while I j ( =  j 51, , 100) independently obey ( ) 1.3, 0.01 . With these configurations and without adaptive control, the coupled FitzHugh-Nagumo models in (9) display synchronized spiking dynamics separated into two clusters, which can be seen from the two types of the strip patterns in figure 11 for Î [ ) t t 0, on . Clearly, it is the two types of the internal noises that produce some groups of the neurons as well as their cluster synchronization. Now, when the adaptive control (10) is switched on at = t t on , the order parameter | | r for the whole oscillators rapidly transits to a low-level with a small quantity of fluctuations. However, the order parameters | | r 1,2 , corresponding respectively to the two types of the neurons, fluctuate within a relatively large range. This implies that the synchronization of neither group is not completely but partly eliminated. Also it is the anti-phase-like dynamics of r 1,2 that result in the low-level and mild fluctuations of the order parameter | | r . Here, the performance of the adaptive control is unlike the effectiveness shown in section 4, where the two groups of the neurons are attributed to the two groups of e j . As a mater of fact, a change of either I j or e j leads to a change of the frequency of the oscillators. However, because of the specific forms of the FitzHugh-Nagumo models, a change of I j also alters the location of the equilibrium circled by the oscillating limit cycle, but yet a change of e j does not. Since we use an approximation r  1 j during our design of the adaptive control scheme in section 3, this assumption actually requires an unchanging or slightly changing location of the equilibrium circled by the oscillating limit cycle, which therefore is essential to maximize the efficiency of our designed adaptive scheme. Figure 10. Dynamics of the 500 coupled Kuramoto oscillators (a) and their order parameters (b) when the Heaviside step function is not included in the adaptive scheme (7). All the parameters and control on-off configurations are the same as those in figure 2(a). 6.3. Adaptive scheme for time delay As discussed in section 5, the stability region for L always vanishes occasionally for a particular feedback delay τ and disappears forever when τ is sufficiently large. In real applications, one way to make the designed scheme efficient is to reduce τ as much as possible. However, τ is unavoidable because it does take a certain amount of time to transmit a signal in a remote distance and process a signal to some required mode. Then, the other way is to design an adaptive scheme also for τ. Thus, τ, as a variable, updates its value based on the mean field signal experimentally recorded and finds its way automatically to the nearest and correct delay which admits a stability region for L. Analogous updating schemes for the feedback delay have been investigated, respectively, in [56][57][58].
However, for the current scheme, there is a ceiling for the feedback delay, such as the critical values t c c 1, 2 derived in section 5. Once τ updates its value close to the ceiling, a replacement of it with an appropriately small value that is acceptable in a real system should be conducted. This adaptive updating procedure is repeated until the synchronization elimination is achieved.

Concluding remarks
In this paper, we have designed an adaptive control scheme with a feedback delay to achieve elimination synchronization not only in coupled Kuramoto's oscillators, but also in two representative models of neuronal networks. More importantly, we have analytically and numerically investigated how the topological structure of the bimodal natural frequency distribution influences the elimination of synchronization with a feedback delay. These investigations also support the feasibility of the designed adaptive scheme in synchronization elimination.
In addition, in the treatment of some synchronization-induced mental disorders, techniques, such as deep brain stimulation (DBS), have been validated to have some therapeutic effect; however, the mechanism of this effect is still unknown [12,59]. Major DBS is generally designed as an open-loop control, whose robustness against noise is notoriously weak and whose effectiveness is not persistently lasting [60]. Actually, our developed adaptive scheme is of a closed-loop with automatic control, which usually has several advantages including strong robustness, sustainable effect, and less requirement for manual intervention [14,61,62]. So, the designed scheme and our new findings are expected to be used clinically to deepen understanding and even refine the existing DBS techniques. The clinical application as well as improvement of the current scheme becomes one of the major research directions for our present and future works.  . Finally, a substitution of all these results into the adaptive scheme designed in(6) yields the explicit form of the adaptive scheme (7).
The other adaptive control schemes developed in this work for the coupled normal form models and for the coupled analogous models of neuronal networks can be obtained in an analogous computational argument.