- Department of Neurosurgery, Stanford University, Stanford, CA, United States

Excessive neuronal synchrony is a hallmark of neurological disorders such as epilepsy and Parkinson's disease. An established treatment for medically refractory Parkinson's disease is high-frequency (HF) deep brain stimulation (DBS). However, symptoms return shortly after cessation of HF-DBS. Recently developed decoupling stimulation approaches, such as Random Reset (RR) stimulation, specifically target pathological connections to achieve long-lasting desynchronization. During RR stimulation, a temporally and spatially randomized stimulus pattern is administered. However, spatial randomization, as presented so far, may be difficult to realize in a DBS-like setup due to insufficient spatial resolution. Motivated by recently developed segmented DBS electrodes with multiple stimulation sites, we present a RR stimulation protocol that copes with the limited spatial resolution of currently available depth electrodes for DBS. Specifically, spatial randomization is realized by delivering stimuli simultaneously to *L* randomly selected stimulation sites out of a total of *M* stimulation sites, which will be called L/M-RR stimulation. We study decoupling by L/M-RR stimulation in networks of excitatory integrate-and-fire neurons with spike-timing dependent plasticity by means of theoretical and computational analysis. We find that L/M-RR stimulation yields parameter-robust decoupling and long-lasting desynchronization. Furthermore, our theory reveals that strong high-frequency stimulation is not suitable for inducing long-lasting desynchronization effects. As a consequence, low and high frequency L/M-RR stimulation affect synaptic weights in qualitatively different ways. Our simulations confirm these predictions and show that qualitative differences between low and high frequency L/M-RR stimulation are present across a wide range of stimulation parameters, rendering stimulation with intermediate frequencies most efficient. Remarkably, we find that L/M-RR stimulation does not rely on a high spatial resolution, characterized by the density of stimulation sites in a target area, corresponding to a large *M*. In fact, L/M-RR stimulation with low resolution performs even better at low stimulation amplitudes. Our results provide computational evidence that L/M-RR stimulation may present a way to exploit modern segmented lead electrodes for long-lasting therapeutic effects.

## 1. Introduction

Synchronization of coupled oscillators is observed in various fields of the natural sciences, for instance, in neurosciences (Steriade et al., 1990; Haken, 2006), medicine (Tass, 1999), physics (Pikovsky et al., 2001; Haken, 2012), biology (Winfree, 2001), and chemistry (Kuramoto, 2003). In the nervous system, synchronization of neuronal activity is critical for successful motor control (Andres and Gerloff, 1999) and information processing (Singer, 1990). However, excessive synchronization is associated with several neurological disorders, e.g., essential tremor, Parkinson's disease (PD) (Alberts et al., 1969; Nini et al., 1995), epilepsy (Mormann et al., 2000), and chronic subjective tinnitus (Eggermont and Tass, 2015).

High-frequency deep brain stimulation (HF DBS) is the standard of care for patients with advanced PD. HF DBS is delivered to target brain regions, such as the subthalamic nucleus (STN), through implanted lead electrodes. The mechanism of action of DBS is still a matter of debate (Ashkan et al., 2017; Jakobs et al., 2019; Lozano et al., 2019). As PD symptoms return shortly after cessation of stimulation, permanent delivery of HF DBS is required for persistent symptom suppression (Temperli et al., 2003). On the other hand, permanent stimulation increases the risk of side effects such as depression, cognitive decline, speech difficulty, instability, dyskinesia, and gait disorders (Rodriguez-Oroz et al., 2005; Temel et al., 2006). The risk of unwanted side effects may be reduced by a substantial reduction of the delivered stimulation current.

Early studies on desynchronization focused on single pulses delivered to a vulnerable phase of a collective oscillation (Mines, 1914; Winfree, 1977, 1980; Warman and Durand, 1989; Tass, 1999), followed by the development of composite single-channel (Tass, 2001, 2002; Zhai et al., 2005) and multi-channel stimuli (Tass, 2003) to further improve the robustness of the desynchronizing effects. In addition, linear or non-linear delayed feedback was used to desynchronize model networks (Rosenblum and Pikovsky, 2004a,b; Hauptmann et al., 2005a,b,c; Popovych et al., 2005, 2006a,b; Pyragas et al., 2007; Popovych and Tass, 2010). The latter approaches might clinically be applied by using the linear or non-linear delayed feedback as an envelope of pulse trains (Popovych et al., 2017a,b). By estimating phase response curves, researchers also identified well-tuned periodic stimulation (Wilson et al., 2011), and the delivery of well-timed stimulation bursts as possible desynchronizing stimulation approaches (Holt et al., 2016). A desynchronization approach that does not rely on well-timed stimuli is coordinated reset (CR) stimulation (Tass, 2003). During CR stimulation, desynchronization is achieved by delivering phase-shifted stimuli to multiple subpopulations of oscillators.

Originally, the mentioned desynchronization techniques have been developed for networks of oscillators with fixed connection strengths. In the brain, however, neuronal networks are subject to synaptic plasticity that alters synaptic weights according to neuronal activity. A prominent mechanism leading to adaptive connectivity is spike-timing dependent plasticity (STDP), which modifies the coupling strengths based on the relative timing of post- and presynaptic spikes (Markram et al., 1997; Abbott and Nelson, 2000; Caporale and Dan, 2008). In several brain areas, STDP strengthens synapses if the postsynaptic neuron fires shortly after the presynaptic one, otherwise the connections become weaker (Markram et al., 1997; Bi and Poo, 1998). Plasticity mechanisms can stabilize certain activity patterns, such as synchronized activity (Karbowski and Ermentrout, 2002) and may lead to the formation of strongly connected neuronal assemblies (Litwin-Kumar and Doiron, 2014). A recent study showed that in presence of STDP, self-organized clusters could emerge, whereby the clusters divide the networks into groups that are synchronized at different firing frequencies (Röhr et al., 2019). Furthermore, the interplay of network adaptation and collective spiking activity can lead to the coexistence of distinct stable states, such as synchronized, desynchronized, and 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; Röhr et al., 2019; Berner et al., 2020; Yanchuk et al., 2020).

Stimulation-induced changes of synaptic connectivity may drive the network into an attractor of a stable desynchronized state and cause long-lasting desynchronization (Tass and Majtanik, 2006). Such long-lasting desynchronization may follow after coordinated reset stimulation as shown by extensive theoretical (Tass and Majtanik, 2006; Kromer and Tass, 2020; Kromer et al., 2020) and computational studies (Tass and Majtanik, 2006; Hauptmann and Tass, 2009; Popovych and Tass, 2012; Lourens et al., 2015; Manos et al., 2018). Corresponding long-lasting desynchronization and therapeutic effects have been confirmed experimentally (Tass et al., 2009), as well as in preclinical (Adamchic et al., 2014; Wang et al., 2016) and clinical studies (Adamchic et al., 2014).

In preclinical and clinical studies, the frequency of CR stimulation has been adjusted to the frequency of the synchronous rhythm as measured by the dominant peak in the power spectrum of the local field potential (Tass et al., 2012; Adamchic et al., 2014; Wang et al., 2016). This parameter choice is motivated by the original idea of CR stimulation; to excite higher-order modes of the Kuramoto order parameter (Tass, 2003). Additionally, recent computational studies indicate that long-lasting desynchronization effects of CR stimulation are more pronounced for well-adjusted CR frequencies (Manos et al., 2018; Kromer and Tass, 2020; Kromer et al., 2020), this includes adjusting the stimulation frequency to the dominant neuronal rhythm (Tass, 2003; Adamchic et al., 2014; Manos et al., 2018) or the STDP time scale (Kromer et al., 2020). This may limit the clinical applicability of CR stimulation as a treatment for PD as different symptoms are associated with pathological synchrony in different frequency bands. In more detail, excessive synchronization of basal ganglia activity in the theta band (3−10 Hz) is associated with symptoms such as dyskinesia and tremor (Brown, 2003; Steigerwald et al., 2008; Tass et al., 2010; Contarino et al., 2012), while synchronized activity in the beta band (13−30 Hz) is associated with rigidity and bradykinesia (Kühn et al., 2006; Weinberger et al., 2006).

In order to increase the robustness of long-lasting effects with respect to stimulation parameters, such as the stimulation frequency, Kromer and Tass (2020) suggested a Random Reset (RR) stimulation protocol. RR refers to the delivery of stimuli in a temporally and spatially randomized manner. In their theoretical work, temporal randomization was realized by delivering stimuli at random times, with exponentially distributed interstimulus intervals. Thus, stimulation times followed a Poisson spike train. Spatial randomization was realized by randomly selecting 50% of the neurons for stimulus application at each stimulation time—irrespective of the neurons' locations relative to realistic spatial stimulation profiles. This implicitly assumed “microscopic control,” i.e., that even nearby neurons can be stimulated independently, which is not possible in DBS-like setups. Remarkably, the suggested RR stimulation method led to robust, long-lasting desynchronization effects after stimulation ceases even though the neurons remained partially synchronized during the entire stimulation period. This was achieved by a pronounced stimulation-induced decoupling of the neurons. Therefore, decoupling stimulation was suggested as the primary goal in order to weaken synaptic connections rather than counteracting synchronization as in previous approaches (Kromer and Tass, 2020).

Segmented depth electrodes for DBS enable spatially selective steering of stimulation current (Krack et al., 2003; Buhlmann et al., 2011; Steigerwald et al., 2019; Krauss et al., 2020). However, so far, a pressing question regarding a possible implementation of RR stimulation using available DBS electrodes is that to which extent the observed decoupling effects rely on the spatial randomization. In particular, whether it is really necessary to deliver randomly timed stimuli to individual neurons or whether delivery to macroscopic neuronal subpopulations is sufficient. Experimentally using segmented DBS electrodes that allow for independent activation of multiple stimulation contacts, one can deliver stimuli to neuronal subpopulations. Traditionally DBS is delivered through a flexible cylinder with 4 stimulation contacts (Gielen, 2001; Butson and McIntyre, 2005). In order to deliver stimuli to individual neuronal subpopulations, modern electrodes are capable of directional current steering (Buhlmann et al., 2011; Steigerwald et al., 2019; Krauss et al., 2020)

Directional steering allows for realizing spatiotemporal current profiles by superposition of stimuli delivered to individual stimulation contacts (Buhlmann et al., 2011; Steigerwald et al., 2019). To improve spatial selectivity, recent research in electrode design is devoted to segmented multisite electrodes with increasing numbers of stimulation contacts (Buhlmann et al., 2011; Steigerwald et al., 2019; Krauss et al., 2020). For instance, certain designs allow for selective activation of up to 32 stimulation contacts (Contarino et al., 2014; Steigerwald et al., 2019).

In the present paper, we study a new implementation of RR stimulation using available DBS electrodes, where we introduce a version of RR stimulation in which individual stimuli are simultaneously delivered to *L* out of *M* randomly selected spatially coherent neuronal subpopulations, here called L/M-RR stimulation. This approach only requires “mesoscopic control,” i.e., independent stimulation of neuronal subpopulations, and not “microscopic control,” i.e., independent stimulation of single neurons, as in the approach of Kromer and Tass (2020). For our analysis, we use a combination of theoretical predictions and simulations of networks of leaky integrate-and-fire (LIF) neurons with STDP. While our theory predicts efficient decoupling for a wide range of stimulation parameters, it reveals qualitative differences between low and high-frequency L/M-RR stimulation, section 3.1. While low-frequency L/M-RR stimulation yields parameter robust decoupling and related long-lasting effects, the performance of high-frequency L/M-RR stimulation is limited to small numbers of simultaneously stimulated subpopulations. These qualitative differences are present for a wide range of stimulation amplitudes and render strong high-frequency L/M-RR stimulation ineffective in terms of long-lasting after-effects, see section 3.2. Finally, in section 3.3, we analyze how the size of individual subpopulations influence the long-lasting effects. Remarkably, we find that stimulation of large subpopulations yields better results for weak stimulation, see section 3.3. Simultaneous stimulation of large neuronal subpopulations may hence be advantageous for possible realizations of L/M-RR for DBS.

## 2. Models and Methods

### 2.1. Neuronal Network Model

We consider a network of *N* conductance-based LIF neurons with STDP previously presented in Kromer et al. (2020). Throughout the paper, we fix the network size to *N* = 1, 000. Neurons are organized along the *x*-axis. Individual neurons' center locations *x*_{i} are uniformly distributed in the interval *x*_{i} ∈ [−2.5, 2.5] mm, which is motivated by the width used in an ellipsoidal volume approximation of the STN in detailed computational studies of STN DBS (Ebert et al., 2014). Random excitatory synaptic connections are added such that the total connectivity is 7%. The probability for two neurons to form a synaptic connection depends on the distance between the neurons as *p*∝exp((|*x*_{j} − *x*_{i}|)/0.5mm) (Ebert et al., 2014).

The subthreshold dynamics of the membrane potential *v*_{i} of neuron *i* obeys

*C*_{i} is the membrane capacitance, *v*_{rest} the resting potential, *g*_{leak} = 0.02 mS/cm^{2} the leakage conductance, ${I}_{i}^{\mathrm{\text{syn}}}(t)$ the synaptic input current, *S*_{i}(*t*) the stimulation current, and ${I}_{i}^{\mathrm{\text{noise}}}(t)$ the noisy input current.

Spiking occurs when *v*_{i} crosses a dynamic threshold potential ${v}_{i}^{\mathrm{\text{th}}}$ given by

Here, τ_{th} is the threshold time constant and ${v}_{\mathrm{\text{rest}}}^{\mathrm{\text{th}}}$ the stationary threshold potential. Artificial spikes are realized by setting the membrane potential *v*_{i} → *v*_{spike} for a time period of *t*_{spike} after a threshold crossing. Afterwards, a reset is performed by setting *v*_{i} → *v*_{reset} and ${v}_{i}^{\mathrm{\text{th}}}\to {v}_{\mathrm{\text{spike}}}^{\mathrm{\text{th}}}$.

Throughout the paper, we use the same parameters as in Kromer et al. (2020): *v*_{rest} = −38 mV, ${v}_{\mathrm{\text{rest}}}^{\mathrm{\text{th}}}=-40$ mV, *t*_{spike} = 1 ms, τ_{th} = 5 ms, *v*_{spike} = 20 mV, *V*_{reset} = −67 mV. The *C*_{i}'s follow a normal distribution with a mean value of 〈*C*_{i}〉 = 3 μF/cm^{2} and a standard deviation of 0.05〈*C*_{i}〉.

Excitatory synaptic input ${I}_{i}^{\mathrm{\text{syn}}}(t)$ to neuron *i*, is given by

where ${g}_{i}^{\mathrm{\text{syn}}}$ is the synaptic conductance, κ = 8 mS/cm^{2} is the coupling strength, *v*_{syn} = 0 mV is the synaptic reversal potential, τ_{syn} = 1 ms the synaptic time scale, *w*_{ji} ∈ [0, 1] is the synaptic weight between presynaptic neuron *j* and postsynaptic neuron *i*. The first sum runs over all presynaptic neurons, and the second sum runs over the spikes of the presynaptic neuron *j*. *G*_{i} is the set of indices of all presynaptic neurons to neuron *i*. ${t}_{{l}_{j}}^{j}$ is the *l*_{j}th spike time of neuron *j*. We consider homogeneous synaptic delays of *t*_{d} = 3 ms.

In addition to presynaptic input from other neurons in the network, each neuron *i* receives noisy input, e.g., from other brain regions. The resulting input current ${I}_{i}^{\mathrm{\text{noise}}}$ is obtained by feeding independent presynaptic Poisson spike trains with firing rate *f*_{noise} = 20 Hz into excitatory synapses on each neuron *i* (Ebert et al., 2014)

Where ${g}_{i}^{\mathrm{\text{noise}}}(t)$ is the synaptic conductance, τ_{syn} = 1 ms, and *v*_{syn} = 0 mV. The noise intensity is controlled by the parameter *D* = 0.026 *mS*/*cm*^{2} scaling the strength of the Poisson input.

### 2.2. Spike-Timing Dependent Plasticity (STDP)

During ongoing spiking, synaptic weights *w*_{ij} evolve according to a nearest-neighbor STDP scheme (Burkitt et al., 2004). Following previous studies on CR stimulation (Popovych and Tass, 2012), we consider a scheme where each arrival of a presynaptic spike at a postsynaptic neuron *j* (at time *t* = *t*_{i} + *t*_{d}) and each postsynaptic spike (at time *t* = *t*_{j}) cause an update of the synaptic weight, i.e., *w*_{ij} → *w*_{ij} + *W*(*t*_{j} − (*t*_{i} + *t*_{d})). Here, *t*_{i} denotes the spike time of the presynaptic spike. *W*(*t*) is the STDP function and is given by two exponentials (Song et al., 2000; Kromer and Tass, 2020)

Δ*t* = *t*_{j} − (*t*_{i} + *t*_{d}) is the time lag between postsynaptic spike times and presynaptic spike arrival times. η = 0.02 scales the weight update per spike, τ_{+} = 10 ms and τ_{−} = τ_{+}τ_{R} are the STDP decay times for long-term potentiation (LTP) and long-term depression (LTD), respectively, τ_{R} = 4 yields asymmetry in STDP decay times, and β = 1.4 scales the ratio of overall LTD to LTP. These STDP parameters lead to bistability between a strongly connected state with synchronized neuronal activity and a weakly connected state with asynchronous neuronal activity (Kromer and Tass, 2020; Kromer et al., 2020).

### 2.3. Quantification of Synchronization

In order to quantify the degree of in-phase synchronization, we calculate the time-averaged Kuramoto order parameter (Kuramoto, 1984)

Here, ψ_{k}(*t*) is a phase function that increases linearly in time during individual interspike intervals of neuron *k*, i.e., ${\psi}_{k}(t)=2\pi ((t-{t}_{l}^{k})/({t}_{l+1}^{k}-{t}_{l}^{k})+l)$ for $t\in [{t}_{l}^{k},{t}_{l+1}^{k})$ (Rosenblum et al., 2001). ρ_{Δ} quantifies the degree of synchronization of a population of *N* neurons during the time interval Δ = 10 s. Perfect in-phase synchronization results in ρ_{Δ} = 1 whereas ρ_{Δ} ≈ 0 refers to absence of in-phase synchronized neuronal activity.

### 2.4. L/M-Random Reset Stimulation

Throughout the paper, we deliver RR stimulation to a randomly selected group of neuronal subpopulations. RR stimulation is characterized by the delivery of temporally and spatially randomized stimulus patterns. Temporal randomization is realized by delivering stimuli at random times *s*_{k}. Interstimulus intervals *S*_{k} = *s*_{k+1} − *s*_{k} are distributed according to an exponential distribution with minimum interstimulus interval τ_{Λ} (Kromer and Tass, 2020)

Where τ_{Λ} = 1/130 s is the minimum interstimulus interval which corresponds to a maximal stimulation frequency of 130 Hz. This frequency is often used in clinical DBS studies (Krauss et al., 2020).

Θ(*t*) is the Heaviside step function, τ_{RR} determines the average stimulation frequency *f*_{RR} by *f*_{RR}: = 1/(τ_{Λ} + τ_{RR}). Figure 1A shows the distribution of interstimulus intervals for *f*_{RR} = 30 Hz.

**Figure 1**. RR stimulation: **(A)** the probability distribution of interstimulus intervals for *f*_{RR} = 30 Hz. The dashed line shows the minimum interstimulus interval, ${\tau}_{\Lambda}={(130\text{}\mathrm{\text{Hz}})}^{-1}$. **(B)** The schematic of RR stimulation for *M* = 4 total sites where 2 sites are activated simultaneously at random times. The excitatory pulse has a duration of ν_{e} = 0.5 ms and an amplitude of ${{A}}_{\mathrm{\text{e}}}={A}_{\mathrm{\text{stim}}}\mu /{\nu}_{\mathrm{\text{e}}}$; separated by a gap of 0.2 ms, the inhibitory pulse has a duration of ν_{i} = 1.5 ms and an amplitude of ${{A}}_{\mathrm{\text{i}}}={A}_{\mathrm{\text{stim}}}\mu /{\nu}_{\mathrm{\text{i}}}$.

In the present study, spatial randomization is realized by delivering each stimulus (at time *s*_{k}) to a randomly selected group of *L* out of *M* neuronal subpopulations. Neurons are assigned to subpopulations according to their centers' locations. In particular, neuron *i* is considered to be part of subpopulation *l* if ${x}_{i}\in [-2.5+\frac{5(l-1)}{M},-2.5+\frac{5l}{M})$ mm, *l* = 1, 2, ..., *M*. Hence, *M* scales the required spatial resolution. This setup is motivated by the shape of commonly used cylindrical DBS electrodes with equidistantly placed stimulation contacts (Gielen, 2001; Krauss et al., 2020). Figure 2 illustrates the division into subpopulations for *M* = 4 sites. At each stimulation time *s*_{k}, a stimulus is delivered to *L* out of *M* randomly selected subpopulations without replacement. Thus, at each stimulation time the first electrode is selected uniformly at random with a probability of $\frac{1}{M}$, and for the second one we select one of the remaining electrodes with the probability of $\frac{1}{M-1}$ and so on. If not stated otherwise, we use *M* = 32 throughout the paper. This mimics a recently developed DBS lead electrode with up to 32 stimulation contacts that can be activated independently (Steigerwald et al., 2016).

**Figure 2**. Schematic of intrapopulation and interpopulation synapses for a total of *M* = *4* subpopulations. The blue arrows represent the intrapopulation synapses where both pre- and postsynaptic neurons are within the same group. The red arrows show the interpopulation synapses where the pre- and postsynaptic neurons are in different subpopulations. Note that neurons are distributed along the *x*-axis and the vertical axis is merely symbolic.

Individual stimuli are charge-balanced and consist of an excitatory and an inhibitory rectangular pulse. The excitatory pulse has a duration of ν_{e} = 0.5 ms and an amplitude of ${{A}}_{\mathrm{\text{e}}}={A}_{\mathrm{\text{stim}}}\mu /{\nu}_{\mathrm{\text{e}}}$; separated by a gap of 0.2 ms, the inhibitory pulse has a duration of ν_{i} = 1.5 ms and an amplitude of ${{A}}_{\mathrm{\text{i}}}=-{A}_{\mathrm{\text{stim}}}\mu /{\nu}_{\mathrm{\text{i}}}$ where μ = (*V*_{th,spike} − *V*_{reset})/〈*C*_{i}〉. Thus, a stimulus of stimulation strength *A*_{stim} = 1 will elevate the membrane potential above the spiking threshold no matter when the last spike of the neuron occurred. Figure 1B shows a schematic of our L/M-RR stimulation protocol for the case where each stimulus is delivered to *L* = 2 (out of *M* = 4) stimulation sites.

### 2.5. Calculation of Mean Rate of Weight Change

The development of RR stimulation was originally based on theoretical predictions of the stimulation-induced synaptic weight dynamics during randomized stimulation. The corresponding theoretical study can be found in Kromer and Tass (2020). In the following, we briefly state the main steps and expand their results to our L/M-RR stimulation protocol.

The approach of Kromer and Tass (2020) and others (Kempter et al., 1999; Burkitt et al., 2004) is based on the mean rate of weight change ${{J}}_{ij}$ for a single synaptic weight *w*_{ij} with presynaptic neuron *i* and postsynaptic neuron *j*. Averaging over a long time interval and a large number of realizations of the spike train, and assuming that post- and presynaptic neurons have the same mean firing rate *r*, the mean rate of weight change can be written as

For the case of different firing rates of neurons *i* and *j*, we refer to Kromer and Tass (2020). Here *G*_{ij}(*t*) is the average number of time lags *t* per spike in the infinitesimal interval *t* ∈ [*t, t* + *dt*) that contribute to weight updates for the chosen STDP scheme. Note that *G*_{ij}(*t*) is in general not normalized to one, but to the mean number of time lags per spike that contribute to weight updates.

Kromer and Tass (2020) derived results for *G*_{ij}(*t*) in the limit of strong and fast stimulation, where each stimulus triggers a spike, and all spikes are triggered by stimuli. Here *strong* refers to large stimulation amplitudes *A*_{stim} ≈ 1 and *fast* to stimulation frequencies that are large compared to the firing rate in the synchronous state.

In the limit of strong and fast stimulation, the neurons' firing rates equal the frequency of stimulus administration *r*_{i} = *r*_{j} = *f*_{RR}*L*/*M*, and the distribution of time lags that lead to weight updates due to STDP results from the distribution of interstimulus intervals between stimuli delivered to the post- and presynaptic neurons, respectively (Kromer and Tass, 2020). We distinguish between two contributions to these time lags *t* by setting *t* = *S* + ξ. Here, *S* is the difference between post- and presynaptic stimulation times and ξ characterizes the difference in delayed spiking responses to stimuli. Note that only time lags that lead to weight updates are considered. Therefore, *S* does not denote any interstimulus interval between stimuli delivered to the post- and presynaptic neurons, but only those for which the time lag between the triggered post- and presynaptic spike results in a weight update. Then, *G*_{ij}(*t*) results from the distributions of ξ and *S*. First, we consider the distribution of ξ. In the limit of strong and fast stimulation, each stimulus triggers a spike of a stimulated neuron and each spike is triggered by a stimulus. We assume that the time lag *t*′ between stimulus delivery and spiking response of the stimulated neuron is distributed according to a distribution λ(*t*′). Further, assuming that spiking responses of the pre- and postsynaptic neuron follow this distribution, we find that ξ is distributed according to $Z(\xi )={\int}_{-\infty}^{\infty}d{t}^{\prime}\text{}\lambda ({t}^{\prime})\lambda ({t}^{\prime}+\xi )$. Second, we denote the distribution of *S* as *p*_{ij}(*S*, ξ). Note that, in general, *p*_{ij}(*S*, ξ) depends on the realization of ξ. This is because certain realizations of ξ may change the order of spike and spike arrival times of the pre- and postsynaptic neurons' spikes. The latter may affect which time lags contribute to weight updates, and therefore, which interstimulus intervals *S* need to be considered.

In order to calculate *p*_{ij}(*S*, ξ), we follow the approach of Kromer and Tass (2020) and consider possible pairings of a presynaptic (postsynaptic) spike that is triggered by the *n*th stimulus with postsynaptic (presynaptic) spikes. We distinguish between two scenarios. In the first scenario, both post- and presynaptic neuron are stimulated simultaneously, while in the second one either the postsynaptic or the presynaptic neuron receives the *n*th stimulus.

In Kromer and Tass (2020) results for *p*_{ij}(*S*, ξ) for either case were derived. In the first case, the distribution of interstimulus intervals is given by

with

Here the *k*th summand contains the (*k* − 1)th convolution with the zeroth convolution referring to *P*(*S*) itself. The latter is given in Equation (7). δ(*x*) denotes the Dirac delta distribution.

As can be seen in Equation (9), the realization of ξ determines which interstimulus intervals *S* are considered for weight updates resulting from the current presynaptic spiking event. If ξ is larger than the synaptic delay *t*_{d} the presynaptic spike arrives before the postsynaptic spike and the arrival time is paired with the current postsynaptic spike (positive update) and the latest postsynaptic spike, triggered by an earlier stimulus delivered to the postsynaptic neuron (negative update). *F*(−*S*) denotes the probability that the latest stimulus was delivered to the postsynaptic neuron with an interstimulus interval *S* relative to the current stimulus. In contrast, if ξ is smaller than the delay time, the presynaptic spike arrives after the postsynaptic one and its arrival time is paired with the current postsynaptic spike (negative update) and the next postsynaptic spike that results from the next stimulus delivered to the postsynaptic neuron (positive update).

In the second case, *p*_{ij}(*S*, ξ) does not depend on ξ, if synaptic delays are short compared to interstimulus intervals and λ(*t*) is narrow compared to interstimulus intervals. *p*_{ij}(*S*, ξ) given by

Here, the current presynaptic spike arrival time is paired with postsynaptic spikes triggered by the latest (negative update) and the next stimulus delivered to the postsynaptic neuron (positive update), respectively.

Using Equations (9), (11), and (10), the distribution of time lags *G*_{ij}(*t*) results from

We separate synapses into two groups based on the probability of receiving stimuli simultaneously during L/M-RR stimulation. The first group consists of synapses that connect neurons within the same subpopulation. This group will be referred to as *intrapopulation* synapses. The second group consists of synapses between neurons belonging to different subpopulations. These synapses will be referred to as *interpopulation* synapses. The classification of intra- and interpopulation synapses is illustrated for a toy network in Figure 2.

First, we consider intrapopulation synapses. Corresponding quantities will be marked by the suffix “intra” in the following. Neurons connected by intrapopulation synapses always receive stimuli simultaneously. In consequence, the distribution of time lags is given by

Here the average is taken over all intrapopulation synapses. Using *G*^{intra}(*t*) in Equation (8) yields the expected rate of weight change for intrapopulation synapses

Accordingly, we mark quantities related to interpopulation synapses by the suffix “inter”. Neurons that are connected by interpopulation synapses belong to different subpopulations and only receive stimuli simultaneously when both subpopulations are selected for stimulus delivery. Given that one of the subpopulations is already selected, the probability to select the other one as well is (*L* − 1)/(*M* − 1), which yields

The average is taken over all interpopulation synapses. We introduce the expected rate of weight change for interpopulation synapses as

### 2.6. Simulation Details

Numerical integration of the LIF model presented in section 2.1 is performed using an explicit Euler integration scheme with an integration time step of 0.1 ms.

Equation (10) is evaluated numerically using a discretization of the time axis with bin size of 5*10^{−3} ms (Figures 4A,B,E,F,G). Plotted histograms in Figures 4C,D were obtained using a bin size of 0.5 ms. The sum is truncated after *k*_{max} summands, where *k*_{max} is the first integer larger than 500.0ms/τ_{Λ}. Convolutions are calculated using the *python* method *numpy.convolve* of *numpy* version 1.16.2.

## 3. Results

To study long-lasting desynchronization by L/M-RR stimulation, we perform numerical simulations of networks of LIF neurons with STDP. For the chosen parameters, a stable synchronized state with strong synaptic connections coexists with a stable desynchronized state with weak connections; see (Kromer et al., 2020). Networks with high initial mean weight approach the synchronized state, while networks with low initial mean weight approach the desynchronized state.

To prepare the network in the synchronized state, we choose a high initial mean weight, 〈*w*〉(*t* = 0) = 0.5. This was realized by selecting 50% of the synapses at random and setting their weights to one, while the other synaptic weights were set to zero. Then, the network is simulated for 500 s in order to reach the synchronous state, see (Kromer et al., 2020) for details.

After preparation, L/M-RR stimulation is delivered for 500 s. Afterward, we continue the simulation for another 1, 000 s to explore potential long-lasting effects of L/M-RR stimulation. To quantify the effect of L/M-RR stimulation, we analyze its acute and long-lasting effects by evaluating the mean synaptic weight 〈*w*〉 and the time-averaged Kuramoto order parameter ρ_{Δ}, Equation (6). To quantify acute effects, we evaluate 〈*w*〉 and ρ_{Δ} when the stimulation ceases (*t* = 995 s and Δ = 10 s in Equation 6). Furthermore, to quantify long-lasting effects, 〈*w*〉 and ρ_{Δ} are evaluated 1, 000 s after cessation of stimulation (*t* = 1, 995 s and Δ = 10 s in Equation 6).

Figure 3 shows representative time traces of the mean synaptic weight and the time-averaged Kuramoto order parameter obtained from simulations for different stimulation periods *T*. While the Kuramoto order parameter decreases within seconds after stimulation onset (Figure 3B) the dynamics of the mean synaptic weight is significantly slower, as shown in Figure 3A. As a consequence, the system might not reach the attractor of the stable desynchronized state for short stimulation periods due to insufficient decoupling. Throughout the paper, we fix the stimulation time to 500 s, which turned out to be sufficient for most parts of the parameter space.

**Figure 3**. Time traces for different stimulation durations. **(A)** The mean weights, 〈*w*〉, and **(B)** the order parameter, ρ_{Δ}, for four different stimulation durations *T*. The dashed line in **(A)** characterizes the weights in the stable synchronized state. Here, *f*_{RR} = 10 Hz, *A*_{stim} = 1, *L* = 5, and *M* = 32. The Kuramoto order parameter in **(B)** is calculated every 10 s by averaging over non-overlapping time windows Δ = 10 s.

Throughout the paper, we time-average the order parameter over Δ = 10 s. Furthermore, results are averaged over three network realizations. We find that results agree qualitatively among these random network realizations. Therefore, and due to the high computational costs, averaging over a large ensemble of network realizations was not performed.

### 3.1. Robust and Long-Lasting Desynchronization for Strong Stimulation

First, we consider the case of strong stimulation *A*_{stim} = 1 in which neuronal spikes follow the stimulus pattern. Theoretical predictions for the mean rates of weight changes of intrapopulation synapses, ${{J}}^{\mathrm{\text{intra}}}$, and interpopulation synapses, ${{J}}^{\mathrm{\text{inter}}}$, are shown in Figures 4A,B for a wide range of stimulation frequencies and fractions of simultaneously stimulated subpopulations *L*/*M*. Negative rates of weight change, ${{J}}^{\mathrm{\text{intra}}}<0$, and ${{J}}^{\mathrm{\text{inter}}}<0$, indicate a weakening of corresponding synapses during stimulation.

**Figure 4**. Theoretical predictions for synaptic weight dynamics. **(A,B)** Expected mean rate of weight change for intrapopulation synapses ${{J}}^{\mathrm{\text{intra}}}$, Equation (14) **(A)**, and interpopulation synapses ${{J}}^{\mathrm{\text{inter}}}$ **(B)**, Equation (16). Colored curves separate regions with expected strengthening from those with expected weakening of synaptic weights. **(C,D)** Predicted distributions of time lags *G*^{intra}(*t*) **(C)**, Equation (13), and *G*^{inter}(*t*) **(D)**, Equation (15) (dotted black curves), compared to simulation results (colored histograms). **(E–G)** Theoretical predictions (lines) of mean synaptic weights for intra- (blue) and interpopulation (orange) synapses compared to simulation results (circles) for different stimulation frequencies *f*_{RR} and fractions of simultaneously stimulated subpopulations *L*/*M*. Parameter combinations correspond to the three possible qualitative effects of L/M-RR stimulation: (i) weakening of all synapses **(E)**, (ii) weakening of intrapopulation synapses **(F)** and (iii) strengthening of all synapses **(G). (H–K)** Overall mean synaptic weight **(H,J)** and time-averaged Kuramoto order parameter **(I,K)** as function of *f*_{RR} and *L*/*M*. Acute values **(H,I)** are compared to long-lasting ones **(J,K)**. Colored curves from **(A,B)** are shown to separate parameter regions (i)–(iii) with predicted strengthening and weakening of intra (blue) and interpopulation (orange) synapses, respectively. Parameters: *M* = 32 with *L* = 15 and *f*_{RR} = 60 Hz **(C–E)**; *L* = 15 and *f*_{RR} = 100 Hz **(F)**; and *L* = 25 and *f*_{RR} = 100 Hz **(G)**. We used λ(*t*) = δ(*t*) for theoretical predictions. “ac” and “ll” refer to acute and long-lasting effects, respectively.

A detailed comparison between simulation results and theory is presented in Figures 4C–G. We find an excellent quantitative agreement for the distributions of time lags, Figures 4C,D, and trajectories of the mean synaptic weights of both intra- and interpopulation synapses after the onset of stimulation, Figures 4E–G.

In Figures 4H,I, we show results for the mean synaptic weight, Figure 4H, and the Kuramoto order parameter, Equation (6), Figure 4I, shortly before stimulation ceases. Our theory separates the parameter space into three regions: (i) a region with ${{J}}^{\mathrm{\text{intra}}}<0$ and ${{J}}^{\mathrm{\text{inter}}}<0$ where L/M-RR stimulation decouples all neurons, (ii) a region where only intrapopulation synapses are weakened (${{J}}^{\mathrm{\text{intra}}}<0$ and ${{J}}^{\mathrm{\text{inter}}}>0$), and (iii) a region where all synapses are strengthened (${{J}}^{\mathrm{\text{intra}}}>0$ and ${{J}}^{\mathrm{\text{inter}}}>0$), see Figure 4H. Our simulation results show that in region (i) L/M-RR stimulation causes long-lasting desynchronization, while the system returns to the synchronized state in regions (ii) and (iii), see Figures 4J,K.

Next, we analyze the impact of the fraction of simultaneously active subpopulations *L*/*M* on long-lasting effects. Two qualitatively different frequency ranges can be found in Figures 4A,B. For low stimulation frequencies, *f*_{RR} < 60 Hz, stimulation improves as *L*/*M* increases. In contrast, most pronounced long-lasting effects for high stimulation frequencies, *f*_{RR} > 60 Hz, are observed at finite values of *L*/*M*.

### 3.2. Moderate Stimulation Yields Most-Pronounced Long-Lasting Effects

We study the impact of the stimulation amplitude *A*_{stim} on acute and long-lasting effects of L/M-RR stimulation for low and high stimulation frequencies. Results are shown in Figure 5. We find that the long-lasting effects of low-frequency stimulation differ significantly from those of high-frequency stimulation. For low stimulation frequencies, stronger stimulation yields better results, see Figures 5A–C. Furthermore, for a smaller fraction of *L*/*M*, stronger stimulation is required to get sustained long-lasting effects. In contrast, for high stimulation frequencies, moderate stimulation yields pronounced long-lasting desynchronization, see Figure 5F, while weak and strong stimulation does not induce long-lasting desynchronization for a wide range of fractions *L*/*M*. Remarkably, we also find that the system returns to the synchronized state after L/M-RR stimulation with very low and very high ratios *L*/*M*, Figure 5F. Therefore, L/M-RR stimulation is most efficient for moderate stimulation amplitudes and intermediate ratios of simultaneously stimulated neuronal subpopulations.

**Figure 5**. Acute and long-lasting effects of L/M-RR stimulation. Effects of low-frequency **(A,B)** and high-frequency stimulation **(D,E)** shortly before the cessation of stimulation. Mean synaptic weight **(A,D)** and time-averaged Kuramoto order parameter **(B,E)** are shown. Additionally, we show results for the time-averaged Kuramoto order parameter evaluated 1, 000 s after cessation of stimulation in **(C)** (low-frequency stimulation) and **F** (high-frequency stimulation) to quantify long-lasting desynchronization effects. Parameters: *f*_{RR} = 30 Hz (low-frequency stimulation, **A–C**) and *f*_{RR} = 100 Hz (high-frequency stimulation, **D–F**). Stimulation duration is 500 s and *M* = 32. Results were averaged over three network realizations. The acute effects are measured shortly before cessation of stimulation and long-lasting effects 1, 000 s after cessation of stimulation. “ac” and “ll” refer to acute and long-lasting effects, respectively.

Next, we explore the dependence of long-lasting effects on the stimulation amplitude, *A*_{stim} ≤ 1, and the stimulation frequency, *f*_{RR}. We find that weak and moderate stimulation are most efficient for intermediate stimulation frequencies, as shown in Figure 6. In contrast, strong stimulation only entails long-lasting desynchronization if applied at low stimulation frequencies.

**Figure 6**. Acute and long-lasting effects of L/M-RR stimulation. Acute mean weight **(A)** and time-averaged Kuramoto order parameter **(B)** as function of stimulation strength *A*_{stim} and stimulation frequency *f*_{RR}. **(C,D)** Same as **(A,B)** but evaluated 1, 000 s after cessation of stimulation. Parameters: stimulation is delivered for 500 s, *L* = 16, and *M* = 32. “ac” and “ll” refer to acute and long-lasting effects, respectively.

### 3.3. Stimulation With Low Spatial Resolution Performs Better at Low Stimulation Amplitudes

Next, we vary the spatial resolution *M*, scaling the distance between adjacent stimulation sites (Figure 2). To this end, we fix the fraction of simultaneously stimulated subpopulations *L*/*M* and consider different spatial resolutions *M*. Results for low and high stimulation frequencies *f*_{RR} are shown in Figure 7. We find that the spatial resolution strongly impacts weight reduction for weak and strong stimulation. Low resolutions (small *M*) seem to be advantageous for weak stimulation, where smaller amplitudes are sufficient to achieve a pronounced weakening of synapses. In contrast, a large *M* leads to more synaptic weakening for moderate stimulation and a low fraction of simultaneously stimulated sites (Figures 7A,C). Results for strong low and high-frequency stimulation differ significantly. For low stimulation frequencies, we observe pronounced decoupling for all considered spatial resolutions. Contrastingly, for high-frequency stimulation, the mean synaptic weight shortly before stimulation ceases possesses a complex dependence on the fraction of simultaneously activated sites and the spatial resolution. Here, low resolutions result in higher mean weights for small fractions (Figure 7C). For high fractions, however, high spatial resolution may even strengthen synaptic weights, rendering stimulation not suitable for inducing long-lasting desynchronization. This can be seen in Figure 7D, where strong high-frequency stimulation results in large values of the mean synaptic weight. The latter increases with increasing resolution *M*. Hence, stimulation approaches with low spatial resolution may be advantageous for weak and very strong stimulation.

**Figure 7**. Acute effect of L/M-RR stimulation on the mean synaptic weight for different values of total numbers of subpopulation, *M*. Results for the mean synaptic weight shortly before cessation of stimulation for low-frequency **(A,B)** and high-frequency **(C,D)** stimulation are shown. Colors indicate different total numbers of subpopulations, *M*. Columns contain results for different fractions of simultaneously stimulated neuronal subpopulations *L*/*M*. Parameters: Stimulation is delivered for 500 s. “ac” and “ll” refer to acute and long-lasting effects, respectively.

## 4. Discussion

We studied desynchronization in networks of leaky integrate-and-fire (LIF) neurons with spike-timing dependent plasticity (STDP) by Random Reset (RR) stimulation, a decoupling stimulation technique (Kromer and Tass, 2020). RR stimulation was designed to specifically target synaptic weights and drive plastic neuronal networks into the attractor of a stable desynchronized state with weak synaptic connections. This stabilizes desynchronized activity after cessation of stimulation and may lead to long-lasting desynchronization effects (Kromer and Tass, 2020). The original RR stimulation paradigm suggests the delivery of temporally and spatially randomized stimulus patterns. To realize spatial randomization, (Kromer and Tass) delivered each stimulus to 50% of the neurons. These neurons were selected at random, irrespective of their locations in space and their distance to stimulation contacts (Kromer and Tass, 2020). In a DBS setup, however, such a microscopic selection process is not possible due to limited spatial resolution. There, each stimulus affects a finite tissue volume, while the approach of Kromer and Tass (2020) would require that even nearby neurons can be stimulated independently.

Here, we present a version of RR stimulation that copes with limited spatial resolution. Specifically, each stimulus is delivered to a spatially coherent group of *L* out of *M* randomly selected stimulation sites, denoted as L/M-RR stimulation. This setup mimics the delivery of DBS through modern segmented lead electrodes with multiple stimulation contacts (Steigerwald et al., 2016). L/M-RR stimulation does not require single-neuron stimulation as the approach presented in Kromer and Tass (2020), where stimuli were administered to 50% of the neurons that were randomly selected without considering their distribution and location in space.

In order to analyze the performance of L/M-RR stimulation, we apply a recently developed theoretical framework to predict the mean rate of the stimulation-induced reshaping of intrapopulation and interpopulation synapses. The latter connect neurons in the same and different subpopulations, respectively (Kromer and Tass, 2020). We find an excellent agreement between theoretical predictions and numerical simulations for strong stimulation amplitudes, *A*_{stim} ≈ 1, where neuronal spiking follows the stimulus pattern.

L/M-RR stimulation causes parameter-robust long-lasting desynchronization effects. We find stimulation-induced decoupling and related long-lasting desynchronization in the major part of the parameter space spanned by the stimulation frequency and the fraction of simultaneously stimulated subpopulations *L*/*M* (Figure 4). Only for high stimulation frequencies and large fractions, L/M-RR stimulation does not entail long-lasting desynchronization, i.e., it is ineffective. Here, stimuli are delivered at a high pace, which causes high neuronal firing rates and short time lags between post- and presynaptic spiking events. As these time lags become of the order of the STDP decay time for LTP, τ_{+}, synaptic weights start increasing, and stimulation becomes ineffective. This effect causes a qualitative difference between low and high frequency stimulation that is well-described by our theory. Based on our results, we would expect a qualitatively similar outcome for the original RR stimulation protocol presented in Kromer and Tass (2020). There, the percentage of simultaneously stimulated neurons might be comparable to the fraction *L*/*M* of simultaneously stimulated subpopulations in the present paper. However, this percentage was not varied systematically in Kromer and Tass (2020). We further expect our results to be robust with respect to the specific choice of the distribution of interstimulus intervals, Equation (7), as long as the resulting distributions of time lags lead to a sufficiently negative rate of weight change, Equation (8).

Qualitative differences between low- and high-frequency stimulation, observed for strong stimulation, are also present for moderate stimulation amplitudes. While strong high-frequency stimulation is ineffective, at moderate stimulation amplitudes, it leads to long-lasing effects. This is because neurons do not spike in response to stimuli of moderate strength that are delivered shortly after spiking. This leads to longer time lags between post- and presynaptic spikes, which reduces the contribution of LTP to the synaptic weight dynamics and supports stimulation-induced decoupling. This effect leads to a trade-off between strong and weak stimulation. For the former, stimuli hardly impact neuronal spiking, whereas the latter induces short time lags that lead to LTP. As a consequence, L/M-RR stimulation performs best at intermediate stimulation amplitudes. This trade-off differs significantly from those presented in earlier studies reporting optimal performance of coordinated reset (CR) stimulation at intermediate stimulation amplitudes (Lysyansky et al., 2011; Popovych and Tass, 2012; Ebert et al., 2014; Zeitler and Tass, 2015). These studies considered spatial stimulation profiles, where strong stimulation affects larger neuronal populations (Butson and McIntyre, 2008). The latter reduces the decoupling effects of CR stimulation. In contrast, the performance of L/M-RR stimulation becomes worse at strong stimulation amplitudes due to an increased contribution of LTP to the synaptic weight dynamics, due to shorter time lags.

That high stimulation frequencies can lead to qualitatively different synaptic weight dynamics was also observed in our recent study on multisite CR stimulation (Kromer et al., 2020). There, high stimulation frequencies could lead to time lags between post- and presynaptic spikes that are even shorter than the synaptic transmission delay. This led to highly non-linear weight dynamics as a function of the stimulation frequency and the number of stimulation sites, i.e., the spatial resolution (Kromer et al., 2020). In contrast, in the present paper we limited the minimum interstimulus interval to τ_{Λ} > *t*_{d}. This ensures non-overlapping stimuli. Therefore, observed time lags are always larger than the delay time. Nevertheless, we observe synaptic strengthening for strong high-frequency stimulation and large numbers of stimulation sites, which is in accordance with the results presented in Kromer et al. (2020) for high-frequency CR stimulation. By the way, very similar results were obtained for a minimum interstimulus interval τ_{Λ} = 1/250 s.

An interesting question is to which extent our results for strong high-frequency L/M-RR stimulation translate to HF DBS. It is widely observed that symptoms return shortly after cessation of HF DBS therapy in Parkinson's patients (Temperli et al., 2003); hence HF DBS may not stabilize physiological activity after cessation of stimulation. RR stimulation represents a temporally and spatially randomized stimulation approach (Kromer and Tass, 2020). This raises the question of whether long-lasting effects may be caused by sufficient randomization of HF DBS. So far, temporally randomized versions of HF DBS were analyzed in a few experimental studies; however, there are mixed results on its efficacy. Furthermore, to our knowledge, all studies were intraoperative and limited to acute effects during stimulation. In Dorval et al. (2010), a randomized HF DBS pattern was used to treat bradykinesia in PD patients. In their study, interpulse intervals were distributed according to a gamma distribution. The authors hypothesized that regular HF DBS leads to symptom alleviation by reducing the firing irregularities in the basal ganglia; they argue that randomized HF DBS fails to regularize the firing and is therefore inferior to regular HF DBS. Birdno et al. (2012) and Brocker et al. (2013) considered five different irregular types of HF DBS. Two deterministic patterns in which a regular pulse train was either interrupted by periods of silence or periods of high-frequency bursts (Birdno et al., 2012; Brocker et al., (2013); and three randomized pulse trains in which interpulse intervals where distribution according to log-uniform distributions of two different widths (Birdno et al., 2012; Brocker et al., 2013) and according to a bimodal distribution where the inverse interpulse intervals were either shorter or longer than the range of therapeutic frequencies (Birdno et al., 2012). Birdno et al. (2012) found that irregular HF DBS is inferior to the regular one in treating tremor. They argue that pathological activity may propagate during long interpulse intervals. In contrast, Brocker et al. (2013) studied the performance of PD patients in a simple motor task (finger tapping) and reported improved performance during irregular HF DBS. Moreover, their computational model showed that these randomized DBS patterns significantly suppressed beta band power. However, as acknowledged by the authors, the intraoperative setting of their clinical trials limits the duration of the experiment. Therefore, possible effects might not fully develop, and the therapeutic effects of irregular HF DBS might be underestimated.

It is unclear whether experimental results on irregular HF DBS apply to the L/M-RR stimulation protocol suggested in the present study. Furthermore, the L/M-RR stimulation differs in two ways from the studies on irregular DBS. First, L/M-RR stimulation targets pathological connectivity rather than pathological neuronal activity. Long-lasting changes due to randomized HF DBS, however, have yet to be studied. As shown in Figure 3, L/M-RR stimulation needs to be administered for a sufficient amount of time to drive the network into the attractor of the stable desynchronized state. For too short stimulation time, the full potential of long-lasting effects might not be released (Figure 3). Furthermore, as pointed out by our theory for strong stimulation, the stimulation-induced weight dynamics is closely related to the statistics of the interstimulus interval, see Equation (8). Whether the irregular HF DBS protocols studied in Dorval et al. (2010), Birdno et al. (2012), and Brocker et al. (2013) cause a reduction of synaptic weights further depends on the underlying plasticity mechanism. Once the latter has been explored, Equation (8) might be used to predict the potential long-lasting outcome. Second, L/M-RR stimulation combines temporal and spatial randomization, while only temporally randomized HF DBS was considered. In fact, we find that high-frequency L/M-RR stimulation is ineffective for large fractions of simultaneously activated stimulation sites, while low fractions result in pronounced long-lasting effects (Figures 5D–F).

The fraction of simultaneously activated subpopulations *L*/*M* controls the degree of spatial randomization. A fraction of *L*/*M* = 1 corresponds to spatially regular single-site stimulation, while low fractions result in a high degree of spatial randomization. This fraction also impacts the frequency at which individual subpopulations receive stimuli. In particular, neurons receive stimuli at higher rates if this fraction is increased for fixed stimulation frequency. However, our computational results show that this only improves long-lasting effects for low stimulation frequencies (Figures 5A,B). In contrast, increasing this fraction yields worse performance for high-frequency stimulation (Figures 5C,D). Improvement of long-lasting effects at low stimulation frequencies results from the so-called *decoupling through synchrony*, which occurs for asymmetric Hebbian plasticity functions if the distribution of spike times within collective spiking events becomes narrow compared to axonal delays (Lubenov and Siapas, 2008; Knoblauch et al., 2012). In the present paper, we use short stimuli that cause such narrow distributions of spike times during collective spiking events (Kromer and Tass, 2020). This supports decoupling when stimuli are simultaneously delivered to a large number of subpopulations, i.e., *L*/*M* ≈ 1, see Figure 5A. For high stimulation frequencies, however, this effect is balanced by LTP due to short time lags, as described in detail above. Therefore, high-frequency stimulation is ineffective for *L*/*M* ≈ 1. As a result, spatial randomization resulting from intermediate fractions *L*/*M* increases the robustness of long-lasting desynchronization by L/M-RR stimulation with respect to changes in the stimulation frequency.

Of particular interest with respect to a possible implementation of L/M-RR stimulation using a DBS setup is the impact of the number of stimulation sites, represented by the spatial resolution *M*. Commonly used DBS electrodes possess 4–8 stimulation contacts that are arranged equidistantly along the electrode axes (Gielen, 2001; Butson and McIntyre, 2005). Recently, Steigerwald et al. (2019) presented an electrode with up to 32 stimulation contacts that can be activated independently. For strong stimulation, the impact of the number of stimulation sites *M* can be analyzed theoretically. Our theory predicts that the dynamics of intrapopulation weights solely depends on the fraction of simultaneously activated subpopulations *L*/*M*, see Equations (9), (10), and (13). In contrast, the probability of simultaneous activation of different subpopulations, and therefore the dynamics of interpopulation weights depends on (*L* − 1)/(*M* − 1), see Equations (15) and (16). Thus, changes of the spatial resolution, *M*, while keeping the fraction *L*/*M* constant, only affect the dynamics of interpopulation weights. As a consequence, after sufficiently long stimulation different mean synaptic weights can be attained, depending on the combination of signs of the mean rates of weight change ${{J}}^{\mathrm{\text{intra}}}$ and ${{J}}^{\mathrm{\text{inter}}}$, i.e., whether stimulation reduces all weights (${{J}}^{\mathrm{\text{intra}}}<0$ and ${{J}}^{\mathrm{\text{inter}}}<0$), increases all weights (${{J}}^{\mathrm{\text{intra}}}>0$ and ${{J}}^{\mathrm{\text{inter}}}>0$), or whether only one type of weights is reduced while the other one increases. While this determines the outcome of strong stimulation, our simulation results provide evidence that low spatial resolution is favorable for weak stimulation, as shown in Figure 7. We hypothesize that this is because simultaneously stimulated subpopulations consist of bigger localized groups of neurons and such groups are more likely to follow the applied stimulus pattern. Since sufficient decoupling is already achieved for weak stimulation with low spatial resolution, L/M-RR stimulation may be suitable for implementation in a conventional DBS setup. There, weak stimulation reduces the risk of side-effects (Rodriguez-Oroz et al., 2005). Furthermore, low spatial resolution, i.e., a smaller number of stimulation sites, is advantageous because it allows the usage of common DBS electrodes and reduces the time for preparation and parameter tuning. In addition, using a smaller number of simulation sites may increase feasibility given the anatomical and topological constraints of currently available chronically implantable depth electrodes and the related consequences on therapeutic effects and side effects (Volkmann et al., 2006; Krauss et al., 2020).

So far, segmented electrodes were mainly used intraoperatively for directional HF DBS stimulation. During directional DBS, the current flow can be directed in both the vertical and in the horizontal plane by activation of a select number of individual contacts. This results in high spatial selectivity (Contarino et al., 2014; Pollo et al., 2014; Steigerwald et al., 2019). These intraoperative studies suggest that directional DBS may lower the threshold current for beneficial HF DBS (Contarino et al., 2014; Pollo et al., 2014). However, the full potential of segmented electrodes for multisite stimulation methods such as CR or L/M-RR stimulation is yet to be explored. In the present paper, we suggest L/M-RR stimulation as one way to realize the recently developed RR stimulation paradigm by means of segmented electrodes. Our theoretical and computational results indicate that separate stimulation of a large number of neuronal subpopulations, which may correspond to a large number of stimulation contacts, improves long-lasting desynchronization by L/M-RR stimulation for low stimulation frequencies, Figure 7A.

In future studies, we anticipate exploring L/M-RR stimulation in a detailed biophysical model of the subthalamic nucleus, a major target region for therapeutic HF DBS stimulation in Parkinson's patients (Benabid et al., 1991) using more detailed models of stimulation contacts accounting for spatial current profiles. Furthermore, we plan to use large-scale neuronal network models, to study how STN stimulation affects interactions between different nuclei in the basal ganglia region. We hope that our encouraging results motivate experimental and preclinical studies on RR stimulation as a potential treatment for neurological disorders that exploits recently developed segmented electrodes (Steigerwald et al., 2019).

## Data Availability Statement

The original contributions presented in the study are included in the article/supplementary materials, further inquiries can be directed to the corresponding author/s.

## Author Contributions

AKN, JK, and PT conceived the idea, designed the study, interpreted the results, and wrote the manuscript. AKN performed the numerical simulations and analysis. JK developed the theory. All authors contributed to the article and approved the submitted version.

## Conflict of Interest

PT works as consultant for Boston Scientific Neuromodulation. JK and PT are co-inventors of a patent for random reset stimulation owned by Stanford University.

The remaining author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Acknowledgments

We gratefully acknowledge support of this study by Boston Scientific Neuromodulation. Some of the computing for this project was performed on Stanford's Sherlock Computing cluster. We are grateful to Stanford University and the Stanford Research Computing Center for computational resources and support that contributed to these research results.

## References

Abbott, L. F., and Nelson, S. B. (2000). Synaptic plasticity: taming the beast. *Nat. Neurosci*. 3, 1178–1183. doi: 10.1038/81453

Adamchic, I., Hauptmann, C., Barnikol, U. B., Pawelczyk, N., Popovych, O., Barnikol, T. T., et al. (2014). Coordinated reset neuromodulation for Parkinson's disease: proof-of-concept study. *Mov. Disord*. 29, 1679–1684. doi: 10.1002/mds.25923

Alberts, W. W., Wright, E. W., and Feinstein, B. (1969). Cortical potentials and Parkinsonian tremor. *Nature* 221, 670–672. doi: 10.1038/221670a0

Andres, F. G., and Gerloff, C. (1999). Coherence of sequential movements and motor learning. *J. Clin. Neurophysiol*. 16:520. doi: 10.1097/00004691-199911000-00004

Aoki, T., and Aoyagi, T. (2009). Co-evolution of phases and connection strengths in a network of phase oscillators. *Phys. Rev. Lett*. 102:034101. doi: 10.1103/PhysRevLett.102.034101

Ashkan, K., Rogers, P., Bergman, H., and Ughratdar, I. (2017). Insights into the mechanisms of deep brain stimulation. *Nat. Rev. Neurol*. 13:548. doi: 10.1038/nrneurol.2017.105

Benabid, A. L., Pollak, P., Hoffmann, D., Gervason, C., Hommel, M., Perret, J., et al. (1991). Long-term suppression of tremor by chronic stimulation of the ventral intermediate thalamic nucleus. *Lancet* 337, 403–406. doi: 10.1016/0140-6736(91)91175-T

Berner, R., Sawicki, J., and Schöll, E. (2020). Birth and stabilization of phase clusters by multiplexing of adaptive networks. *Phys. Rev. Lett*. 124:088301. doi: 10.1103/PhysRevLett.124.088301

Bi, G.-Q., and Poo, M.-M. (1998). Synaptic modifications in cultured hippocampal neurons: dependence on spike timing, synaptic strength, and postsynaptic cell type. *J. Neurosci*. 18, 10464–10472. doi: 10.1523/JNEUROSCI.18-24-10464.1998

Birdno, M. J., Kuncel, A. M., Dorval, A. D., Turner, D. A., Gross, R. E., and Grill, W. M. (2012). Stimulus features underlying reduced tremor suppression with temporally patterned deep brain stimulation. *J. Neurophysiol*. 107, 364–383. doi: 10.1152/jn.00906.2010

Brocker, D. T., Swan, B. D., Turner, D. A., Gross, R. E., Tatter, S. B., Koop, M. M., et al. (2013). Improved efficacy of temporally non-regular deep brain stimulation in parkinson's disease. *Exp. Neurol*. 239, 60–67. doi: 10.1016/j.expneurol.2012.09.008

Brown, P. (2003). Oscillatory nature of human basal ganglia activity: Relationship to the pathophysiology of Parkinson's disease. *Mov. Disord*. 18, 357–363. doi: 10.1002/mds.10358

Buhlmann, J., Hofmann, L., Tass, P. A., and Hauptmann, C. (2011). Modeling of a segmented electrode for desynchronizing deep brain stimulation. *Front. Neuroeng*. 4:15. doi: 10.3389/fneng.2011.00015

Burkitt, A. N., Meffin, H., and Grayden, D. B. (2004). Spike-timing-dependent plasticity: the relationship to rate-based learning for models with weight dynamics determined by a stable fixed point. *Neural Comput*. 16, 885–940. doi: 10.1162/089976604773135041

Butson, C. R., and McIntyre, C. C. (2005). Role of electrode design on the volume of tissue activated during deep brain stimulation. *J. Neural Eng*. 3:1. doi: 10.1088/1741-2560/3/1/001

Butson, C. R., and McIntyre, C. C. (2008). Current steering to control the volume of tissue activated during deep brain stimulation. *Brain Stimul*. 1, 7–15. doi: 10.1016/j.brs.2007.08.004

Caporale, N., and Dan, Y. (2008). Spike timing-dependent plasticity: a Hebbian learning rule. *Annu. Rev. Neurosci*. 31, 25–46. doi: 10.1146/annurev.neuro.31.060407.125639

Contarino, M. F., Bour, L. J., Bot, M., van den Munckhof, P., Speelman, J. D., Schuurman, P. R., et al. (2012). Tremor-specific neuronal oscillation pattern in dorsal subthalamic nucleus of parkinsonian patients. *Brain. Stimul*. 5, 305–314. doi: 10.1016/j.brs.2011.03.011

Contarino, M. F., Bour, L. J., Verhagen, R., Lourens, M. A., de Bie, R. M., van den Munckhof, P., et al. (2014). Directional steering: a novel approach to deep brain stimulation. *Neurology* 83, 1163–1169. doi: 10.1212/WNL.0000000000000823

Dorval, A. D., Kuncel, A. M., Birdno, M. J., Turner, D. A., and Grill, W. M. (2010). Deep brain stimulation alleviates parkinsonian bradykinesia by regularizing pallidal activity. *J. Neurophysiol*. 104, 911–921. doi: 10.1152/jn.00103.2010

Ebert, M., Hauptmann, C., and Tass, P. A. (2014). Coordinated reset stimulation in a large-scale model of the stn-gpe circuit. *Front. Comput. Neurosci*. 8:154. doi: 10.3389/fncom.2014.00154

Eggermont, J. J., and Tass, P. A. (2015). Maladaptive neural synchrony in tinnitus: origin and restoration. *Front. Neurol*. 6:29. doi: 10.3389/fneur.2015.00029

Gielen, F. L. (2001). *U.S. Patent No. 6,253,109. Washington, DC: U.S. Patent and Trademark Office*. Available online at: https://patents.google.com/patent/US6253109B1/en

Haken, H. (2006). *Brain Dynamics: Synchronization and Activity Patterns in Pulse-Coupled Neural Nets With Delays and Noise*. Berlin; Heidelberg; New York, NY: Springer Science Business Media.

Haken, H. (2012). *Advanced Synergetics: Instability Hierarchies of Self-Organizing Systems and Devices*, Vol. 20. Berlin; Heidelberg; New York, NY; Tokyo: Springer Science Business Media.

Hauptmann, C., Popovych, O., and Tass, P. (2005a). Multisite coordinated delayed feedback for an effective desynchronization of neuronal networks. *Stoch. Dyn*. 5, 307–319. doi: 10.1142/S0219493705001420

Hauptmann, C., Popovych, O., and Tass, P. A. (2005b). Delayed feedback control of synchronization in locally coupled neuronal networks. *Neurocomputing* 65, 759–767. doi: 10.1016/j.neucom.2004.10.072

Hauptmann, C., Popovych, O., and Tass, P. A. (2005c). Effectively desynchronizing deep brain stimulation based on a coordinated delayed feedback stimulation via several sites: a computational study. *Biol. Cybern*. 93, 463–470. doi: 10.1007/s00422-005-0020-1

Hauptmann, C., and Tass, P. (2009). Cumulative and after-effects of short and weak coordinated reset stimulation: a modeling study. *J. Neural Eng*. 6:016004. doi: 10.1088/1741-2560/6/1/016004

Holt, A. B., Wilson, D., Shinn, M., Moehlis, J., and Netoff, T. I. (2016). Phasic burst stimulation: a closed-loop approach to tuning deep brain stimulation parameters for Parkinson's disease. *PLoS Comput. Biol*. 12:e1005011. doi: 10.1371/journal.pcbi.1005011

Jakobs, M., Fomenko, A., Lozano, A. M., and Kiening, K. L. (2019). Cellular, molecular, and clinical mechanisms of action of deep brain stimulation–a systematic review on established indications and outlook on future developments. *EMBO Mol. Med*. 11:e9575. doi: 10.15252/emmm.201809575

Karbowski, J., and Ermentrout, G. B. (2002). Synchrony arising from a balanced synaptic plasticity in a network of heterogeneous neural oscillators. *Phys. Rev. E*. 65:031902. doi: 10.1103/PhysRevE.65.031902

Kempter, R., Gerstner, W., and van Hemmen, J. (1999). Hebbian learning and spiking neurons. *Phys. Rev. E* 59, 4498–4514. doi: 10.1103/PhysRevE.59.4498

Knoblauch, A., Hauser, F., Gewaltig, M.-O., Körner, E., and Palm, G. (2012). Does spike-timing-dependent synaptic plasticity couple or decouple neurons firing in synchrony? *Front. Comput. Neurosci*. 6:55. doi: 10.3389/fncom.2012.00055

Krack, P., Batir, A., Van Blercom, N., Chabardes, S., Fraix, V., Ardouin, C., et al. (2003). Five-year follow-up of bilateral stimulation of the subthalamic nucleus in advanced Parkinson's disease. *N. Engl. J. Med*. 349, 1925–1934. doi: 10.1056/NEJMoa035275

Krauss, J. K., Lipsman, N., Aziz, T., Boutet, A., Brown, P., Chang, J. W., et al. (2020). Technology of deep brain stimulation: current status and future directions. *Nat. Rev. Neurol*. doi: 10.1038/s41582-020-00426-z

Kromer, J. A., Khaledi-Nasab, A., and Tass, P. A. (2020). Impact of number of stimulation sites on long-lasting desynchronization effects of coordinated reset stimulation. *Chaos* 30:083134. doi: 10.1063/5.0015196

Kromer, J. A., and Tass, P. A. (2020). Long-lasting desynchronization by decoupling stimulation. *Phys. Rev. Res*. 2:033101. doi: 10.1103/PhysRevResearch.2.033101

Kühn, A. A., Kupsch, A., Schneider, G.-H., and Brown, P. (2006). Reduction in subthalamic 8-35 hz oscillatory activity correlates with clinical improvement in parkinson's disease. *Eur. J. Neurosci*. 23, 1956–1960. doi: 10.1111/j.1460-9568.2006.04717.x

Kuramoto, Y. (1984). *Chemical Oscillations, Waves, and Turbulence*. Berlin; Heidelberg; New York, NY; Tokyo: Springer. doi: 10.1007/978-3-642-69689-3

Litwin-Kumar, A., and Doiron, B. (2014). Formation and maintenance of neuronal assemblies through synaptic plasticity. *Nat. Commun*. 5, 1–12. doi: 10.1038/ncomms6319

Lourens, M. A., Schwab, B. C., Nirody, J. A., Meijer, H. G., and van Gils, S. A. (2015). Exploiting pallidal plasticity for stimulation in Parkinson's disease. *J. Neural Eng*. 12:026005. doi: 10.1088/1741-2560/12/2/026005

Lozano, A. M., Lipsman, N., Bergman, H., Brown, P., Chabardes, S., Chang, J. W., et al. (2019). Deep brain stimulation: current challenges and future directions. *Nat. Rev. Neurol*. 15, 148–160. doi: 10.1038/s41582-018-0128-2

Lubenov, E. V., and Siapas, A. G. (2008). Decoupling through synchrony in neuronal circuits with propagation delays. *Neuron* 58, 118–131. doi: 10.1016/j.neuron.2008.01.036

Lysyansky, B., Popovych, O. V., and Tass, P. A. (2011). Desynchronizing anti-resonance effect of m: n ON-OFF coordinated reset stimulation. *J. Neural Eng*. 8:036019. doi: 10.1088/1741-2560/8/3/036019

Maistrenko, Y. L., Lysyansky, B., Hauptmann, C., Burylko, O., and Tass, P. A. (2007). Multistability in the Kuramoto model with synaptic plasticity. *Phys. Rev. E* 75:066207. doi: 10.1103/PhysRevE.75.066207

Manos, T., Zeitler, M., and Tass, P. A. (2018). How stimulation frequency and intensity impact on the long-lasting effects of coordinated reset stimulation. *PLoS Comput. Biol*. 14:e1006113. doi: 10.1371/journal.pcbi.1006113

Markram, H., Lübke, J., Frotscher, M., and Sakmann, B. (1997). Regulation of synaptic efficacy by coincidence of postsynaptic APs and EPSPs. *Science* 275, 213–215. doi: 10.1126/science.275.5297.213

Masuda, N., and Kori, H. (2007). Formation of feedforward networks and frequency synchrony by spike-timing-dependent plasticity. *J. Comput. Neurosci*. 22, 327–345. doi: 10.1007/s10827-007-0022-1

Mines, G. R. (1914). On circulating excitations in heart muscle and their possible relation to tachycardia and fibrillation. *Trans. R. Soc. Can*. 8, 43–52.

Mormann, F., Lehnertz, K., David, P., and Elger, C. E. (2000). Mean phase coherence as a measure for phase synchronization and its application to the EEG of epilepsy patients. *Phys. D* 144, 358–369. doi: 10.1016/S0167-2789(00)00087-7

Nini, A., Feingold, A., Slovin, H., and Bergman, H. (1995). Neurons in the globus pallidus do not show correlated activity in the normal monkey, but phase-locked oscillations appear in the MPTP model of parkinsonism. *J. Neurophysiol*. 74, 1800–1805. doi: 10.1152/jn.1995.74.4.1800

Pikovsky, A., Rosenblum, M., and Kurths, J. (2001). *Synchronization: A Universal Concept in Nonlinear Sciences*. Cambridge, UK: Cambridge University Press. doi: 10.1017/CBO9780511755743

Pollo, C., Kaelin-Lang, A., Oertel, M. F., Stieglitz, L., Taub, E., Fuhr, P., et al. (2014). Directional deep brain stimulation: an intraoperative double-blind pilot study. *Brain* 137, 2015–2026. doi: 10.1093/brain/awu102

Popovych, O. V., Hauptmann, C., and Tass, P. A. (2005). Effective desynchronization by nonlinear delayed feedback. *Phys. Rev. Lett*. 94:164102. doi: 10.1103/PhysRevLett.94.164102

Popovych, O. V., Hauptmann, C., and Tass, P. A. (2006a). Control of neuronal synchrony by nonlinear delayed feedback. *Biol. Cybern*. 95, 69–85. doi: 10.1007/s00422-006-0066-8

Popovych, O. V., Hauptmann, C., and Tass, P. A. (2006b). Desynchronization and decoupling of interacting oscillators by nonlinear delayed feedback. *IJBC* 16, 1977–1987. doi: 10.1142/S0218127406015830

Popovych, O. V., Lysyansky, B., Rosenblum, M., Pikovsky, A., and Tass, P. A. (2017a). Pulsatile desynchronizing delayed feedback for closed-loop deep brain stimulation. *PLoS ONE* 12:e0173363. doi: 10.1371/journal.pone.0173363

Popovych, O. V., Lysyansky, B., and Tass, P. A. (2017b). Closed-loop deep brain stimulation by pulsatile delayed feedback with increased gap between pulse phases. *Sci. Rep*. 7:1033. doi: 10.1038/s41598-017-01067-x

Popovych, O. V., and Tass, P. A. (2010). Synchronization control of interacting oscillatory ensembles by mixed nonlinear delayed feedback. *Phys. Rev. E*. 82:026204. doi: 10.1103/PhysRevE.82.026204

Popovych, O. V., and Tass, P. A. (2012). Desynchronizing electrical and sensory coordinated reset neuromodulation. *Front. Hum. Neurosci*. 6:58. doi: 10.3389/fnhum.2012.00058

Pyragas, K., Popovych, O., and Tass, P. (2007). Controlling synchrony in oscillatory networks with a separate stimulation-registration setup. *EPL* 80:40002. doi: 10.1209/0295-5075/80/40002

Rodriguez-Oroz, M. C., Obeso, J., Lang, A., Houeto, J.-L., Pollak, P., Rehncrona, S., et al. (2005). Bilateral deep brain stimulation in Parkinson's disease: a multicentre study with 4 years follow-up. *Brain* 128, 2240–2249. doi: 10.1093/brain/awh571

Röhr, V., Berner, R., Lameu, E. L., Popovych, O. V., and Yanchuk, S. (2019). Frequency cluster formation and slow oscillations in neural populations with plasticity. *PLoS ONE* 14:e0225094. doi: 10.1371/journal.pone.0225094

Rosenblum, M., and Pikovsky, A. (2004a). Delayed feedback control of collective synchrony: an approach to suppression of pathological brain rhythms. *Phys. Rev. E* 70:041904. doi: 10.1103/PhysRevE.70.041904

Rosenblum, M., Pikovsky, A., Kurths, J., Schäfer, C., and Tass, P. A. (2001). “Phase synchronization: from theory to data analysis,” in *Neuro-Informatics and Neural Modelling, Handbook of Biological Physics*, Volume 4, eds S. Gielen and F. Moss (Amsterdam: Elsevier), 279–321. doi: 10.1016/S1383-8121(01)80012-9

Rosenblum, M. G., and Pikovsky, A. S. (2004b). Controlling synchronization in an ensemble of globally coupled oscillators. *Phys. Rev. Lett*. 92:114102. doi: 10.1103/PhysRevLett.92.114102

Seliger, P., Young, S. C., and Tsimring, L. S. (2002). Plasticity and learning in a network of coupled phase oscillators. *Phys. Rev. E* 65:041906. doi: 10.1103/PhysRevE.65.041906

Singer, W. (1990). Search for coherence: a basic principle of cortical self-organization. *Concept Neurosci*. 1, 1–26.

Song, S., Miller, K. D., and Abbott, L. F. (2000). Competitive hebbian learning through spike-timing-dependent synaptic plasticity. *Nat. Neurosci*. 3, 919–926. doi: 10.1038/78829

Steigerwald, F., Matthies, C., and Volkmann, J. (2019). Directional deep brain stimulation. *Neurotherapeutics* 16, 100–104. doi: 10.1007/s13311-018-0667-7

Steigerwald, F., Müller, L., Johannes, S., Matthies, C., and Volkmann, J. (2016). Directional deep brain stimulation of the subthalamic nucleus: a pilot study using a novel neurostimulation device. *Mov. Disord*. 31, 1240–1243. doi: 10.1002/mds.26669

Steigerwald, F., Potter, M., Herzog, J., Pinsker, M., Kopper, F., Mehdorn, H., et al. (2008). Neuronal activity of the human subthalamic nucleus in the parkinsonian and nonparkinsonian state. *J. Neurophysiol*. 100, 2515–2524. doi: 10.1152/jn.90574.2008

Steriade, M., Jones, E. G., and Llinás, R. R. (1990). *Thalamic Oscillations and Signaling*. New York, NY: John Wiley & Sons.

Tass, P., Smirnov, D., Karavaev, A., Barnikol, U., Barnikol, T., Adamchic, I., et al. (2010). The causal relationship between subcortical local field potential oscillations and Parkinsonian resting tremor. *J. Neural Eng*. 7:016009. doi: 10.1088/1741-2560/7/1/016009

Tass, P. A. (1999). *Phase Resetting in Medicine and Biology: Stochastic Modelling and Data Analysis*. Berlin: Springer. doi: 10.1007/978-3-540-38161-7

Tass, P. A. (2001). Effective desynchronization by means of double-pulse phase resetting. *EPL* 53:15. doi: 10.1209/epl/i2001-00117-6

Tass, P. A. (2002). Desynchronization of brain rhythms with soft phase-resetting techniques. *Biol. Cybern*. 87, 102–115. doi: 10.1007/s00422-002-0322-5

Tass, P. A. (2003). A model of desynchronizing deep brain stimulation with a demand-controlled coordinated reset of neural subpopulations. *Biol. Cybern*. 89, 81–88. doi: 10.1007/s00422-003-0425-7

Tass, P. A., and Majtanik, M. (2006). Long-term anti-kindling effects of desynchronizing brain stimulation: a theoretical study. *Biol. Cybern*. 94, 58–66. doi: 10.1007/s00422-005-0028-6

Tass, P. A., Qin, L., Hauptmann, C., Dovero, S., Bezard, E., Boraud, T., et al. (2012). Coordinated reset has sustained aftereffects in Parkinsonian monkeys. *Ann. Neurol*. 72, 816–820. doi: 10.1002/ana.23663

Tass, P. A., Silchenko, A. N., Hauptmann, C., Barnikol, U. B., and Speckmann, E.-J. (2009). Long-lasting desynchronization in rat hippocampal slice induced by coordinated reset stimulation. *Phys. Rev. E* 80:011902. doi: 10.1103/PhysRevE.80.011902

Temel, Y., Kessels, A., Tan, S., Topdag, A., Boon, P., and Visser-Vandewalle, V. (2006). Behavioural changes after bilateral subthalamic stimulation in advanced Parkinson disease: a systematic review. *Parkinsonism Relat. Disord*. 12, 265–272. doi: 10.1016/j.parkreldis.2006.01.004

Temperli, P., Ghika, J., Villemure, J.-G., Burkhard, P., Bogousslavsky, J., and Vingerhoets, F. (2003). How do Parkinsonian signs return after discontinuation of subthalamic DBS? *Neurology* 60, 78–81. doi: 10.1212/WNL.60.1.78

Volkmann, J., Moro, E., and Pahwa, R. (2006). Basic algorithms for the programming of deep brain stimulation in parkinson's disease. *Mov. Disord*. 21, S284–S289. doi: 10.1002/mds.20961

Wang, J., Nebeck, S., Muralidhara, A., Johnson, M. D., Vitek, J. L., and Baker, K. B. (2016). Coordinated reset deep brain stimulation of subthalamic nucleus produces long-lasting, dose-dependent motor improvements in the 1-methyl-4-phenyl-1,2,3,6-tetrahydropyridine non-human primate model of Parkinsonism. *Brain. Stimul*. 9, 609–617. doi: 10.1016/j.brs.2016.03.014

Warman, E., and Durand, D. (1989). “Desynchronization of epileptiform activity by phase resetting,” in *Images of the Twenty-First Century. Proceedings of the Annual International Engineering in Medicine and Biology Society* (Seattle, WA: IEEE), 1286–1287. doi: 10.1109/IEMBS.1989.96197

Weinberger, M., Mahant, N., Hutchison, W. D., Lozano, A. M., Moro, E., Hodaie, M., et al. (2006). Beta oscillatory activity in the subthalamic nucleus and its relation to dopaminergic response in Parkinson's disease. *J. Neurophysiol*. 96, 3248–3256. doi: 10.1152/jn.00697.2006

Wilson, C. J., Beverlin, I. I. B, and Netoff, T. (2011). Chaotic desynchronization as the therapeutic mechanism of deep brain stimulation. *Front. Syst. Neurosci*., 5:50. doi: 10.3389/fnsys.2011.00050

Winfree, A. T. (1977). Phase control of neural pacemakers. *Science* 197, 761–763. doi: 10.1126/science.887919

Winfree, A. T. (1980). *The Geometry of Biological Time*. New York, NY: Springer. doi: 10.1007/978-3-662-22492-2

Winfree, A. T. (2001). *The Geometry of Biological Time*, Volume 12. New York, NY: Springer Science & Business Media.

Yanchuk, S., Berner, R., and Schöll, E. (2020). “Frequency clusters in adaptive networks,” in *2020 European Control Conference (ECC)* (Saint Petersburg: IEEE), 313–316. doi: 10.23919/ECC51009.2020.9143609

Zanette, D. H., and Mikhailov, A. S. (2004). Dynamical clustering in oscillator ensembles with time-dependent interactions. *EPL* 65, 465–471. doi: 10.1209/epl/i2003-10124-1

Zeitler, M., and Tass, P. A. (2015). Augmented brain function by coordinated reset stimulation with slowly varying sequences. *Front. Syst. Neurosci*. 9:49. doi: 10.3389/fnsys.2015.00049

Keywords: random reset stimulation, spike-timing dependent plasticity (STDP), desynchronization, segmented electrodes, long-lasting effects

Citation: Khaledi-Nasab A, Kromer JA and Tass PA (2021) Long-Lasting Desynchronization of Plastic Neural Networks by Random Reset Stimulation. *Front. Physiol.* 11:622620. doi: 10.3389/fphys.2020.622620

Received: 28 October 2020; Accepted: 23 December 2020;

Published: 05 February 2021.

Edited by:

Plamen Ch. Ivanov, Boston University, United StatesReviewed by:

Serhiy Yanchuk, Technical University of Berlin, GermanySonya Bahar, University of Missouri System, United States

Copyright © 2021 Khaledi-Nasab, Kromer and Tass. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Ali Khaledi-Nasab, khaledi@stanford.edu