All known closed-loop techniques assume one can monitor the activity in question and stimulate at least a large part of the oscillator population. The population is supposed to be highly interconnected, so the mean-filed approximation is justified. If continuous-time stimulation is possible, one achieves the desired goal by feeding back a delayed or phase-shifted measurement (Rosenblum and Pikovsky, 2004b; Tukhlina et al., 2007). The explanation is simple: elements of the population maintain synchrony due to forcing by the mean field. If stimulation compensates for the mean field, the oscillators naturally desynchronize due to frequency inhomogeneity and individual noisy perturbations. It means the stimulation shall be approximately in antiphase to the measurement and have roughly the same amplitude. Thus, the problem reduces to determining the proper phase shift and amplification; this can be done by trial and error or by an adaptive control scheme (Montaseri et al., 2013). The paramount property of this technique is the stimulation tends to zero, or practically to the level of noise, as soon as the control goal is achieved and the undesired rhythm declines. Therefore, the scheme is denoted as the vanishing stimulation. In a practical situation, one measures a mixture of the activity to be suppressed with other rhythms and noise, so filtration shall be a part of the feedback scheme.

Another way to describe the mechanism of the continuous feedback stimulation is to assume that the collective mode appears via the Hopf bifurcation and write the corresponding normal form equation. Linear delayed or non-delayed feedback (Rosenblum and Pikovsky, 2004b; Tukhlina et al., 2007) changes the linear term and thus stabilizes the otherwise unstable origin. On the level of individual oscillators, this corresponds to the asynchronous state. The normal form representation clarifies the effect of the nonlinear feedback (Popovych et al., 2005): it changes the nonlinear term in the normal form equation, thus reducing the collective mode’s amplitude without causing the bifurcation.

In many applications, especially in neuroscience, continuous-time stimulation is not feasible and has to rely on a pulsatile feedback scheme. The most straightforward solution is to employ the continuous feedback signal to modulate a high-frequency pulse train (Popovych et al., 2017). Another approach is to modulate the pulses’ amplitude following the measurement but stimulate only at the vulnerable phases (Rosenblum, 2020); in this way, one can essentially reduce the intervention in the system. In this paper, we follow this path and simulate the desynchronizing feedback accounting for such practical issues as extracting the signal of interest from its mixture with noise and real-time estimation of the signal’s phase and amplitude. We suggest a network model with time-dependent connectivity to mimic the time course of data registered from Parkinsonian subjects and demonstrate a successful synchrony control of such a network by appropriately timed charge-balanced pulses. Next, we propose and discuss an algorithm for automatically tuning feedback parameters.

Before modeling the pulsatile feedback-based desynchronization, we discuss why certain oscillation phases are vulnerable. Consider a limit-cycle motion subject to an external stimulus that acts along some unknown direction. This direction depends on the system’s equations and how the stimulus enters these equations. Since we do not know the system’s equations, we cannot predict how the given stimulus influences the oscillation amplitude. However, we can say that the stimulus pushes the system off the limit cycle. If we apply the stimulus at such a phase that it pushes the system towards the unstable steady state inside the cycle, this action is optimal for desynchronization. Indeed, the ensemble’s synchronous state corresponds to the collective activity’s limit-cycle oscillation, while the asynchronous state corresponds to a fixed-point solution. Thus, there exists an optimal, vulnerable phase θ ₀. Obviously, stimulation at θ ₀ + π with a stimulus of opposite polarity is also optimal. Unfortunately, we cannot guess θ ₀. We can only apply stimuli at different phases and look at the response.

A commonly used model of collective dynamics in a highly interconnected oscillatory network is a mean-field coupled population. In computational neuroscience, the models frequently incorporate two groups of units, modeling excitatory and inhibitory neurons; in this case, the coupling organization may be more complex. The common feature of the models from this class is that they produce a nearly periodic collective mode. The deviation from the periodicity is due to the finite-size effect and, sometimes, weak collective chaos that is more pronounced if individual units are chaotic. However, the envelope of the rhythm to be suppressed does not vanish, and currently available techniques successfully treat this case. We concentrate here on the case when the activity of our interest is waxing and waning, as frequently observed in neuroscience measurements [for an example, see the time plots of band-passed beta-band activity from a patient with Parkinson’s disease in (Rosenblum et al., 2021)].

We use a phenomenological model to simulate the waxing and waning patterns when the activity bursts are interrupted by quiescent epochs. We take N globally coupled heterogeneous ¹ Bonhoeffer–van der Pol oscillators: x ̇ k = x k − x k 3 / 3 − y k + I k + ε t X + cos ⁡ ψ ⋅ P t , y ̇ k = 0.1 x k − 0.8 y k + 0.7 + sin ⁡ ψ ⋅ P t , (1) where k = 1, …, N is the oscillator index, and the coupling is organized via the mean field X = N ⁻¹ ∑ _k x _k . Function P(t) describes external stimulation, and the parameter ψ quantifies how the stimuli affect the system and determines the phase beneficial for stimulation. For the following, we fix ψ = π/4. In the simulation of the feedback loop, this parameter is considered unknown, and a search algorithm is used to find an appropriate stimulation phase.

A new feature of the model is the time-dependent coupling strength ɛ(t). We assume that it fluctuates around some central value ɛ _c . For simulations, we use a simple algorithm; namely, we let ɛ(t) be a piece-wise constant function of time, ɛ(t) = ɛ _n = const for t _n ≤ t < t _n+1 = t _n + τ _n , where ɛ _n and τ _n are random numbers, chosen from a uniform distribution between ɛ _c ∓Δɛ and τ _min, τ _max, respectively. Figure 1 illustrates the synchronization transition in the ensemble of N = 1,000 elements for Δɛ = 0.015, τ _min = 200 and τ _max = 500. Here, for the order parameter, we take the oscillation amplitude ² , averaged over a long time interval, ⟨a(t)⟩. We see, that for 0.015 ≲ɛ _c ≲ 0.03 we obtain strongly modulated mean field X. The time plots are given in the Results section below. For comparison, we also show the synchronization transition for the usual globally-coupled model with ɛ = const (formally, it corresponds to Δɛ = 0). For simulations, we exploited the fourth-order Runge-Kutta scheme with the time step 0.1.

The suggested model is certainly purely phenomenological. However, it is natural to assume that, e.g., neuronal networks in the brain are not frozen but permanently altered, e.g., due to plasticity (Manos et al., 2021; Asl et al., 2022; Khaledi-Nasab et al., 2022). Real-world networks are nonautonomous, adaptive (Gross and Blasius, 2008; Berner et al., 2023), and subject to internal fluctuations and external perturbations. Though the account of all details is impossible, a phenomenological model can help to describe the dynamics close to the synchronization transition point.

We model stimulation by short stimuli. Thus, P(t) = ∑ _n p(t − t _n ), where t _n are the instants of stimulus’s application and p(t) describes the pulse shape. The stimuli have a finite length T _s , i.e., p(t) = 0 for t∉[0, T _s ]. In many neuroscience applications, an additional requirement applies: the stimuli must be charge-balanced to avoid charge accumulation in the live tissue. Thus, we use bipolar stimuli consisting of two rectangular pulses of opposite polarity and equal area to fulfill the condition ∫ 0 T s p ( t ) d t = 0 . Notice that, generally, there is a gap between two monopolar pulses. We apply them when the collective oscillation phase is close to the optimal one (for the moment, we assume that this optimal phase θ ₀ is known). Namely, we check whether the circular distance between θ and θ ₀ is smaller than tolerance α. We also introduce the minimal interval Δ between the pulses. In dependence on α, Δ, there can be one or several stimuli when phase θ(t) crosses the α-vicinity of θ ₀. Next time, the stimulation is turned on when θ ≈ θ ₀ + π, within the tolerance α. We illustrate the algorithm in Figure 2.

Next, we discuss the intensity of the stimulation. We implement a feedback algorithm with the feedback factor ɛ _fb < 0. It means the signal’s instantaneous amplitude a(t) determines the stimulation amplitude A _n , i.e., the height of the narrow rectangular pulse of the nth stimulus: A _n = A(t _n ) = max(ɛ _fb a(t), − A ₀) for stimulation around θ ₀ and A _n = A(t _n ) = −max(ɛ _fb a(t), − A ₀) for stimulation around θ ₀ + π. ³ The parameter A ₀ determines the maximal allowed stimulation amplitude.

Until now, we assumed that the optimal (vulnerable) phase θ ₀ and the feedback factor ɛ _fb are known. In an actual experiment, we must determine these parameters, i.e., we need an adaptive control scheme. A natural approach is to try stimulation with different θ ₀, ɛ _fb and look at whether the signal’s amplitude is enhanced or suppressed. The known algorithms (Montaseri et al., 2013; Rosenblum, 2020) work well for signals without strong amplitude modulation, both for time-continuous and pulsatile stimulation. However, our tests find that these adaptation schemes become unstable for strong modulation, i.e., in the case treated in this paper.

We solve the problem by exploiting a trial-and-error algorithm. It means we try to stimulate at various phases and look when the stimulation is most efficient. However, this is not that trivial for the case of a burst-like, amplitude-modulated signal because we cannot tell the amplitude reduction due to stimulation from that due to internal dynamics. Indeed, the data’s essential feature is that the amplitude decreases practically to zero from time to time without any stimulation. On the other hand, even if we stimulate at the proper phase, the amplitude can grow due to increased coupling within the oscillator population. We approach the problem as follows.

In the first step, we do not stimulate but calculate the autonomous system’s average amplitude a _aut . Next, in the learning phase, we start with some initial values for ɛ _fb and then gradually change θ ₀ from zero to 2πN _cycl . Namely, we stimulate as described for m = 5 oscillation periods and then check whether a _curr < a _min. Here, a _curr is the current amplitude value, averaged over m periods, and a _min is the minimal amplitude value (Initially, we set a _min = 0.3a _aut .) If a _curr < a _min, we set a _min = a _curr and choose the current value for the optimal phase as θ _opt = θ ₀. Then, we stimulate the next m cycles and check the condition a _curr < a _min again. If a _curr > a _min but decreases, we make no changes and continue stimulation; otherwise we adjust the parameters as θ 0 ↦ θ 0 + max Δ θ , Δ θ ⋅ a curr / a aut , where Δθ is the maximal allowed step, e.g., Δθ = 2π/25, and ε f b ↦ ε f b − f ε f b , (2) where f is a decreasing function of |ɛ _fb |. Using such a function, we avoid feedback that is too strong. For real-time applications, it is beneficial to choose a computationally efficient function; we choose f ( x ) ∼ ( 1 + c x 2 ) − 1 , where c is a parameter. This described step implements the trial-and-error search: if amplitude reduction accompanies the stimulation, we keep the stimulation parameters unchanged. Otherwise, we adjust them. We proceed with the learning epoch unless θ ₀ reaches 2πN _cycl and then take the current values θ _opt and ɛ _fb for the following. (If we never achieved a _min < 0.3a _aut during the learning epoch, we force the adaptation algorithm to perform one more learning cycle.) Thus, the length of the learning epoch is determined by the number N _cycl of complete cycles of variation of θ ₀; Figure 7 in the next Section illustrates the dependence of the results on N _cycl . After the learning epoch, we keep optimal phase constant, continue tracing the signal’s amplitude, and adapt the feedback factor according to Eq. 2 if a _curr > 2a _min.