Long-Lasting Desynchronization of Plastic Neural Networks by Random Reset Stimulation

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.

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, 1977Winfree, , 1980Warman and Durand, 1989;Tass, 1999), followed by the development of composite single-channel (Tass, 2001(Tass, , 2002Zhai et al., 2005) and multichannel 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., 2005Popovych et al., , 2006aPyragas 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; 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 , 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 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;, 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 . 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,  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 . 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 . For our analysis, we use a combination of theoretical predictions and simulations of networks of leaky integrateand-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.

Neuronal Network Model
We consider a network of N conductance-based LIF neurons with STDP previously presented in . 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 (1) C i is the membrane capacitance, v rest the resting potential, g leak = 0.02 mS/cm 2 the leakage conductance, I syn i (t) the synaptic input current, S i (t) the stimulation current, and I noise i (t) the noisy input current.
Spiking occurs when v i crosses a dynamic threshold potential v th i given by Here, τ th is the threshold time constant and v th rest 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 th i → v th spike . Throughout the paper, we use the same parameters as in : v rest = −38 mV, v th rest = −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 syn i (t) to neuron i, is given by where g syn i 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 j l 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 noise i 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 noise i (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.

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., 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;) 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 longterm 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 .

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 . ρ 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.

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.
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 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 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 e = A stim µ/ν e ; separated by a gap of 0.2 ms, the inhibitory pulse has a duration of ν i = 1.5 ms and an amplitude of 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 1 M , and for the second one we select one of the remaining electrodes with the probability of 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).
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 e = A stim µ/ν e ; separated by a gap of 0.2 ms, the inhibitory pulse has a duration of ν i = 1.5 ms and an amplitude of 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.

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 . In the following, we briefly state the main steps and expand their results to our L/M-RR stimulation protocol.
The approach of  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 . 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.  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 . 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 . 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  and consider possible pairings of a presynaptic (postsynaptic) spike that is triggered by the nth 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 nth stimulus.
In  results for p ij (S, ξ ) for either case were derived. In the first case, the distribution of interstimulus intervals is given by Here the kth 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

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.0 ms/τ . Convolutions are calculated using the python method numpy.convolve of numpy version 1.16.2.

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 . 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  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.
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.

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 intra , and interpopulation synapses, J 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 intra < 0, and J inter < 0, indicate a weakening of corresponding synapses during stimulation.
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 intraand 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 intra < 0 and J inter < 0 where L/M-RR stimulation decouples all neurons, (ii) a region where only intrapopulation synapses are weakened (J intra < 0 and J inter > 0), and (iii) a region where all synapses are strengthened (J intra > 0 and J 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.

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. . We used λ(t) = δ(t) for theoretical predictions. "ac" and "ll" refer to acute and long-lasting effects, respectively. 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.
frequencies, moderate stimulation yields pronounced longlasting 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. 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.

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.

DISCUSSION
We studied desynchronization in networks of leaky integrateand-fire (LIF) neurons with spike-timing dependent plasticity (STDP) by Random Reset (RR) stimulation, a decoupling stimulation technique . 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 . 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 . 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  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 singleneuron stimulation as the approach presented in , 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 . 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 . 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 . 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 postand 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 . 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 . 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  for highfrequency 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 . 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 loguniform 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 highfrequency 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 . 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 longlasting 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 intra and J inter , i.e., whether stimulation reduces all weights (J intra < 0 and J inter < 0), increases all weights (J intra > 0 and J 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.