Feedback control of collective dynamics in an oscillator population with time-dependent connectivity

We present a numerical study of pulsatile feedback-based control of synchrony level in a highly-interconnected oscillatory network. We focus on a nontrivial case when the system is close to the synchronization transition point and exhibits collective rhythm with strong amplitude modulation. We pay special attention to technical but essential steps like causal real-time extraction of the signal of interest from a noisy measurement and estimation of instantaneous phase and amplitude. The feedback loop’s parameters are tuned automatically to suppress synchrony. Though the study is motivated by neuroscience, the results are relevant to controlling oscillatory activity in ensembles of various natures and, thus, to the rapidly developing field of network physiology.


Introduction
The merger of control science and nonlinear dynamics ideas led to essential achievements in controlling chaos, coherence of noisy dynamics, and noise-induced motion; see (Schöll and Schuster, 2008) for a review.Another field of application is controlling the level of synchrony.As is known, synchronization is a general and frequently encountered phenomenon that may be beneficial or harmful (Winfree, 1967;Kuramoto, 1984;Pikovsky et al., 2001;Strogatz, 2003;Osipov et al., 2007).So, synchronizing low-power generators helps to sum their outputs and thus create a high-power source; see, e.g., (Tiberkevich et al., 2009).Synchronization is vital for the stable operation of power grids (Dörfler and Bullo, 2012) and cardiac pacemakers (Jalife, 1984).Examples of an opposite value stem from neuroscience, where pathologically increased synchrony is believed to be a cause of, e.g., Parkinson's disease (Tass, 1999;Tass, 2000) and pedestrian bridge dynamics (Dallard et al., 2001;Eckhardt et al., 2007).This paper addresses the case of a negative impact.It discusses suppressing undesired synchrony in a highly interconnected population of active, self-sustained oscillators -a problem of general interest for network physiology (Ivanov, 2021).
techniques can be categorized as continuous-time and pulsatile methods.Here, we concentrate on the latter case, suitable for neuroscience applications such as deep brain stimulation (DBS) (Benabid et al., 1991;2009;Kühn and Volkmann, 2017), where only stimulation by pulses is possible.
The main features of our modeling study are as follows.1) We suggest a network model with time-dependent global connectivity to imitate brain activity data with strong amplitude modulation.2) We consider that one cannot measure the rhythm of interest directly but shall extract it from its mixture with other rhythms and noise.3) We implement a realistic scheme for real-time estimating the signal's phase and amplitude.4) We suggest an adaptive control approach to tune the stimulation on the fly.Our main result is that we automatically find two vulnerable phases per oscillation period and achieve suppression of the strongly modulated activity by stimulation at these phases only.
In the rest of this section, we briefly review the main results in the feedback-based manipulation of the ensemble synchrony and specify our problem in more detail.

Feedback control of collective synchrony: state of the art
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, continuoustime 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 limitcycle oscillation, while the asynchronous state corresponds to a fixed-point solution.Thus, there exists an optimal, vulnerable phase θ 0 .Obviously, stimulation at θ 0 + π with a stimulus of opposite polarity is also optimal.Unfortunately, we cannot guess θ 0 .We can only apply stimuli at different phases and look at the response.

Modeling the closed-loop desynchronization 2.1 Simulating the amplitudemodulated activity
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 1 Bonhoeffer-van der Pol oscillators: where k = 1, . .., N is the oscillator index, and the coupling is organized via the mean field X = N −1 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 piecewise constant function of time, 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 2 , 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.

Real-time preprocessing 2.2.1 Measurement and filtering
Measuremental noise is inevitable.This noise and some irrelevant rhythms are often intense and contaminate the useful signal, so filtration becomes necessary before phase and amplitude estimation.For real-time applications, the filter can use only the current and previous measurement points, i.e., it must be causal.Thus, we implement a finite impulse response bandpass filter of length 2M + 1 points.The filter length M determines its quality: a longer digital filter provides better signal attenuation for the frequencies outside the desired band.Finite impulse response filters have linear phase response, which means they do not distort the signal but delay it by M points.Thus, the filter's output is delayed by Mδ, where δ is the sampling rate.Hence, we can obtain only delayed values of the instantaneous phase and amplitude.If the delay is small, we can estimate the actual phase as θ(t) ≈ θ(t − Mδ) + ωMδ, where ω is the average oscillation frequency.However, the larger the delay, the worse the estimation because the system is noisy and exhibits phase diffusion.Thus, the filter shall not be too long.On the other hand, short filters are not efficient.Hence, there is an optimal filter semilength M.

Causal phase and amplitude estimation
Here, we implement an algorithm exploiting a non-resonant linear oscillator described in (Rosenblum et al., 2021).The idea is as follows.Suppose, for a moment, our signal is harmonic, and we use it to drive a linear damped oscillator; we use this oscillator as a virtual "measuring device."The phase and amplitude of the oscillator are related to those of the driving force by well-known resonance curve formulas.Since we know the state of the linear oscillator, we invert the resonance relations to obtain the phase and amplitude of the force, i.e., the phase and amplitude of the analyzed signal.The crucial step is to choose the proper parameters for the "device."By choosing the linear oscillator's frequency much larger than the signal's frequency and selecting appropriate damping, we can ensure that the amplitude and frequency response is approximately flat in a large interval of input frequencies.We need two oscillators to achieve this: a strongly-damped one yields a nearly constant amplitude response for frequencies much smaller than the resonant one, while a weakly-damped one provides a nearly constant phase response.The algorithm also works for the force with slowly varying amplitude and frequency.Practically, we substitute each linear oscillator with the corresponding differential equation and solve it numerically for the discretely-spaced input.Reference (Rosenblum et al., 2021) provides an efficient numerical scheme.The approach has demonstrated its efficiency in tests with bandpassed beta activity of Parkinsonian subjects (Busch et al., 2022).

Stimulation
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 Ts 0 p(t)dt 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 θ 0 is known).Namely, we check whether the circular distance between θ and θ 0 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 θ 0 .Next time, the stimulation is turned on when θ ≈ θ 0 + π, 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 0 ) for stimulation around θ 0 and A n = A(t n ) = −max(ε fb a(t), − A 0 ) for stimulation around θ 0 + π. 3 The parameter A 0 determines the maximal allowed stimulation amplitude.
Until now, we assumed that the optimal (vulnerable) phase θ 0 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 θ 0 , ε 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 θ 0 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 = θ 0 .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 where Δθ is the maximal allowed step, e.g., Δθ = 2π/25, and This schematic figure illustrates the algorithm of stimulation at vulnerable phases.Let such phase be θ 0 .It means that a negative pulse applied at this phase pushes the point on the limit cycle (blue closed curve) towards the unstable steady state inside the cycle (filled circle), and so does a positive stimulus applied at θ 0 + π.Practically, we check whether the circular distance between the current phase θ(t) and optimal phase θ 0 is smaller than the tolerance parameter α.So, the system is stimulated if the phase is in the vicinity of θ 0 or θ 0 + π, as indicated by double-headed arrows.Generally, several stimuli (shown by red bars) fulfill this condition; for this illustration, we assume that this number is two.The right part of the figure shows the actual stimulus shape used in simulations.Here, two negative stimuli are shown.Each stimulus consists of a narrow but high rectangular pulse followed after some gap by a compensating pulse that is low but wide.
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 + cx 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 θ 0 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 θ 0 ; 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 .

Results
In this section, we show the results of feedback-based suppression of the noise-contaminated mean field X(t) generated by the globally coupled population described by Eq. 1.For the observable used by the control scheme, we take X N = X + σξ, where ξ is white Gaussian noise with zero average and standard deviation one.We fix σ = 3, i.e., the noise is strong.
As an auxiliary step, we demonstrate in Figure 3 the effect of the digital filter.4Above, we argued that there is an optimal filter semilength M. Here, we show the suppression degree in dependence on M. To concentrate on this issue, we consider the deterministic case of Eq. 1 with ε = const = 0.03.We scan M and θ 0 and show by color coding the suppression factor S = S(θ 0 , M).We define the suppression factor as the ratio of the autonomous system's standard deviation over the stimulated system's standard deviation (the latter is computed after a sufficiently long transient).As expected, short filters are inefficient, while for long filters, the stimulation setup becomes more sensitive to the choice of θ 0 .For the following, we choose M = 350; for this filter, good suppression is achieved for 0.64π ≲θ 0 ≲ 0.84π.
Figure 4 presents the main results.Now, we exploit the model with the time-dependent coupling with parameters ε c = 0.025 and Δε = 0.015. 5We run the unstimulated system for t < 2 • 10 4 , then turn on the stimulation and successfully suppress the collective oscillation.We take N cycl = 1 for the adaptation epoch.The estimated optimal phase θ opt = 0.67π corresponds to the domain of efficient suppression for the chosen filter.
Figure 5 and Figure 6 provide further illustrations for this test.In Figure 5, we show a short epoch of the data from Figure 4 for the suppressed state magnified so that one sees individual stimuli.Figure 6 proves that the decrease of the amplitude is indeed due to stimulation and not internal dynamics.To demonstrate this, we run the system twice, with the same initial conditions and the same realization of the random time-dependent coupling ε(t).However, we switch on the stimulation in the first case while the system remains autonomous in the second.
The natural question is whether the suggested algorithm always finds the proper stimulation parameters.The answer is no.The results in Figure 4 are typical, but sometimes the algorithm fails.Since it traces the amplitude variation in response to stimulation, the errors are inevitable because the amplitude can occasionally reduce due to internal dynamics, not stimulation.To estimate the probability of a failure, we performed 1,000 runs with different random realizations of the time-dependent coupling and computed the suppression factor S. We repeated this test with two and three learning cycles, letting θ 0 in the learning phase grow to 4π and 6π, respectively.The results shown in Figure 7 demonstrate that the stimulation algorithm fails in ≈ 10% of runs (S < 1, the amplitude increases) and provides essential suppression (S > 2) in ≈ 70% of trials.We also see that adaptation with several learning cycles is not helpful: the performance in the case of three cycles is only slightly better than for one cycle.

Discussion
We presented a detailed simulation of the pulsatile feedbackbased control of collective synchrony, concentrating on the case when the observed dynamics exhibit strong amplitude modulation; this case is relevant, e.g., for neuroscience.To simulate such modulated activity, we introduced a model of an oscillatory network with global randomly fluctuating coupling.We assumed Suppression factor (color-coded) as a function of the stimulation phase θ 0 and causal filter semilength M. Short filters are not efficient and, therefore, yield poor suppression.Large M ensures much better filtration but increases the delay introduced by the filter and thus spoils the overall performance of the feedback desynchronization scheme.The optimal M is about 350.For larger M, the maximal suppression factor remains nearly the same, but the algorithm becomes more sensitive to the choice of θ 0 .
that the input to the feedback loop is a mixture of the ensemble mean field and broad-band noise of high intensity.We demonstrated an efficient network desynchronization by incorporating a bandpass filter into the feedback loop and using an adaptation algorithm to tune the stimulation setup.We emphasize that we stimulate at vulnerable phases only and reduce the stimulation amplitude as soon as the goal is achieved, thus minimizing the intervention into the system.This property may be helpful for applications in life sciences and network physiology (Ivanov, 2021).
The results show that contamination of the rhythm in question by other spectrally separated rhythms and noise is not an obstacle.The feedback loop compensates for the delay introduced by a causal bandpass digital filter.Substituting the simple filter used in our simulations with a more advanced one [see, e.g.(Smetanin et al., 2020)] may improve the overall performance.
The algorithm finding the optimal stimulation phase and feedback factor can possibly be improved.One option to enhance its performance is to restart the learning epoch after some time if the suppression factor after the first search for optimal parameters is too small.An alternative approach is to substitute the adaptation algorithm by inference of the amplitude response of the system (Duchet et al., 2020;Cestnik et al., 2022) using specially designed test stimulation; the inferred information can then be used for suppressing stimulation.Illustration of the feedback-based pulsatile desynchronization.The top panel presents the mean field X of system (1) and its mixture with noise, X N .The latter is used as the input to the feedback loop.The second panel from the top presents X f , that is, the filtered series X N , and the stimulation that is turned on at t = 2 • 10 4 .At this instant, the learning epoch begins.During this epoch, the stimulation phase θ 0 grows from zero to 2π.After that, the stimulation continues with the optimal value of θ 0 .The bottom panel presents the time-dependent coupling.Here, we magnify the short epoch of the data from Figure 4 (second panel from the top).In this magnification, one can see individual stimuli.Negative polarity pulses appear near the optimal phase θ opt , while positive pulses appear near θ opt + π.There are one or two pulses in the neighborhood of the vulnerable phase for the chosen parameters.

FIGURE 1
FIGURE 1Synchronization transition in system (1) illustrated by the dependence of the time-averaged amplitude 〈a(t)〉 on the coupling strength.Black dashed line presents the standard case when the coupling is time-independent, i.e., Δε = 0. Blue, red, and magenta solid lines show the averaged, minimal, and maximal amplitude, respectively.Vertical dotted-dashed line indicates the value ε c = 0.025 used in the following; this value corresponds to the strongly modulated mean field with vanishing envelope, see time plots in Figure4.