The human body can be viewed as a complex network with various interacting physiological systems. Stimulation of one system might have a strong and potentially even delayed impact on others (Bashan et al., 2012; Ivanov et al., 2016). A deeper understanding of the complex interactions between physiological systems through various signaling pathways and how they lead to the emergence of new physiological states might lead to the development of novel treatments for several diseases (Bartsch et al., 2015; Ivanov et al., 2016). For instance, a recent study on vibrotactile fingertip coordinated reset (CR) stimulation for the therapy of Parkinson’s disease (PD) revealed clinically and statistically significant improvement of motor control along with a significant decrease of abnormal, PD-related high-beta power in the sensorimotor cortex after 3 months of treatment (Pfeifer et al., 2021). Notably, sensory stimuli delivered to only a small part of the body were able to cause pronounced bilateral motor improvement of the entire body (Pfeifer et al., 2021).

PD patients typically suffer from pronounced motor symptoms, as well as several non-motor symptoms such as sensory impairments (Conte et al., 2013), impairments of sensorimotor integration (Abbruzzese and Berardelli, 2003), mood disorder, sleep disorder, and cognitive decline (Weiner et al., 2001; Reich et al., 2022). A network of strongly interconnected brain areas is involved, including the basal ganglia, the thalamus, and the sensorimotor cortex (McGregor and Nelson, 2019). Individual symptoms are associated with abnormally strong neuronal synchrony in these brain areas (Nini et al., 1995; Hammond et al., 2007).

An established treatment for medically refractory PD is high-frequency deep brain stimulation (HF DBS). Chronic delivery of HF DBS to target regions such as the subthalamic nucleus (STN) may suppress PD symptoms during stimulation (Krack et al., 2003; Krauss et al., 2021), where symptoms return shortly after cessation of stimulation (Temperli et al., 2003). Hence, permanent stimulation is required for persistent symptom relief, limiting battery life and, more importantly, increasing the risk of unwanted side effects, e.g., caused by stimulation of neighboring tissue due to current spread as well as stimulation of the very target region (Krack et al., 2002).

To substantially reduce the integral stimulation current, based on computational studies, it was suggested to develop DBS approaches that aim at specifically counteracting pathological synchrony by desynchronzation (Tass, 1999; Tass, 2000). Theory-based desynchronization techniques harness the nonlinear response of ensembles of coupled oscillators to external stimuli. For instance, an ensemble of synchronized oscillators can be desynchronized by delivering a single stimulus pulse at a vulnerable phase of the collective rhythm (Mines, 1914; Winfree, 1977; Warman and Durand, 1989; Tass, 1999). To reliably desynchronize an ensemble of coupled oscillators, it was suggested to deliver such a desynchronizing pulse shortly after a strong phase-resetting pulse (Tass, 2001; Tass, 2002; Zhai et al., 2005). A stimulation technique that does not rely on delivering stimuli at specific phases of the target rhythm is CR stimulation (Tass, 2003a). During CR stimulation, desynchronization is achieved by delivering phase-shifted stimuli to individual neuronal subpopulations. Other studies analyzed desynchronization by linear and nonlinear delayed feedback stimulation (Rosenblum M. and Pikovsky A., 2004; Rosenblum M. G. and Pikovsky A. S., 2004; Hauptmann et al., 2005a; Hauptmann et al., 2005b; Hauptmann et al., 2005c; Popovych et al., 2005; Popovych et al., 2006a; Popovych et al., 2006b; Pyragas et al., 2007; Popovych and Tass, 2010). Periodic stimulation at a suitable frequency may counteract neuronal synchrony by chaotic desynchronization (Wilson et al., 2011). Another study suggested closed-loop phasic burst stimulation during which burst stimuli are delivered at certain phases of the synchronous targeted rhythm (Holt et al., 2016). These phases were calculated based on online estimations of the phase response curve of the collective rhythm (Holt et al., 2016).

These stimulation techniques were originally developed for networks with fixed synaptic connections. The brain, however, is subject to synaptic plasticity and reorganizes network connectivity continuously (Liu et al., 2015; Van Ooyen and Butz-Ostendorf, 2017). One of these plasticity mechanisms is spike-timing-dependent plasticity (STDP), where the synaptic strengths are modulated according to the relative timings of post- and presynaptic spikes (Markram et al., 1997; Abbott and Nelson, 2000; Caporale and Dan, 2008). In many brain regions, STDP leads to a strengthening of synapses if the postsynaptic spike follows the presynaptic one and to a weakening of synapses if the presynaptic spike follows the postsynaptic one (Markram et al., 1997; Bi and Poo, 1998). STDP may enable the formation of neuronal assemblies (Litwin-Kumar and Doiron, 2014) and other network motifs (Ocker et al., 2015; Madadi Asl et al., 2018). It may also stabilize neuronal activity patterns, such as synchronous activity (Karbowski and Ermentrout, 2002), and lead to the coexistence of different stable states characterized by different activity patterns, such as desynchronized, synchronized, or cluster states (Seliger et al., 2002; Zanette and Mikhailov, 2004; Tass and Majtanik, 2006; Maistrenko et al., 2007; Masuda and Kori, 2007; Aoki and Aoyagi, 2009; Berner et al., 2020; Yanchuk et al., 2020).

Computational studies in plastic neuronal networks found that CR stimulation may induce not only acute desynchronization during stimulation but also long-lasting desynchronization that outlasts stimulation (Tass and Majtanik, 2006). The authors found that plastic synapses weakened during stimulation, which drove the network from a stable synchronous state, with strong synaptic connections, into the attractor of a stable desynchronized state, with weak synaptic connections. Based on these computational findings, CR stimulation was suggested as a possible therapy for inducing long-lasting therapeutic effects in movement disorders and epilepsies (Tass and Majtanik, 2006). In PD, corresponding long-lasting desynchronization and therapeutic effects were later confirmed by preclinical studies (Tass et al., 2009; Tass et al., 2012b; Wang et al., 2016; Wang et al., 2022) and clinical studies in PD patients delivering CR to the STN (Adamchic et al., 2014) or using vibrotactile fingertip CR stimulation (Tass, 2017; Syrkin-Nikolau et al., 2018; Pfeifer et al., 2021; Tass, 2022).

A deeper understanding of the parameter dependence of long-lasting desynchronization effects is critical for future clinical applications. In an earlier computational study, the long-lasting desynchronization effects of CR stimulation were studied in networks of Hodgkin-Huxley neurons with STDP (Manos et al., 2018). CR stimulation with two spatio-temporal stimulus patterns was delivered: CR with rapidly varying sequence (CR RVS) and CR with slowly varying sequence (CR SVS). For CR RVS the stimulus pattern is varied after each subpopulation received one stimulus, whereas the pattern remains fixed for longer time intervals during CR SVS. The authors reported that long-lasting desynchronization effects of strong CR RVS and CR SVS stimulation were sensitive to the ratio between the stimulation frequency and the frequency of the pathological synchronous rhythm. In contrast, for weak stimulation amplitudes, long-lasting desynchronization effects became more robust with respect to changes of the stimulation frequency, especially for CR RVS, whereas CR SVS did not induce long-lasting desynchronization for certain unfavorable frequency ratios. A similar observation was made in a later computational study analyzing long-lasting effects of CR RVS stimulation in networks of leaky integrate-and-fire (LIF) neurons with STDP (Kromer and Tass, 2020). There, it was hypothesized that certain unfavorable stimulation frequencies might lead to resonances and stabilize synchronous neuronal activity. Another computational study using networks of LIF neurons with STDP found that the stimulation frequency also needs to be adjusted to the time scales of synaptic long-term depression (LTD) and long-term potentiation (LTP) (Kromer et al., 2020). From a theoretical standpoint, these frequency dependences might limit clinical applicability as different PD symptoms have been associated with excessive neuronal synchrony in different frequency bands. Specifically, dyskinesia and tremor have been associated with synchronous neuronal activity in the theta band (3–10 Hz) (Brown, 2003; Steigerwald et al., 2008; Tass et al., 2010; Contarino et al., 2012). In contrast, rigidity and bradykinesia have been associated with synchronized beta-band activity (13–30 Hz) (Kühn et al., 2006; Weinberger et al., 2006). Furthermore, multiple central oscillators were found to underlie the generation of tremor (Raethjen et al., 2000).

To improve the robustness of long-lasting desynchronization effects with respect to changes of the stimulation frequency, recent computational studies suggested temporal randomization of the CR stimulus pattern. Typically CR stimulation is delivered with a fixed cycle period, such that each stimulation site is activated exactly once per cycle (Tass, 2003a). A recent series of computational and theoretical studies found that the frequency robustness of long-lasting desynchronization effects increased if stimulus delivery times were randomized. In particular, stimulus deliveries according to a Poisson pulse train (Kromer and Tass, 2020; Khaledi-Nasab et al., 2021a) and the addition of random jitters to the individual stimulus delivery times within individual CR cycles were considered (Khaledi-Nasab et al., 2021b).

In the present study, we suggest CR stimulation with randomized stimulus amplitudes to increase the robustness of long-lasting desynchronization effects of CR stimulation with respect to the stimulation frequency. We perform computer simulations of networks of LIF neurons with STDP and study long-lasting effects of regular CR (Popovych and Tass, 2012; Zeitler and Tass, 2015) and temporally uncorrelated CR (referred to as uncorrelated multichannel noisy stimulation in Zeitler and Tass (2018)). Based on results from previous studies, we expect long-lasting effects of regular CR to show a pronounced frequency dependence (Manos et al., 2018; Kromer et al., 2020) and long-lasting effects of temporally uncorrelated CR to be more robust with respect to changes of the stimulation frequency (Kromer et al., 2020; Kromer and Tass, 2020; Khaledi-Nasab et al., 2021a; Khaledi-Nasab et al., 2021b). We compare long-lasting desynchronization effects for both patterns with and without randomized stimulus amplitudes. Remarkably, double-random CR stimulation, i.e., CR stimulation with randomization of stimulus amplitudes as well as temporally uncorrelated stimulus times, substantially improves the long-term desynchronization outcome for high CR stimulation frequencies (compared to the intrinsic frequency of the abnormal target rhythm) and large numbers of separately stimulated subpopulations.

This present paper is organized as follows: In section 2, we present details on the plastic neuronal network model used throughout the paper. Then, the different stimulation patterns and stimulus randomizations are introduced. We also present the detailed derivation of our theoretical approximations for the synaptic weight dynamics during ongoing stimulation. In section 3, we compare results for CR stimulation patterns with constant stimulus amplitudes and with randomized stimulus amplitudes. In particular, we consider uniformly distributed stimulus amplitudes and binarily distributed stimulus amplitudes. For the latter, a random fraction of stimuli is removed from the pattern (by setting its amplitudes to zero) while the remaining stimuli possess a constant amplitude. Finally, in section 4, we discuss our results in the context of the current literature and point out possible consequences for clinical studies.

Simulations were performed for networks of N = 10³ excitatory LIF neurons with STDP. All parameters were chosen according to Kromer et al. (2020), Khaledi-Nasab et al. (2021a), and Khaledi-Nasab et al. (2021b) such that a stable desynchronized state, with weak synaptic connections, and a stable synchronized state, with strong synaptic connections, coexist. All equations and parameters are given in Kromer et al. (2020), Khaledi-Nasab et al. (2021a), and Khaledi-Nasab et al. (2021b) and are presented in Supplementary Material for the reader’s convenience.

We study the consequence of randomized stimulus amplitudes on the acute and long-lasting effects of CR and tCR stimulation. To this end, we perform simulations of networks of 10³ LIF neurons with STDP (see Methods) and compare the results to our theoretical predictions (see Methods). Simulated networks were prepared in the synchronized state (see Methods). If not stated otherwise, stimulation was delivered for 1,000 s.

Motivated by therapeutic brain stimulation in the context of Parkinson’s disease, we focus on the effect of stimulation on neuronal synchrony. We distinguish between acute effects and long-lasting effects on synchronization (see Methods). Acute effects were quantified by the acute Kuramoto order parameter, ρ ̄ ac , quantifying the degree of neuronal synchrony during stimulation. Long-lasting effects were quantified by the long-lasting Kuramoto order parameter, ρ ̄ ll , quantifying the degree of neuronal synchrony 1,000 s after cessation of stimulation. Additionally, we analyze the mean synaptic weight shortly before the cessation of stimulation ⟨w⟩_ac. ⟨w⟩_ac quantifies the effect of stimulation on the synaptic connectivity. In particular, stimulation led to an overall weakening of synaptic connections if ⟨w⟩_ac < w ₀ and to an overall strengthening of synaptic weights if ⟨w⟩_ac > w ₀. w ₀ ≈ 0.38 is the mean synaptic state in the stable synchronized state. In an earlier study, we found that ⟨w⟩_ac is highly predictive of the long-term outcome of stimulation. For ⟨w⟩_ac⪅0.27, the network typically approached the stable desynchronized state, whereas it approached the synchronized state for larger ⟨w⟩_ac (see Figure 3 in Kromer et al. (2020)). Furthermore, the immediate after-effect on synchronization is quantified by the Kuramoto order parameter ρ ̄ af , evaluated shortly after cessation of stimulation (see Methods for details). We are particularly interest in stimulation techniques that cause long-lasting desynchronization, i.e., desynchronization effects that persist after cessation of stimulation.

In the present paper, we studied the effect of isolated stimulus amplitude randomization, isolated stimulus timing randomization as well as combinations thereof on the parameter robustness of long-lasting desynchronization effects of CR stimulation in networks of LIF neurons with STDP. Long-lasting desynchronization is observed when a sufficient weakening of synaptic connections occurs during the stimulation. Such weakening may drive the network into a stable desynchronized state, which allows desynchronized activity to persist after cessation of stimulation (Tass and Majtanik, 2006; Kromer et al., 2020). We considered regular CR and temporally uncorrelated CR stimulation, i.e., CR with stimulus timings that are randomized within CR cycles and uncorrelated between stimulus channels. In addition, two distributions of stimulus amplitudes were considered: a uniform distribution and a binary distribution. For binarily distributed stimulus amplitudes, a random fraction of stimuli was removed from the stimulus pattern (A _stim = 0), while the remaining stimuli had amplitude A _stim = 1. We find that randomization of stimulus amplitudes extends the parameter region of stimulation-induced long-lasting desynchronization effects towards high stimulation frequencies. However, it slightly reduces desynchronization effects for low stimulation frequencies.

First, we studied regular CR and temporally uncorrelated CR with constant stimulus amplitudes. We set A _stim = 1 which corresponds to a strong stimulation regime where stimuli trigger neuronal spikes (Kromer et al., 2020). The regular CR pattern was used in previous computational studies (Popovych and Tass, 2012; Zeitler and Tass, 2015), preclinical studies (Tass et al., 2009; Tass et al., 2012b; Wang et al., 2016), and clinical studies in Parkinson’s disease patients where it was either delivered using DBS electrodes (Adamchic et al., 2014) or using vibrotactile fingertip stimulation (Syrkin-Nikolau et al., 2018; Pfeifer et al., 2021). It was also delivered to suppress binge alcohol drinking in mice (Ho et al., 2021) and as a treatment for tinnitus (Tass et al., 2012a; Munjal et al., 2021). Temporally uncorrelated CR was introduced in an earlier computational study, where it was referred to as uncorrelated multichannel noisy stimulation (Zeitler and Tass, 2018). In networks of LIF neurons, both stimulation protocols yield pronounced long-lasting desynchronization effects for stimulation frequencies in the range from one to five times the frequency of the targeted synchronous rhythm. This frequency range is independent of the number stimulated subpopulations for temporally uncorrelated CR (Figure 3B). Contrastingly, it may shrink substantially for unfavorable numbers of subpopulations for regular CR (Figure 3A). This shrinking was first reported in Kromer et al. (2020) (note that shorter stimulation durations were used in Kromer et al. (2020)). In that study, the pronounced dependence on the number of subpopulations was explained as a delay-induced effect. More specifically, transitions between parameter regions with long-lasting desynchronization and parameter regions with long-lasting synchronization occurred when integer multiples of the minimal interstimulus interval, 1/f _CR N _s, equaled the delay time, τ (white dashed curves in Figure 3A). Thus, stimulus-triggered spikes cannot arrive at their postsynaptic neurons before the next stimulus is delivered (Kromer et al., 2020). This led to a sudden transition from stimulation-induced LTP (the presynaptic spike arrives shortly before the postsynaptic neuron fires a spike triggered by a later stimulus) to stimulation-induced LTD (the presynaptic spike arrives shortly after the postsynaptic neuron fires a spike triggered by a later stimulus). Remarkably, this effect does not occur for temporally uncorrelated CR rendering its long-lasting desynchronization effects more robust with respect to parameter changes (Figure 3B). This is in line with results from previous studies suggesting that a reduction of temporal correlations between stimuli improves the frequency robustness of long-lasting desynchronization effects (Kromer and Tass, 2020; Khaledi-Nasab et al., 2021a; Khaledi-Nasab et al., 2021b). Specifically, we added temporal jitter and found that this smoothes out these sudden transitions and leads to overall stronger stimulation-induced LTD (see Figure 3 in Khaledi-Nasab et al. (2021b)). However, we find that temporally uncorrelated CR performs slightly worse than CR at low stimulation frequencies (lower than the frequency of the targeted synchronous rhythm). This is in line with the results of Zeitler and Tass (2018). Zeitler et al. considered a network of coupled Hodgkin-Huxley neurons with excitatory and inhibitory synaptic connections and adjusted the stimulation frequencies to that of the targeted synchronous rhythm. They reported on average weaker long-lasting effects for temporally uncorrelated CR than for the regular CR stimulation.

In the present paper, we considered random stimulus amplitudes. First, we distributed stimulus amplitudes uniformly between A _stim = 0 (no stimulus is delivered) and the maximum stimulus amplitude A _stim = 1. For A _stim = 1, the integral current over the excitatory part of a single stimulus is strong enough to drive the LIF neurons’ membrane potentials over the spiking threshold (Kromer et al., 2020). We find that this type of stimulus amplitude randomization substantially extends the frequency range for long-lasting desynchronization effects towards higher stimulation frequencies. For uCR stimulation, this mainly affects the parameter region of large N _s (Figure 4B). In contrast, if combined with temporal randomization, stimulation with uniformly distributed stimulus amplitudes induces long-lasting desynchronization effects for stimulation frequencies at least up to twelve times the frequency of the synchronous target rhythm and arbitrary values of N _s (see results for utCR stimulation in Figure 4). Long-lasting desynchronization can be achieved for all considered numbers of subpopulations. For regular CR, we find that stimulus amplitude randomization reduces the delay-induced effects mentioned above. This substantially extends the parameter region in which long-lasting desynchronization can be achieved. However, this occurs mainly in the range of large numbers of subpopulations (Figure 4B).

To put our results for larger numbers of subpopulations N _s in perspective, we note that the clinical proof of concept data obtained with CR stimulation delivered to the STN in Parkinson’s patients were obtained with only two or three active circular stimulation contacts being located in the STN (Adamchic et al., 2014; Ebert et al., 2014). In general, in clinical applications, the number of subpopulations N _s is limited by the electrode design, the number and dimension of stimulation contacts and the size of the target brain region (Krauss et al., 2021). While recent multisite stimulation electrodes possess large numbers of stimulation contacts (Steigerwald et al., 2019), e.g., up to 32 (Contarino et al., 2014), it remains to be shown whether such electrodes can be used to deliver stimuli to large numbers of separate neuronal subpopulations in target brain areas, such as the STN, independently. Larger numbers of smaller, e.g., segmented rather than circular, stimulation contacts might enable more focal stimulation by selecting groups of active stimulation contacts, hence, providing more favorable outcome. However, larger numbers of stimulation contacts do not imply a larger number of subpopulations N _s, since, due to the non-homogeneous anatomy of DBS targets, activation of particular subsets of stimulation contacts typically does not yield any therapeutic effects or even cause side effects. Further studies are required to understand the anatomy and physiology of the spatial dimensions of DBS target subpopulations to infer realistic ranges and, specifically, maximum values of N _s. With these caveats in mind, in Figures 4, 5 we used the range N _s⪅8 (pink dashed line) (Krauss et al., 2021) to illustrate a plausible clinically relevant range of N _s achievable with recent multi-contact DBS electrodes. In that range, however, uniform stimulus amplitude randomization alone does not improve the frequency robustness of long-lasting desynchronization effects compared to regular CR stimulation (Figures 4A,B,A”,B”) while binary stimulus amplitude randomization even worsens the frequency robustness of long-lasting desynchronization effects compared to regular CR (Figures 5A,B,A”,B”). By the same token, the advantage of stimulus timing randomization alone is rather limited (Figures 3A,B,A‴,B‴). However, double-random CR, i.e., CR with the combination of (uniform or binary) amplitude randomization as well as stimulus timing randomization, leads to more favorable long-term desynchronization for values of the stimulation frequency f _CR that exceed the intrinsic frequency of the abnormally synchronized target rhythm by a factor of five or more (Figures 4A,C,A”,C” as well as Figures 5A,C,A”,C”). If future clinical studies showed that a large number of neuronal subpopulations can be stimulated independently in related brain areas, double-random CR stimulation might induce more favorable long-lasting desynchronization effects, specifically for larger numbers of subpopulations N _s.

Since weak stimuli may not trigger neuronal spikes, stimulus amplitude randomization may lead to more sparse stimulation-induced spike patterns than stimulation with constant amplitudes. Intuitively, this may lead to more long time lags between arrivals of stimulus-triggered presynaptic spikes and stimulus-triggered postsynaptic spikes, which may favor LTD for STDP rules where LTD dominates over LTP for long time lags (τ ₊ < τ ₋ and β > 1 in Eq. 1). Previous studies analyzed the firing rate dependence of the synaptic weight dynamics in neuronal networks with STDP. They found that lower firing rates typically result in weaker synaptic weights for asymmetric STDP functions with strong LTP for short positive and LTD for negative time lags (Equation 1) (Sjöström et al., 2001; Burkitt et al., 2004). A sparsening of the stimulation-induced spike pattern might therefore reduce synaptic potentiation. This may explain the observed extension of the frequency range for long-lasting desynchronization towards higher stimulation frequencies. To test this hypothesis, we considered a second type of stimulus randomization where only a random fraction p of stimuli was delivered: binary stimulus randomization. We find that binary stimulus randomization leads to qualitatively similar effects as uniform stimulus randomization (Figures 4, 5). This supports the hypothesis that the observed improvement of long-lasting effects at high stimulation frequencies results from the effective sparsening of the stimulation-induced spike pattern.

For the same mean stimulus amplitude, we find that stimulation with uniform stimulus amplitude randomization performs better for high stimulation frequencies than stimulation with binary stimulus amplitude randomization. We suggest that this arises from differences in the statistics of time lags between stimulation-induced post- and presynaptic spikes. These time lags determine synaptic weight updates via the STDP function (Eq. 1). In the present paper, we considered an asymmetric STDP function, τ _R > 1, for which synaptic LTP dominates for time lags that are short compared to the STDP decay time τ ₊ (Eq. 1). Synaptic LTD dominates for long time lags. For uniform stimulus amplitude randomization, a large portion of stimuli only triggers spikes when the LIF neurons’ membrane potentials are close to the spiking threshold, i.e., after it left the refractory period following its previous spike. Such stimuli typically lead to rather large interspike intervals. In contrast, for binary randomization, spikes are triggered independently of the neurons’ membrane potential. Thus, for high stimulation frequencies, binary randomization results in more short interspike intervals. We therefore expect more short time lags between post- and presynaptic spikes during stimulation with binary stimulus amplitude randomization than during stimulation with uniform stimulus amplitude randomization, which would then translate into more pronounced stimulation-induced synaptic LTP and less pronounced LTD (compare Figures 4B-B”,C-C”, and Figures 5B-B”,C-C”).

We derived theoretical approximations for the mean rate of weight change during CR and temporally uncorrelated CR with binary distributed stimulus amplitudes. Our theoretical results extend earlier results from Kromer and Tass (2020), Kromer et al. (2020), and Khaledi-Nasab et al. (2021b) to the case of temporally uncorrelated CR stimulation, and, for the first time, incorporate stimulation with binary stimulus amplitude randomization. Our theoretical approach builds on the assumption that each stimulus triggers a neuronal spike, and each spike is triggered by a stimulus. Thus, neuronal spiking due to the intrinsic dynamics, noise, and synaptic input is neglected. Previous studies found that corresponding results well approximated the synaptic weight dynamics during ongoing strong (A _stim ≈ 1) and fast stimulation, i.e., when the stimulation frequency is higher than that of the synchronous target rhythm.

Our theory predicts that the synaptic weight dynamics during ongoing binary randomized stimulation depends on the statistics of stimulus trains delivered to the postsynaptic and the presynaptic neuron. We distinguish between intrapopulation synapses, where both neurons belong to the same subpopulation, and interpopulation synapses, where both neurons belong to different subpopulations. For intrapopulation synapses, the post- and presynaptic neuron receive stimuli simultaneously. We find that these synapses typically weaken. In contrast, for interpopulation synapses, the post- and presynaptic neuron typically receive stimuli at different times. Our theory predicts that interpopulation synapses strengthen in a large part of the (p, f _CR)-parameter space (Figure 6). We find that our simulation results are well-approximated by the zeros of the predicted mean rate of weight change of interpopulation synapses. In particular, long-lasting synchronization is observed when interpopulation synapses strengthen during stimulation. In contrast, long-lasting desynchronization occurs when interpopulation synapses weaken during stimulation (Figure 6). We find that deviations between theoretical predictions and simulations mostly occur for low values of p. Thus, when a significant portion of stimuli is removed from the stimulus pattern and the probability for long interstimulus intervals increases (Figure 6).

Our analysis of the CR patterns with binary stimulus amplitude randomization also yields information on the impact of stimulation parameters on long-lasting desynchronization effects. In particular, we studied the impact of the fraction of removed stimuli p and the stimulation frequency for fixed numbers of subpopulations. For temporally uncorrelated CR, we find that the higher the stimulation frequency, the more stimuli have to be removed to ensure long-lasting desynchronization (Figure 6C and D). In contrast, for regular CR stimulation with binary stimulus amplitude randomization, there is an intermediate frequency range where no long-lasting desynchronization is possible for certain unfavorable numbers of subpopulations (compare Figures 6A”,B”). This effect has a similar origin as the delay-induced effect described above. It occurs for stimulation frequencies where the presynaptic spikes arrive shortly before the next stimulus is delivered (Figure 5 and Kromer et al. (2020)).

Besides supporting long-lasting desynchronization for high stimulation frequencies, stimulus amplitude randomization also reduces the delivered stimulation current. This may be advantageous for possible clinical applications, as the integral modulus of the stimulation current is strongly related to the risk of unwanted side effects (Krack et al., 2002). We compared the stimulation-induced mean synaptic weight for the different stimulation patterns as a function of the integral stimulation current. We find that uniform randomization performs best for intermediate and high stimulation frequencies. The best results were obtained for double-random CR, i.e., temporally uncorrelated CR with uniform stimulus amplitude randomization (Figure 7). This indicates that uniform stimulus amplitude randomization does perform better for a fixed stimulation duration, and requires less integral stimulation current. This suggests that temporally uncorrelated CR with uniform stimulus amplitude randomization may be superior to other stimulation patterns if the stimulation frequency is substantially higher than that of the synchronous target rhythm. Using a basal ganglia model without synaptic plasticity, an earlier computational study analyzed the acute effects of CR stimulation with noise stimuli instead of charge-balanced pulses with random amplitudes (Liu et al., 2022). The authors found that such noisy stimulation is also efficient in inducing acute desynchronization. Here, we solely focused on long-lasting effects.

To the best of our knowledge, the present paper presents the first study that analyzes long-lasting desynchronization effects of stimulation patterns with randomized stimulus amplitudes. Some experimental studies considered DBS with temporally randomized stimulus patterns. However, these studies analyzed the acute effect of temporally randomized DBS and provided mixed results: Dorval et al. (2010) and Birdno et al. (2012) reported that temporally randomized DBS was inefficient in providing symptom alleviation (Dorval et al., 2010; Birdno et al., 2012), whereas Brocker et al. (2013) reported that irregular DBS led to improved performance of Parkinson’s disease patients in a fingertip task. It is currently not known, whether DBS with randomized stimulus amplitudes is safe and feasible.

Our promising results suggest that temporally randomized CR with uniform stimulus amplitude randomization might be suitable for desynchronizing brain rhythms across a wide frequency range. This might be advantageous in PD, where different symptoms are associated with abnormal synchrony in different frequency bands (Brown, 2003; Kühn et al., 2006; Weinberger et al., 2006; Steigerwald et al., 2008). Assuming that the interaction between rhythms is weak, the stimulation frequency could be adjusted to the fastest rhythm. Our computational results then suggest that such stimulation might also induce long-lasting desynchronization in subnetworks with pathological synchrony in lower frequency bands.

In the present study, we used a simple model of 10³ excitatory LIF neurons with STDP. Its low computational costs enabled us to perform detailed scans of the parameter space (Figures 3–6). In future studies, we anticipate working with more detailed neuron models specifically fit to target brain areas for HF DBS in PD, such as the basal ganglia (Terman et al., 2002; Fountas and Shanahan, 2017). Note that the more detailed conductance-based STN neuron model presented in (Terman et al., 2002) fired one spike per DBS pulse during HF DBS, similar to our LIF model (Rubin and Terman, 2004; Dorval et al., 2010). We anticipate studying the effect of stimulation in systems with multiple synchronous rhythms that interact. An interesting question for future studies would also be to which extend non-neural cells such as glia cells contribute to the therapeutic effects of stimulation. Growing evidence suggests that such cells play an important role in both PD pathogenesis (Booth et al., 2017) and the therapeutic mechanism of HF DBS (Vedam-Mai et al., 2012). We also aim at including physiological input from other brain areas. In a network model with additional, e.g., sensory input, after CR stimulation the synaptic connectivity pattern converged to a physiological connectivity pattern (Hauptmann and Tass, 2010). Given the complexity of such high-dimensional models with large numbers of parameters, predictions generated in simple models, e.g., the LIF model, may guide and inform the analysis of high-dimensional models, in a similar way as they have enabled the pre-clinical development (animal testing; (Tass et al., 2012b; Wang et al., 2016)) and clinical development (Adamchic et al., 2014) so far. Furthermore, we hope that our promising results inspire future studies using detailed computational models of target brain regions, e.g., for DBS in PD and other brain disorders, and preclinical and clinical studies on long-lasting therapeutic effects of randomized brain stimulation.