Long-lasting desynchronization by decoupling stimulation

Several brain disorders are characterized by abnormally strong synchronization of neuronal activity. In Parkinson’s patients, permanent high-frequency deep brain stimulation is used to suppress symptoms. To speciﬁcally counteract synchronized neuronal activity with a substantially reduced amount of stimulation current, theory-based desynchronizing stimulation techniques were developed, e.g., coordinated reset stimulation. Desynchronizing stimulation may shift adaptive networks from attractors with strong synchronization and strong synaptic coupling to attractors with weak synchronization and weak coupling. This is to cause stimulation effects that persist after cessation of stimulation. Corresponding preclinical and clinical studies reported long-lasting desynchronization and related symptom relief. However, desynchronizing stimulation requires parameters to be adapted to characteristics of the synchronized neuronal activity. Furthermore, desynchronization does not guarantee long-lasting change of network activity. We here present a qualitatively different approach to induce long-lasting, sustained changes of neuronal network dynamics: decoupling stimulation. Instead of primarily desynchronizing neuronal activity, decoupling stimulation employs synaptic plasticity mechanisms to speciﬁcally decouple neuronal networks. In this way, neuronal networks get robustly shifted to attractors with desynchronized neuronal activity. We present a theoretical framework that explains how neuronal responses to single stimuli as well as to spatiotemporally correlated stimulus sequences impact on network connectivity. This provides a theoretical base for designing effectively decoupling stimulation protocols. To overcome limitations of primarily desynchronizing stimulation, we present a random reset stimulation protocol, which uses spatiotemporal stimulus randomization to effectively decouple neurons. Theoretical predictions of random reset-induced decoupling as opposed to desynchronization-induced decoupling achieved by coordinated reset stimulation are compared to simulations of networks of integrate-and-ﬁre neurons with spike-timing-dependent plasticity. Decoupling and related long-lasting desynchronization effects achieved by random reset stimulation are more robust with respect to parameter changes than those for coordinated reset stimulation. For both random reset and coordinated reset stimulation, simulation results and theoretically predicted decoupling rates show good quantitative agreement for sufﬁciently strong stimulation amplitudes. Intriguingly, single stimulus-related mechanisms may have a stronger decoupling impact than stimulus sequence-related mechanisms. We discuss scope and limitations of our decoupling approach for different types of synaptic plasticity and its application to deep brain stimulation.


I. INTRODUCTION
Synchronization processes abound in various fields of the natural sciences [1]. Prominent examples include the synchronous flashing of fireflies, the synchronization of circadian rhythms with the daylight cycle, and the synchronized action of heart cells [1]. In many systems synchrony of rhythmic activity is critical for optimum functioning. In other systems, however, abnormally strong synchrony is associated with impaired function. For instance, in the brain pathological synchronization of neuronal activity is related to several brain disorders, such as Parkinson's disease [2] and epilepsy [3]. Experimentally observed synchronization phenomena motivated theoretical studies on coupled oscillators and led to a deep theoretical understanding of synchronization processes [4][5][6]. The state of an oscillator is typically described by its phase on a limit cycle. Networks of such phase oscillators have been used to study synchronization of populations of active elements, e.g., metronomes, heart cells, or neurons [4][5][6][7][8][9][10][11].
Possible applications in biology and medicine led to extensive research on stimulation techniques that allow for actively synchronizing or desynchronizing the collective dynamics of coupled oscillators. Early studies were devoted to inducing temporary desynchronization by delivering a single stimulus pulse of moderate amplitude at a vulnerable phase of the collective rhythm [12][13][14]. Repetitive or demand-controlled stimulus delivery can counteract pronounced synchronization. Stimulation-induced desynchronization has been demonstrated computationally [14], in experiments on coupled electrochemical oscillators [14,15], and in patients with essential tremor stimulated through implanted electrodes [16]. Stimuli with sufficiently high amplitude can cause a phase reset, so that the phases of simultaneously stimulated oscillators align regardless of their individual phases before stimulus administration [13,14]. Thus, a phase resetting stimulus provides control over the oscillators' phases and can be used to restart a network of oscillators independently of its initial dynamic state. A second, weaker pulse following a phaseresetting pulse may desynchronize a neuronal network by hitting it in a vulnerable state [17]. Periodic stimulation may also cause desynchronization. The effect of periodic stimulation on a population of synchronized oscillators depends on the stimulation frequency. While synchronous oscillations become phase locked to the stimulation for certain stimulation frequencies, other frequencies desynchronize the oscillators. This effect is known as chaotic desynchronization [18,19]. Suitable frequencies were predicted from estimated phase response curves of synchronous oscillations. Furthermore, real time estimates of phase response curves can be used to time stimulus delivery such that phase differences between individual oscillators grow. A corresponding close-loop setup for invasive deep brain stimulation was suggested in Ref. [20]. In contrast, coordinated reset (CR) stimulation causes desynchronization without phase-dependent stimulus delivery and without synchronizing stimuli [21]. During CR stimulation, phase resetting stimuli are administered in a coordinated way to different subpopulations of oscillators. Appropriate time shifts between stimuli induce phase shifts between rhythmically active subpopulations and break the dominant rhythm of synchronized oscillators [21]. Other desynchronizing stimulation methods use linear and nonlinear delayed feedback as smooth stimuli or envelop of pulse trains, delivered either permanently or on demand [22][23][24][25][26][27], or specifically target vulnerable nodes and perturb their dynamics such that they destabilize and finally break the global synchronous rhythm [28]. Theoretical studies on stimulation by means of a common white Gaussian noise input, on the other hand, showed that such input can either synchronize or desynchronize a population of oscillators [29].
All those studies were devoted to networks of oscillators with fixed coupling strengths. However, in various systems adaptive coupling mechanisms allow for a reshaping of network connectivity in response to appropriate stimuli or environmental changes [30,31]. Networks with adaptive coupling strengths can be found in different scientific fields and include transport networks, wireless communication networks, social networks or neuronal networks. Often reshaping of network connectivity is realized by local rules that alter weights of individual connections according to the local dynamics [32]. For instance, in the tubular network formed by the slime mold Physarum polycephalum the thickness of links is regulated according to the flow through the network [33], or in wireless communication networks where paths from initial sender to final receiver have to be established by self-organization of nodes in a way that balances the needs for network connectivity, robustness against perturbations and data-traffic performance [34]. In the nervous system, a prominent adaptive coupling mechanism is spike-timing-dependent plasticity (STDP). It modifies synaptic coupling strengths according to the relative timing of pre-and postsynaptic discharges, or spikes [35][36][37][38]. In several brain regions STDP causes a strengthening of excitatory synaptic connections if the postsynaptic spike follows shortly after the presynaptic one and a weakening in the opposite case. Network adaptation due to STDP has been found to stabilize activity patterns of neuronal networks [39] and is believed to play a crucial role in memory formation and maintenance [40]. In general, adaptive connectivity may increase the self-organizational complexity and can lead to coexisting states with distinct collective dynamics. In adaptive networks of oscillators, these include synchronized, desynchronized and cluster states [41][42][43][44][45][46][47][48].
Of particular interest is the response of plastic networks to external stimuli. In neuronal networks with STDP, stimulation of subpopulations may lead to the formation of strongly connected neuronal clusters encoding information about earlier stimuli [40,49]. Other, studies revealed that responses of adaptive networks to external stimuli may be complex and even counterintuitive. For instance, uncorrelated white Gaussian noise stimulation, a perfect means for desynchronization in networks with fixed coupling strength, may strengthen synaptic connections in neuronal networks with STDP, thereby supporting synchrony [50]. In contrast, as shown computationally, desynchronization by CR stimulation of such networks may result in a desynchronization-induced synaptic decoupling [51]. In multistable networks, desynchronizationinduced decoupling may drive the synchronized network into the basin of attraction of stable desynchronized states. This causes long-lasting desynchronization, which persists after cessation of stimulation [51]. Multistable neuronal networks have been used to model disease-related qualitative changes of neuronal activity. Parkinson's disease, for instance, relates to excessive neuronal synchrony in the subthalamic nucleus (STN) [2]. Corresponding brain regions have therefore been described as a multistable neuronal network in which a pathological state with synchronized activity coexists with a physiological state with desynchronized activity [30,51,52]. This led to the hypothesis that appropriate stimulation may induce long-lasting therapeutic effects, by driving the system into the attractor of a physiological state [51]. Long-lasting, cumulative desynchronization as well as therapeutic effects of CR stimulation were demonstrated in animal experiments as well as clinical proof-of-concept studies [53][54][55][56].
However, desynchronizing stimulation is an indirect means which does not necessarily lead to decoupling [52]. Furthermore, effective desynchronization requires stimulation parameters to be adapted to characteristics of the collective synchronized rhythm, especially its dominant frequency [22,52,57]. The latter requirement may further complicate matters if there is not just one dominant frequency as, for instance, in Parkinson's disease and other brain disorders [58,59].
What might be an alternative to desynchronizing stimulation? Previous work on CR stimulation has shown that spatiotemporal stimulus patterns may not only desynchronize neurons, but also reshape network connectivity and thereby allow for long-lasting therapeutic effects [51,[53][54][55][56]. As stimulation electrodes improve, more complex spatiotemporal stimulus patterns can be delivered to target brain regions. For instance, segmented multisite electrodes enable to target distinct neuronal subpopulations of the stimulated brain region using directional steering [60,61]. Appropriately designed spatiotemporal stimulus patterns may be able to induce 033101-2 long-lasting therapeutic effects during brain stimulation and overcome limitations of primarily desynchronizing stimulation.
We present a stimulation approach that directly weakens or strengthens the adaptive connections. To this end, we deliver spatiotemporal patterns of phase-resetting stimuli to generate activity patterns that trigger a desired reshaping of network connectivity. In the first part of the paper, we present a theoretical framework, Sec. III, that relates reshaping of network connectivity to the response characteristics of neurons to a single stimulus and spatiotemporal correlations in the sequence of applied stimuli, Sec. IV. Our approach is quite general and broadly applicable in that it allows for various modifications of network structure, including strengthening or weakening of synaptic connections between individual neurons or certain groups of neurons, Sec. VI. However, we here emphasize the decoupling aspect because of its significance to neuroscience and neurology. Based on our theoretical results, we design a random reset (RR) protocol in Sec. V A. Delivering stimuli at random times to random locations leads to parameter-robust decoupling in networks with stable desynchronized states. We compare theoretical predictions on decoupling by RR and CR stimulation for different types of STDP in Sec. VII. Finally, we perform simulations for networks of leaky integrate-andfire neurons with STDP and compare simulation results with our theoretical predictions, Sec. VIII. Our results provide a novel theoretical framework, e.g., for designing deep brain stimulation protocols for preclinical and clinical use.

II. SPIKE-TIMING-DEPENDENT PLASTICITY
We consider network adaptation as a result of STDP. In the main text, we focus on canonical Hebbian STDP [36,38], which was used in previous studies on desynchronization by CR stimulation [51,52]. Other forms of STDP, such as symmetric Hebbian and anti-Hebbian STDP, will be discussed in the Appendix. Canonical Hebbian STDP reinforces causual relations ships between post-and presynaptic spiking activity. It is realized by updating the weights of individual synaptic connections w i→ j when spikes of the presynaptic neuron i or postsynaptic neuron j arrive at the synapse. The update is given by the STDP function [63]: W (t ) depends on the time lag t = t post − t pre between latest spike arrival times of post-and presynaptic spikes. Here we consider only axonal delay t d as those are typically much longer than dendritic ones. Thus, a presynaptic spike at time t i arrives at time t pre = t i + t d , while postsynaptic spikes arrive instantaneously, see schematics in Figs. 1(a) and 1(c). The parameter δ scales the update of individual synaptic weights per spike. The STDP decay time τ + is usually in the range of tens of milliseconds [36]. τ R scales the asymmetry in STDP decay times. β > 0 scales the ratio of overall depression 0 −∞ dt W (t ) < 0 to potentiation  (1). Weights of ingoing synapses are updated at spike times of the postsynaptic neuron while outgoing weights are updated at the arrival times of the presynaptic spikes at the synapse [62]. Time lags t post − t pre considered for weight updates are shown by solid (ingoing) and dashed (outgoing) red bars, respectively. (d) One possible realization of a stimulation sequence S causing spiking responses of neurons i and j shown in panel (a). Stimuli are delivered to the four subpopulations (colors) shown in (b). The stimulation sequence illustrated in (d) is generated using the CR stimulation protocol described in Sec. V B.

III. THEORETICAL DESCRIPTION OF STIMULATION-INDUCED SYNAPTIC WEIGHT DYNAMICS
We analyze the weight dynamics resulting from the delivery of spatiotemporal sequences of stimuli to a neuronal network with STDP. To this end, we consider sequences of charge-balanced pulsatile stimuli delivered to neuronal subpopulations. In neuroscience and medicine electrical stimuli are typically charge-balanced to avoid tissue damage [64]. We characterize pulsatile stimulation protocols by means of a stimulation sequence S. Formally, S is a sequence of pairs including a stimulation time s k and a location P k , i.e., Here the index k counts individual stimuli. P k specifies which neuronal subpopulation receives the stimulus that is delivered at time s k . We will consider P k as the set of neurons that receive the stimulus. By assuming that all neurons in P k receive the same stimulus, we neglect the influence of the distance between the neuron and the stimulation contact. One possible stimulation sequence is illustrated in Figs. 1(b) and 1(d).
To study stimulation-induced reshaping of network connectivity, we first consider a single synaptic weight w i→ j between presynaptic neuron i and postsynaptic neuron j. Spiking of either neuron is characterized by their spike trains x l (t ) = n δ(t − t l n ), with l = i, j. The sum runs over all spike times t l n of the neuron. δ(t ) is the Dirac delta distribution.

033101-3
When a phase-resetting stimulus is delivered to neuron l at time s k the neuron is assumed to fire a single spike in direct response to that stimulus. We denote the timing of that spike by t l . As the neuronal response is in general stochastic, we introduce the distribution of spike response times λ l (s k , t l ). It characterizes the probability distribution of spike times t l of neuron l in response to the stimulus delivered at time s k . A phase reset is characterized by the fact that the distribution of spike response times does not depend on the precise timing of the stimulus in relation to the neuron's dynamical state, i.e., λ l (s k , t l ) = λ l (t l − s k ). In the present paper, we further restrict ourselves to the case where the distribution of spike response times is sufficiently homogeneous among the neurons, i.e., λ l (t l − s k ) ≈ λ(t l − s k ).
We consider the case where neurons mainly spike in response to stimulus administration. This case will be referred to as stimulation-controlled spiking in the following. Formally, this requires that stimuli are delivered faster than the neurons intrinsic spiking frequency, i.e., r l S k S < 1, and that the spiking response occurs before the next stimulus is received, i.e., λ l decays fast compared to interstimulus intervals. Here r l is the lth neuron's mean firing rate, i.e., r l = x l (t ), in the absence of stimulation. x refers to time averaging and x S to averaging over the sequence S. For stimulation-controlled spiking, we can approximate the stimulated neurons' spike trains as (3) k := t l k − s k is the time lag between stimulus delivery and evoked spiking response of neuron l. k is distributed according to the distribution of spike response times λ( k ).
Next, we consider the mean rate of weight change J i j , which is given by [62] Here T is the length of a time window starting at t and the sum runs over all pairs of post-and presynaptic spike times to be considered for weight updates in that time window. Note that the ensemble of considered spike pairs depends on the STDP scheme [65]. Using Eq. (3) for pre-and postsynaptic neurons in Eq. (4) enables us to calculate J i j for a given stimulation sequence S. The dynamics of a synaptic weight w i→ j resulting from stimulation with a certain stimulation protocol can be quantified by the expectation value J i j (t, T ) which can be obtained by ensemble averaging over all possible realizations of stimulation sequences S for the respective protocol. For stationary stimulation protocols, where the spatiotemporal correlations in the sequence S do not change over time, that expectation value is independent of time t, i.e., J i j (t, T ) → J i j (T ) . Furthermore, it saturates as a function of T , i.e., J i j (T ) → J ∞ i j for time windows that are long enough to capture all temporal characteristics of the stimulation protocol. From Eq. (4), we find Here N i j (t ) is the average number of time lags t, between post-and presynaptic spike times, per unit time that contribute to weight changes.

IV. STIMULUS-AND SEQUENCE-INDUCED RESHAPING OF NETWORK CONNECTIVITY
Analyzing Eq. (5), we find two fundamental mechanisms of stimulation-induced reshaping of network connectivity: stimulus-induced reshaping and sequence-induced reshaping. These two mechanisms result from two distinct contributions to time lags t between post-and presynaptic spikes. The first one results from the stochastic responses of the neurons to individual stimuli and will be denoted as . It is quantified by the shape of the distribution of spike response times λ. The second one results from interstimulus intervals between stimuli delivered to the pre-and postsynaptic neuron, respectively. The latter will be denoted as S. Using t = S + , we find where ( ) = dt λ(t )λ(t + ) is the distribution of time lags if post-and presynaptic neurons spike in response to the same stimulus. p i j (S| ) is the probability distribution that neurons i and j receive stimuli with an interstimulus interval S and the time lag between the resulting pair of spikes contributes to the weight change in Eq. (5) for a given . S l refers to the mean time interval between subsequent stimuli that arrive at neuron l. From Eqs. (5) and (6), we can infer the contributions of stimulus-and sequence-induced reshaping to the synaptic weight dynamics. Stimulus-induced reshaping is represented by the outer integral in Eq. (6). The latter runs over all possible realizations of time lags resulting from the stochastic spiking responses of post-and presynaptic neurons to a single stimulus. The distribution of those time lags is given by ( ) and does not depend on the realization of interstimulus intervals. Therefore, reshaping of network connectivity due to stimulus-induced reshaping is solely controlled by the shape of λ. The latter may be manipulated by adjusting stimulus parameters, e.g., the stimulus wave form or amplitude.
The contribution of sequence-induced reshaping to the synaptic weight dynamics is represented by the inner integral in Eq. (6). The key quantity characterizing sequence-induced reshaping is p i j (S| ). It captures two aspects: (i) the distribution of interstimulus intervals between stimuli delivered to neuron j and i, respectively, and (ii) whether time lags resulting from the triggered spiking responses contribute to weight updates. The former aspect characterizes spatiotemporal correlations in the stimulation sequence S. The latter aspect characterizes the considered STDP scheme. In the nearest neighbor scheme considered throughout the paper, postsynaptic spikes are paired with the latest presynaptic spike arrival and vice versa. Another choice may be an all-to-all scheme in which time lags between all post-and presynaptic spikes contribute to weight updates, see Discussion. Next, we consider stimulus and sequence-induced reshaping for two special cases.
In the first case, we assume that pre-and postsynaptic neurons are stimulated simultaneously and stimulation is slow compared to the STDP decay times τ + and τ + τ R , i.e., p i j (S| ) ≈ δ(S). We denote the frequency at which stimuli are delivered as f slow stim . In this case, weight changes are dominated by stimulus-induced reshaping, and J ∞ i j can be approximated by With Eq. (7) we calculate J ∞ i j SIR for a Gaussian distribution of spike response times Note that J ∞ i j SIR does not depend on μ λ as only (t ) enters Eq. (7). We will therefore set μ λ = 0 in the following. In Fig. 2, J ∞ i j SIR is shown for STDP functions with balanced overall depression and potentiation, β = 1. We find that the width of the distribution of spike response times σ λ in units of the synaptic transmission delay and the asymmetry of STDP decay times τ R are critical for the net effect of SIR on synaptic weights. In particular, sharp distributions lead to a weakening of the synaptic weight, while strengthening is observed for asymmetric STDP decay times, τ R > 1, and broad distribution, σ λ > t d . Thus, this special case describes the previously observed decoupling through strongly synchronized spiking as studied in Refs. [66,67].
As the second special case, we consider stimulation with independent Poisson spike trains delivered at a stimulation frequency f P stim to pre-and postsynaptic neurons, respectively. Then, we find p(S) = f P stim exp(− f P stim |S − t d |) and can set λ( ) = δ( ). Assuming stimulation-controlled spiking, Eq. (3), we obtain independent Poisson spike trains with firing rate f P stim for both, the pre-and the postsynaptic neuron. Using Eqs. (6) and (5) we find P < 0 guarantees that STDP leads to decoupling of neurons with uncorrelated Poisson spike trains. It is therefore often considered as an indicator for the existence of a stable desynchronized state [68].

V. RESHAPING OF NETWORK CONNECTIVITY DURING RANDOM AND COORDINATED RESET STIMULATION
We study the contributions of stimulus-and sequenceinduced reshaping in more detail for two specific stimulation protocols. The first protocol intends to avoid unfavorable resonant relationships between stimulation frequency and neuronal firing rates. To this end, we obviate spatial and time periodicities in the stimulation pattern and rather deliver stimuli to randomly chosen neuronal subpopulations of finite size at random times. We will denote such a protocol as random reset (RR). We hypothesize that spatiotemporal random stimulation may reduce spatial-temporal correlations between stimulus administrations and, in turn, induce a statistic of time lags between pre-and postsynaptic spikes that is comparable to that obtained for independent Poisson spiking (see previous paragraph). For comparison, as second protocol we consider coordinated reset (CR) stimulation. By design, CR stimulation is a desynchronizing stimulation protocol [21] which is characterized by distinct spatiotemporal correlations in the stimulation sequence.
For both protocols, we study stimulus-and sequenceinduced reshaping of network connectivity. To this end, we calculate N i j (t ) and J ∞ i j according to Eqs. (6) and (5), respectively.

A. Random reset stimulation
Random reset stimulation means to deliver stimuli at random stimulation times s k to randomly selected neuronal subpopulations P k , thereby minimizing spatiotemporal correlations in the stimulation sequence. Interstimulus intervals S k = s k+1 − s k are drawn from a shifted exponential distribution Here T accounts for a finite duration of a single stimulus. (x) is the Heaviside step function and T RR determines the time scale of interstimulus intervals. T RR sets the mean stimulation frequency to f RR stim ≈ (T + T RR ) −1 . In the limit of short stimuli T → 0, Eq. (10) yields a Poisson spike train for the sequence of stimulation times. Each stimulus is delivered to a neuronal subpopulation P k that contains N e = N/2 neurons. Here N is the total number of neurons and N e is the number of neurons per stimulated subpopulation P. For illustration, 033101-5 subpopulations contain neurons with subsequent indices while accounting for periodic boundary conditions. A raster plot of spiking activity resulting from a representative RR stimulation sequence under the assumptions of stimulation-controlled spiking is shown in Fig. 3

(a).
We calculate the distribution of time lags per unit time N i j (t ) during RR stimulation. This distribution will be denoted as N RR i j (t ) in the following. To this end, we consider a presynaptic spike in response to the nth stimulus. According to the nearest neighbor STDP scheme, Sec. II, the arrival time of this presynaptic spike at the synapse is paired with the latest and the next postsynaptic spikes. Using the stimulationcontrolled spiking approximation introduced above, these spikes result from the latest and the next stimuli delivered to the postsynaptic neuron. We distinguish between two cases: (a) only the presynaptic neuron receives the nth stimulus and (b) both post-and presynaptic neurons receive the nth stimulus. For both cases, we calculated the probability distribution p RR i j (S| ) by considering intervals S between stimuli delivered to the post-and presynaptic neuron. In case (a), It results from the distribution of intervals between stimulus deliveries to a single neuron F RR (S). Assuming that delays are shorter than the minimum interstimulus interval, t d < T , we find with The kth summand in Eq. (12) contains the (k − 1)th convolution of ρ(S). The zeroth convolution refers to ρ itself.
The prefactors are the probabilities that the next (latest) stimulus received by the postsynaptic neuron is the (n ± k)th stimulus. In case (b), p RR i j (S| ) does depend on as time lags between different spike pairs are considered when the presynaptic spike arrives before or after the postsynaptic one at the synaptic terminal, i.e., depending on whether > t d or < t d , respectively. This yields Contributions p RR,a i j and p RR,b i j to p RR i j (S| ) are weighted by the probability that both, neuron i and j receive stimulus n. This probability is given by (N e − | j − i| N )/N e , for | j − i| N < N e and 0 otherwise. Here | j − i| N refers to the absolute distance between indices i and j under consideration of periodic boundary conditions. Using the result in Eq. (6) yields the average number of time lags per unit time for RR stimulation Results for N RR i j (t ) are shown in Fig. 3(b). The probability that pre-and postsynaptic neurons i and j, respectively, receive stimuli simultaneously is critical for the shape of N RR i j (t ). In Fig. 3 we consider the two extreme cases of maximum and minimum probability. Those will be denoted as P RR max and P RR min synapses in the following. For RR stimulation, P RR max synapses are synapses that connect neurons with adjacent indices, i.e., | j − i| N = 1, while P RR min synapses are connecting neurons with faraway indices | j − i| N N e . The two cases result in qualitatively different shapes of N RR i j (t ). For P RR max synapses, we find a pronounced peak centered at t = 0. This results from a strong contribution of stimulus-induced reshaping. In contrast, for P RR min synapses, N RR i j (t ) approaches a combination of two exponentials as T decreases. For T → 0 it converges to r 2 RR exp (−r RR |t − t d |), expressing the lack of temporal correlations in the stimulation sequence. Note that N i j (t ) would possess a similar shape if pre-and postsynaptic spike trains were independent Poisson spike trains with firing rates r i = r j = r RR , with r RR = 1/ N N e (T RR + T ), see Eq. (9).
which yields an estimate of the rate of weight change for a synapse between neurons i and j. Results for J ∞ i j for P RR max and P RR min synapses are shown in Fig. 3(c) (top and bottom), respectively. For P RR max synapses, J ∞ i j becomes positive for σ λ > t d indicating strengthening of synaptic connections. Otherwise, we find that RR stimulation results in a weakening of synaptic connections. The transition from weakening to strengthening is a result of the asymmetry in STDP decay times which is expressed by τ R = 1 in the STDP function, Eq. (1). In contrast, we find solely negative J ∞ i j for P RR min synapses. As T RR increases, stronger weakening of synaptic connections per spike is expected, however, less spikes per second are evoked.

B. Coordinated reset stimulation
During CR stimulation, stimuli are administered periodically. The interval between two subsequent stimuli is denoted by T CR which, together with the finite duration of a single stimulus T , yields the mean stimulation frequency f CR stim = 1/(T CR + T ) and stimulation times s k = k(T CR + T ). CR stimulation is administered to M e = 4 disjoint neuronal subpopulations P i ∈ {P 0 , P 1 , P 2 , P 3 }, resembling a set of four stimulation contacts, a typical setup [69,70] enabling pronounced desynchronizing effects [71]. Each subpopulation contains 25% of the neurons. This is illustrated in Fig. 1(b).
In the CR stimulation setup (Fig. 1) synapses may either connect two neurons that are part of the same subpopulation P l and, hence, are always stimulated simultaneously, or synapses may connect neurons belonging to different subpopulations, which are never stimulated simultaneously. These two types of synapses will be denoted as P CR max and P CR min synapses, respectively. For both types of synapses, we calculated p CR i j (S| ) by considering the time intervals between stimulus deliveries to the presynaptic and the latest and the next stimulus deliveries to the postsynaptic neuron. As stimuli are delivered periodically, these time intervals S are integer multiples of T CR + T . As each subpopulation is stimulated exactly once during a CR cycle, intervals S are bounded from above and below by ±2M e (T CR + T ), respectively.
p CR i j (S| ) for P CR max synapses can be obtained by counting all possible sequences that lead to a certain interval S between post-and presynaptic stimuli. Again assuming t d < T , we The result for P CR min synapses is given in the Appendix, see for the two types of synapses are shown in Fig. 3(g). In contrast to RR stimulation, we find that N CR i j (t ) possesses pronounced peaks at integer multiples of T CR + T λ , reflecting strong temporal correlations in the stimulation sequence.
Results for J ∞ i j are shown in Fig. 3(h) for P CR max (top panel) and P CR min synapses (bottom panel). We find qualitatively similar behavior as for corresponding types of synapses during RR stimulation. However, in contrast to RR stimulation, slow CR stimulation leads to weaker decoupling per spike for P CR min synapses.

VI. CONTROL OF NETWORK TOPOLOGY
The probability that post-and presynaptic neurons are stimulated simultaneously is critical for the expected rate of weight change of a synapse. This probability scales the relative contributions of stimulus-and sequence-induced reshaping to the overall synaptic weight dynamics. If the probability is high, stimulus-induced reshaping strongly contributes to the weight dynamics. Furthermore, stimulus-induced reshaping can either strengthen or weaken synaptic connections, see results for P CR max and P RR max synapses in Figs. 3(c) and 3(h). These effects can be harnessed for controlling network structure by using appropriately designed stimulation protocols. First, by frequent simultaneous stimulation of certain groups of neurons, strong contributions of stimulus-induced reshaping can be realized. And second, adjusting the shape of stimuli allows for either weakening or strengthening corresponding synaptic connections.
We illustrate reshaping of network topology by means of RR and CR stimulation in Fig. 3. Spike trains are generated according to the stimulation-controlled spiking approximation, Eq. (3), using the Gaussian approximation, Eq. (8), for the distribution of spike response times, see the Appendix for details. Snapshots of connectivity matrices, containing the synaptic weights w i→ j , during stimulation with σ λ < t d are depicted in Figs. 3(d) and 3(i), for RR and CR stimulation, respectively. Snapshots for σ λ > t d are shown in Figs. 3(e) and 3(j), respectively.

VII. SHAPE OF STDP FUNCTION DETERMINES STIMULATION-INDUCED WEIGHT DYNAMICS
Synaptic reshaping during stimulation is determined by the STDP function, Eq. (1). In the present paper, we focus on the canonical form of STDP, Eq. (1) [63]. This type of STDP was observed in various brain regions including the neocortex [35] and hippocampus [36]. Synaptic reshaping for other types of STDP is discussed in the Appendix. We evaluate the expected rate of weight change J ∞ i j , Eq. (5), for different STDP parameters for RR and CR stimulation. In particular, we vary the ratio between longterm depression and long-term potentiation β. In Fig. 4, we compare the expected weight change per spike J ∞ i j S i /δ during RR and CR stimulation for P CR max and P RR max synapses, Figs. 4(d)-4(f), and P CR min and P RR min synapses, Figs. 4(g)-4(i), and different shapes of the STDP function W (t ), Eq. (1). We find that decoupling during both RR and CR stimulation is most pronounced for depression-dominated STDP for which all synaptic weights can be reduced for sharp distributions of spike response times λ. Furthermore, for such rules the weakly connected state is typically stable [72]. This indicates that stimulation decouples the neurons, but also leads to longlasting aftereffects.
In general weakening of synapses between subpopulations during CR stimulation shows a pronounced dependence on the interstimulus intervals and therefore on the stimulation frequency, see Figs. 4(h) and 4(i). In contrast, weakening due to RR stimulation is more robust to changes of interstimulus intervals. This results from the qualitatively different shapes of N CR i j (t ) and N RR i j (t ) discussed above. While N CR i j (t ) possesses pronounced peaks at multiples of T CR + T λ , which may or may not be adjusted to time lags at which the STDP function yields effective decoupling, i.e., i j (t ) lacks such peaks and therefore allows for averaging over wide ranges of time lags. The latter typically causes decoupling for depressiondominated STDP rules, see Fig. 4(i).
Additionally, we present similar plots for various different types of STDP functions in the Appendix. In general, we observe that RR stimulation causes robust weakening for P RR min synapses for depression-dominated STDP. Well-adjusted CR, however, my perform better in a limited range of stimulation frequencies. The performance of both stimulation protocols for P CR max and P RR max synapses, respectively, depends on the shape of the STDP function at slightly negative time lags (≈−t d ). If W (−t d ) < 0, this may lead to fast decoupling, exploiting a phenomenon called "decoupling by synchronization" [66,67]. On the other hand, if W (−t d ) > 0, then stimulation strengthens these synapses.

VIII. STIMULATION OF PLASTIC NETWORKS OF INTEGRATE-AND-FIRE NEURONS
We compare theoretical predictions for the limit of stimulation-controlled spiking of the previous sections to computer simulations for RR and CR stimulation of networks of oscillatory excitatory leaky integrate-and-fire neurons with STDP. Integrate-and-fire models describe the subthreshold dynamics of the neurons' membrane potentials. Once the latter pass a spiking threshold spikes are generated. Details on the model are given in Sec. X. In particular, we consider slightly depression-dominated STDP for which a stable desynchronized state with weak synaptic connections coexists with a stable synchronized state with strong synaptic connections, see Fig. 8 and Ref. [72]. To study stimulation-induced decoupling, we prepare the network in the synchronized state. This mimics the initial conditions used in previous computational studies on CR stimulation, see for instance Ref. [73]. Then either RR or CR stimulation is applied for 1 h of biological time. At each stimulation time step, a short charge-balanced stimulus of amplitude A is administered, see Sec. X for details. After cessation of stimulation we simulated the network dynamics for additional 1000 s in order to study long-lasting aftereffects.
To study stimulation-induced reshaping of network connectivity as well as acute and long-lasting desynchronization effects, we determine the time-averaged mean synaptic weight w and the time-averaged Kuramoto order parameter [5] ρ t quantifies the average degree of synchronization in a time interval of length T starting at time t. In particular, ρ t ≈ 1 indicates highly synchronous spiking activity and ρ t ≈ 0 desynchronous activity. In order to calculate ρ t we assign phases φ i (t ) to each neuron i, which increase linearly in time during individual interspike intervals by a total amount of 2π [74]. Results for the Kuramoto order parameter and the mean synaptic weight are shown in Fig. 5. Additionally, raster plots of neuronal spiking activity evoked by either stimulation protocol are shown in Fig. 5(a). When a stimulus of moderate amplitude A is applied to a neuronal subpopulation, only some neurons show spiking responses. This is due to a finite refractory period during which the neurons' membrane potentials are far below the spiking threshold, see Sec. X. Nevertheless we find decoupling during RR and CR stimulation as predicted by our theory, see Fig. 5(b). However, decoupling during RR stimulation is significantly stronger so that the mean 033101-8 synaptic weight approaches zero while its reduction during CR stimulation is weaker. Still, both stimulation protocols drive the network into the attractor of a stable desynchronized state, see Fig. 5(b) (right panel) for asymptotic solutions. In that state, ongoing weight updates due to STDP lead to slightly nonzero mean weight (≈0.01). As weight updates depend on the statistics of time lags, stimulation may lead to even lower mean synaptic weights, see Fig. 5

(b) (left panel).
Following the literature on desynchronization stimulation, we distinguish between acute and long-lasting desynchronization effects [51,52]. Acute desynchronization refers to desynchronization while stimuli are delivered. Long-lasting desynchronization refers to the phenomenon that neuronal activity remains desynchronized after the cessation of stimulation. Results for Kuramoto order parameters ρ AC evaluated during RR and CR stimulation are shown in Figs. 5(c) and 5(e), respectively. Corresponding mean synaptic weights w AC are shown in Figs. 5(d) and 5(f), respectively. We find that the decoupling effect of RR stimulation is more robust with respect to changes of the stimulation frequency than that of CR stimulation, compare Figs. 5(d) and 5(f). Most remarkably, RR stimulation causes decoupling even if only weak acute desynchronization is achieved, compare Figs. 5(c) and 5(d). This is in marked contrast to the decoupling effect of CR stimulation which relies on acute desynchronization, see Figs. 5(e), 5(f), and 5(b). For long-lasting desynchronization, stimulation has to reduce synaptic weights such that the network is driven into the basin of attraction of a stable desynchronized state. In consequence, long-lasting desynchronization effects induced by RR stimulation are more robust with respect to changes of the stimulation frequency than those of CR stimulation, compare Figs. 5(g) and 5(i).
In Fig. 6, we compare simulation results for different stimulus amplitudes A to theoretical predictions, Eq. (5). We find good quantitative agreement for strong stimulation. For weaker stimulation our theory overestimates the decoupling effect of RR and CR stimulation. Remarkably, differences between decoupling rates for the different types of synapses remain even for weak stimulation. Furthermore, we find faster decoupling for RR than for CR stimulation for the same stimulation frequency.

IX. DISCUSSION
Desynchronizing stimulation can decouple oscillatory neuronal networks with STDP [51]. Several computational 033101-9 studies predicted different features of stimulation-induced long-term desynchronization and corresponding symptom relief [73,75,76], which were thereupon verified in animal experiments and clinical studies [53][54][55][56]. However, desynchronization does not necessarily lead to decoupling and may be limited by the need to adjust stimulation parameters to dynamic characteristics of the network, typically its dominant frequency [52]. We here presented a qualitatively different approach to specifically decouple oscillatory neuronal networks: decoupling stimulation. Stimulation is not applied to primarily desynchronize oscillators, but to cause activity patterns that specifically decouple adaptively coupled oscillators. We studied decoupling stimulation in neuronal networks with STDP. We presented a theoretical framework that relates the stimulation-induced synaptic weight dynamics to the response characteristics of individual neurons to individual stimuli as well as spatiotemporally correlated sequences of stimuli. Based on our findings, we presented an effective decoupling stimulation technique-RR stimulation-that overcomes the need to adjust parameters to the dominant frequency.
The presented theoretical framework holds for resetting stimuli where each spike results from the stochastic response of individual neurons to an administered stimulus. It provides estimates for the rate of synaptic weight changes during stimulation of neuronal networks with STDP. Our theory can be applied to stationary and nonstationary systems and does not rely on a particular implementation of STDP. A critical assumption is that the stimulation allows for controlling neuronal spike timings. Corresponding high response fidelity of neuronal spikes to administered stimuli has been observed during deep brain stimulation [77,78], an established therapy for Parkinson's patients [79]. Beyond that, our theoretical framework is valid for a large number of systems because phase resetting by delivery of a stimulus is a universal phenomenon. Phase resetting has been demonstrated in a large number of theoretical and experimental studies on hyperpolarizing or depolarizing electrical pulses [80][81][82][83] as well as excitatory or inhibitory postsynaptic potentials [84][85][86][87][88]. Accordingly, stimulusinduced neuronal phase resetting was observed in a variety of different systems, ranging from simple squid axons [89] to central pattern generators for respiratory rhythm [90][91][92], neuronal population activity in brain slices of inferior olive [93] and consummatory licking [94]. Phase resetting of cortical rhythms has been realized using several stimulation techniques, e.g., transcranial magnetic stimulation [95] and sensory stimulation, e.g., visual and auditory stimulation [96][97][98].
Our theoretical framework can be applied to a general class of stimulation protocols and can be used for the design of stimulation techniques which specifically weaken or strengthen plastic synaptic connections. It relates the stimulation-induced weight dynamics to the response characteristics of neurons to a single stimulus and to spatiotemporal correlations between stimuli. Because of its implications to neuroscience and medicine [51,[54][55][56], we here focused on decoupling stimulation. Instead of primarily inducing desynchronization with stimuli tuned to the network's dominant frequency, we directly decouple the neurons by delivering RR stimulation, a stimulation pattern with appropriate spatiotemporal randomization. RR stimulation differs significantly from previous approaches using randomized stimuli, e.g., delivering common white Gaussian noise [29] or uncorrelated noise input [50]. While common noise can be adjusted to cause either acute synchronization or desynchronization [29], uncorrelated noise has been found to strengthen synaptic connections and thereby support synchrony in plastic neuronal networks [50]. Spatiotemporal patterns of phaseresetting stimuli delivered during RR stimulation, on the other hand, synchronize simultaneously stimulated neurons, while subpopulations receiving subsequent stimuli desynchronize temporarily. The degree of resulting partial acute synchrony depends on the number of simultaneously stimulated neurons. The two limiting cases of simultaneous stimulation of all neurons and single neuron stimulation result in perfect acute synchronization and Poisson like spike trains, respectively, in the theoretical limit of strong and fast stimulation.
While our theory does not require a particular type of STDP, the latter has a great impact on the emergence of dynamical states in plastic neuronal networks. It has been found that the interplay of synaptic and spiking dynamics can lead to a rich repertoire of emerging network motifs [72,99] and coexisting dynamical states in networks of oscillatory neurons. These states include desynchronized, synchronized and cluster states, see, for instance, Refs. [41,46,73,100]. Previous theoretical studies managed to connect the shape of the STDP function to the emerging dynamical states. It has been found that antisymmetric STDP [τ R = 1 in Eq. (1)], with balanced depression and potentiation, (β = 1), leads to the emergence of a single synchronous state in networks of oscillatory neurons [39]. A general theory for the dynamics of synaptic weights in weakly connected excitatory recurrent neuronal networks was presented in Ref. [72]. Analyzing the shape of W (t ), the authors reported that a stable state with strong synaptic connections emerge for strongly potentiation-dominated STDP (β 1) and a stable state with weak synaptic connections for strongly depression-dominated 033101-10 STDP (β 1). For values of β that are close to one, different states may coexist. In particular, for STDP functions with asymmetric STDP decay times (τ R > 1) and β/τ R < 1 that are slightly depression-dominated (β 1) a complex motif structure and coexisting states with strong and weak synaptic connections, respectively, can emerge [72]. States with strong and weak synaptic weights, are typically associated with synchronized and desynchronized spiking activity, respectively.
Based on our theoretical framework, RR stimulation was designed to induce decoupling in networks with depressiondominated STDP, see Figs. 4, 9, 10, and 11, by inducing a corresponding statistics of time lags between spikes of postand presynaptic neurons without primarily desynchronizing them. Such STDP favors the existence of a stable state with weak synaptic weights [72] and desynchronized activity. We therefore hypothesize that sufficiently long RR stimulation of such networks induces long-lasting desynchronization that outlasts stimulation. Furthermore, our theory predicts that decoupling effects of RR stimulation are more robust with respect to changes of the stimulation frequency than effects of CR stimulation for depression-dominated STDP. However, CR stimulation may perform better in a limited range of stimulation frequencies, Fig. 4(i).
We compared our theoretical predictions for decoupling effects caused by RR and CR stimulation to simulations of networks of oscillatory excitatory leaky integrate-and-fire models with STDP. In these networks, stable desynchronized and synchronized states coexist for slightly depressiondominated STDP, see Fig. 8. We find good quantitative agreement between theory, which assumes stimulation-controlled spiking, and simulations for sufficiently high stimulus amplitudes enabling stochastic phase resetting of individual stimuli. For moderate stimulus amplitudes decoupling is slower than predicted as nonstimulus locked network dynamics weakens control of the collective dynamics by the stimuli delivered, see Fig. 5(a). However, the theoretically predicted qualitative differences between decoupling rates of synapses within and between separately stimulated populations are still valid for moderate stimulus amplitudes, see Fig. 6. Unlike CR, RR stimulation has pronounced decoupling effects during and after stimulation [Figs. 5(d) and 5(h) and then 5(f) and 5(j)]. The decoupling effects of RR stimulation are robust with respect to variations of the frequency ratio ν = r synch / f stim and do not require the stimulation frequency to be adapted to the mean frequency of the synchronized rhythm r synch . This translates to pronounced long-lasting desynchronization, achieved over a wide range of frequency ratios ν [ Fig. 5(g)].
Remarkably, pronounced desynchronization during stimulation is not required for significant changes of the synaptic connectivity and corresponding long-term desynchronizing outcome. In fact, during RR stimulation synchronization is only moderately reduced [ Fig. 5(c)]. In contrast, by design CR stimulation causes pronounced desynchronization within a considerable range of frequency ratios ν [ Fig. 5(e)]. However, desynchronization during CR stimulation does not necessarily imply reduction of synaptic weights (see the green vertical stripes in [Fig. 5(f)]. Accordingly, the parameter range associated with long-term decoupling [ Fig. 5(j)] and long-term desynchronization [ Fig. 5(i)] is smaller and intersected by vertical stripes in the parameter plane spanned by stimulation amplitude A and frequency ratio ν. The latter occur for values slightly exceeding small integer ratios ν. For those frequency ratios a considerable amount of stimuli gets delivered to neurons during their refractory period which reduces efficiency of CR stimulation.
To enhance long-lasting desynchronization effects, a series of computational studies focused on modifying the sequence by which stimuli are delivery to different stimulation sites. For this, rapidly or slowly varying CR sequences [70,76], periodic flashing sequences [101] as well as various stochastic versions thereof [102] were used. Recently, it was also suggested to de-033101-11 liver neuron-type specific stimuli in order to induce a constant time shift between spikes of excitatory and inhibitory neurons as STDP would lead to a modulation of the synaptic coupling between both types of neurons and thereby lead to transitions between desynchronized and synchronized states [103]. Our results reveal two basic mechanisms for stimulationinduced reshaping of network connectivity: stimulus-induced and sequence-induced reshaping. The former characterizes synaptic weight changes in response to individual stimuli and the latter refers to weight changes caused by spatiotemporal correlations in the sequence of stimuli. By appropriately designing a stimulation protocol, the two mechanisms can be adequately combined to weaken or strengthen synaptic connections within and between individual neuronal subpopulations. Accordingly, the theoretical framework presented here allows for a quantitative understanding of the results from computational studies mentioned above and can be used to optimize the outcome of these multisite stimulation protocols or design novel stimulation techniques. Spatial randomization, as put forward in this study, can technically be achieved by means of segmented stimulation contacts [104] combined with superposition-based spatial steering of stimulation current [105].
Both stimulus-and sequence-induced reshaping can cause a strengthening or weakening of synaptic connections. Their respective impact on network structure depends on the STDP function. In the main text, we focused on the canonical Hebbian form of STDP as presented in Refs. [35,36]. However, various different forms of STDP were observed in experiments [106]. We discuss both types of reshaping for other Hebbian and anti-Hebbian STDP functions in the Appendix. Stimulus-induced reshaping results from time lags between spiking responses triggered by a single stimulus. For phaseresetting stimuli these time lags are short compared to the STDP decay times, which are typically in the range of 10-100 ms [106]. The resulting weight change is determined by the interplay of STDP, synaptic transmission delays, and the shape of the distributions of spike response times. Decoupling occurs if the distributions of spike response times are sharp compared to the effective synaptic transmission delays. In that case decoupling takes place for any type of STDP where longterm synaptic depression occurs when postsynaptic spikes arrive at the synapse shortly before presynaptic ones. This includes the canonical Hebbian STDP with balanced and unbalanced depression and potentiation as, for instance, observed in hippocampal cultures [36] and pyramidal neurons in rat cortical slices [107], as well as anti-Hebbian STDP, as, for instance, observed for excitatory synapses on inhibitory interneurons in neocortex [108]. Corresponding decoupling was reported by studies on highly synchronized population bursts [66,67] and dominates synaptic reshaping of synapses between neurons that are likely to receive stimuli simultaneously during RR and CR stimulation, see Figs. 4(d)-4(f) and Fig. 11 in the Appendix for respective STDP functions. A strengthening of synaptic connections due to stimulus-induced reshaping may occur for STDP functions that cause potentiation for small negative time lags, for instance symmetric Hebbian STDP, see Figs 6) and (5) using p CR i j and p RR i j from Sec. V. Results for P RR max and P CR max synapses (orange) are shown in panels (d)-(f) and results for P RR min and P CR min synapses (blue) are shown in panels (g)-(i). λ(t ) is approximated by a Gaussian distribution with σ λ = 1 ms. Parameters: τ + = 10 ms, τ R = 4. tic spikes that arrive at the synapse shortly before postsynaptic ones lead to strong potentiation. The latter holds for the canonical Hebbian STDP function [36]. On the contrary, sequenceinduced reshaping results from spatiotemporal correlations between stimuli. The latter determine the temporal structure of stimuli received by post-and presynaptic neurons. Resulting time lags between the neurons' spikes are of the order of interstimulus intervals which are of the same order as STDP decay times for stimulation frequencies, e.g., commonly used in CR stimulation [54][55][56]. Therefore, the resulting effect on the network structure is weaker [Figs. 6(c) and 6(d)] and depends on the shape of the STDP function for time lags of the order of the mean interstimulus intervals S i . For stimulation protocols with strong temporal correlations between stimuli, e.g., CR stimulation [69] or periodic flashing CR stimulation [101], our theoretical framework predicts decoupling of stimulated populations for STDP functions with W (− S i − t d ) + W ( S i − t d ) < 0, see panels (h) and (i) in Figs. 4, 9, 10, and 11. The connectivity matrices presented in Ref. [69] indeed confirm decoupling of separately stimulated populations. On the other hand, randomized stimulation protocols, e.g., RR stimulation, lack strong temporal correlations between stimuli. Then, the expected synaptic weight change is similar to the one resulting for independent Poisson spiking of pre-and postsynaptic neurons. The latter typically leads to a weakening of synapses in networks with depressiondominated STDP and stable desynchronized states, see panels (h) and (i) in Figs. 4, 9, 10, and 11. During RR stimulation, synapses weaken when it is unlikely that post-and presynaptic neurons receive stimuli simultaneously, otherwise, stimulusinduced reshaping dominates the synaptic dynamics. The ratio of respective synapses is controlled by the number of simultaneously stimulated neurons. The latter can be varied in order to scale contributions of stimulus-and sequenceinduced reshaping to the overall synaptic weight dynamics. We therefore expect RR stimulation to be suitable to drive such a network into the attractor of a stable desynchronized state.
Here we presented results for nearest neighbor STDP schemes, see Figs. 1(a) and 1(c). Nearest neighbor STDP schemes were used in previous studies on desynchronizing stimulation in plastic neuronal networks, see for instance [52,73]. Other studies considered all-to-all schemes in which all possible pairings between post-and presynaptic spikes are considered [62,68,109]. A detailed discussion on different STDP schemes is presented in Ref. [65]. In fact, experimental studies on cortical neurons reported that the depencence of synaptic weight changes on the firing rate is better described by nearest neighbor than by all-to-all schemes [110]. Our theoretical result, Eq. (5), can be applied to other STDP schemes by considering corresponding spike pairs for the calculation of N i j (t ) and p i j (S| ) in Eq. (6). For all-to-all STDP schemes, N i j (t ) possesses higher counts for negative time lags, because delayed spike arrivals from simultaneous stimulation of pre-and postsynaptic neurons are paired with all prior postsynaptic spikes. We may therefore expect an even stronger decoupling effect for RR stimulation for allto-all STDP schemes with pronounced depression occurring when presynaptic spikes appear shortly after postsynaptic ones. For fast-decaying STDP functions, where STDP decay times τ + and τ R τ + are short compared to typical time lags t post − t pre , all-to-all and nearest neighbor schemes result in similar weight dynamics [111].
Our implementation of the RR protocol assumes that separate stimulation of post-and presynaptic neurons is possible. This assumption may not hold in an experimental setup due to limited spatial resolution. In an experiment, segmented multisite electrodes with several stimulation contacts [60] could be used to implement stimulation patterns similar to our RR protocol. Further stimulation techniques allowing for multisite stimulation include noninvasive vibrotactile stimulation [112] and optical stimulation [113]. Using four stimulation sites, as for instance used in the standard CR setup, each stimulus could be delivered to two out of four randomly selected stimulation sites. Our theory can be applied to this setup in order to predict the synaptic weight dynamics in the limit of stimulation-controlled spiking. The probability to receive stimuli simultaneously scales the contributions of stimulusand sequence-induced reshaping to the overall weight dynamics. For synapses between neurons located in the vicinity of the same stimulation site, this yields very similar decoupling properties as for the implementation of RR stimulation considered here in the case of neurons with nearby indices, see panels (d)-(f) of Figs. 4, 9, 10, and 11, for respective STDP functions. On the other hand, if post-and presynaptic neurons are in the vicinity of different stimulation sites, both Eqs. (11) and (13) contribute to the weight dynamics. Contributions are scaled by the probability that both sites are activated simultaneously. Decoupling properties can be constructed by superposition of the curves shown in panels (d)-(f) and (g)-(i) for RR stimulation in Figs. 4, 9, 10, and 11. Our results provide a theoretical base for the prediction of synaptic reshaping due to stationary or nonstationary spatiotemporal stimulus patterns. Stationary patterns include periodic high-frequency stimulation or coordinated reset stimulation, which are delivered via deep brain stimulation as treatment for neurological disorders such as medically refractory Parkinson's disease [56,114]. Nonstationary patterns include, for instance a stepwise adjustment of stimulation parameters, which might be advantageous in heterogeneous networks where reshaping of certain synapses first might strongly affect the network's response to later stimuli. In future studies, we anticipate to study synaptic reshaping in detailed models that incorporate disease-specific network connectivity. For instance, in the basal ganglia region, which includes major target regions for deep brain stimulation in Parkinson's patients, e.g., the STN, different segregated portions of neurons are responsible for processing different sensory information 033101-13 [115,116]. While to date the exact mechanism leading to motor symptoms is not completely understood, a major hypothesis is a loss of functional segregation due to unwanted up-regulated connectivity between microcircuits [115]. Experimental studies reported stimulation-induced reshaping of excitatory connectivity which may contribute to therapeutic effects [117][118][119]. In our generic model, RR stimulation led to a complete decoupling of neurons. It is an open question whether RR stimulation is also able to restore physiological network connectivity. A previous computational study provided evidence that CR stimulation of the STN specifically down-regulates pathological connectivity and restores physiological synaptic connectivity patterns related to spatially patterned time-correlated input [120]. In this context, it is key to understand how multisite stimulation can be adjusted to down-regulate certain synaptic connections, while not affecting or even up-regulating others. Our results may provide critical insight in this regard and may enable the design of corresponding stimulation protocols. In general, our approach can be applied to various networks with adaptive coupling in different fields of science. Apart from decoupling neural networks with pathology-related strong synaptic couplings (see above), our approach may also be used to probe neuronal synchrony-mediated physiological information processing, e.g., in the context of sensory binding [121,122]. However, our decoupling approach may also be applied to systems beyond neuroscience and neurology, e.g., networks with adaptive connections used to study opinion formation and epidemic spreading [123][124][125]. Spatiotemporal stimulus patterns that trigger a decoupling of nodes or groups of nodes in such networks could present an effective tool to control network dynamics.

A. Leaky integrate-and-fire model
We perform simulations for excitatory networks of N = 1000 conductance-based oscillatory leaky integrate-and-fire neurons. The range of membrane potential oscillations and the single neurons' firing rates were adjusted to match recordings of oscillatory neurons in the subthalamic nucleus, a major target brain region for deep brain stimulation in Parkinson's disease patients [126,127].

Subthreshold dynamics
The dynamics of the ith model neuron's subthreshold membrane potential V i (t ) is given by C i is the membrane capacity, g leak the leakage conductance, and V rest the resting potential. The second term is the excitatory input current with reversal potential V syn , modeled by a time-dependent conductance g syn,i (t ), see below. I stim,i (t ) is the stimulation current and I noise,i (t ) is a noisy input current accounting for Poisson input from other brain regions, see below.

Spike generation
Whenever the membrane potential crosses the dynamic threshold V th,i (t ) a spike is generated and the dynamic threshold value is set to V th,spike . The threshold dynamics during individual interspike intervals is given by We implement a simple rectangular spike shape by setting the membrane potential to V spike for a duration of τ spike . Afterward, V i is set to V reset .

Synaptic dynamics
We consider excitatory coupling between neurons. The synaptic conductances g syn,i (t ) obey Here τ syn is the synaptic timescale, N the number of neurons, t j l j is the timing of the l j th spike of neuron j, and t d the synaptic delay. We restrict ourselves to axonal delays as these are typically longer than dendritic delays. The outer sum runs over all presynaptic neurons G i of neuron i. The inner sum runs over all spike times of the presynaptic neuron j. κ is the maximal coupling strength and w j→i ∈ [0, 1] is the synaptic weight of the synapse between neurons j and i.

Poisson input
In addition to presynaptic input, each neuron receives noisy excitatory input I noise,i (t ). To this end, we generate independent Poisson spike trains with constant firing rates f noise . These are fed into excitatory synapses as presynaptic spike trains such that the noisy input current to neuron i is given by With synaptic conductance g noise,i obeying where κ noise scales the noise intensity and t i k i is the time of the k i th spike of the Poisson spike train fed into neuron i. A time trace of the membrane potential for a single neuron is shown in Fig. 7.
In order to ensure heterogeneity of the single neurons' firing rates, membrane capacitances were Gaussian distributed with mean C i and standard deviation 0.05 C i . All model parameters are given in Table I.

Network topology
For reasons of comparability, the considered network structure mimics that of earlier studies on coordinated reset stimulation of STN neurons Ref. [128]. First, neurons were distributed in an elipsoidal volume with semi-principle axes, a x = 2.5 l scale , a y = 6.0 l scale , a z = 3.0 l scale . These were 033101-14 chosen such that the resulting volume for l scale = 1 mm fits experimental measurements on the subthalamic nucleus presented in Ref. [129]. For a given total number of neurons N, the scaling parameter l scale is used to artificially shrink the volume in a way that the mean distance between connected neurons matches experimental data for rat STN neurons, r d ≈ 0.543 mm [130]. Results presented in the main text are for N = 1000. We performed additional simulations for N = 400, 800, and 1600 which led to qualitatively similar results. N neurons were uniformly distributed in that volume. Synaptic connections were introduced such that the total connectivity amounts to about 7% compared to all-to-all coupling [128], which yields on average 0.07N connections per neuron. The probability for two neurons to connect is distance dependent and given by p(d ) ∝ exp (−d/d c ), with d c = 0.5 l scale [128]. The resulting networks have a broader degree distribution than networks generated with p independent of the distance between neurons. We performed additional simulations for networks with p independent of the distance between neurons and found qualitatively similar result as for the ones presented in the main text. All networks used throughout the manuscript were fully connected. For all simulations that were presented in the main text we set l scale = 0.35 mm, which sets the mean length of synaptic connections to ≈0.545 mm.

Simulation details
Initially, the neurons' membrane potentials are distributed according to a uniform distribution V i ∈ [V reset , V rest ], thresholds V th,i are set to V th,rest , and synaptic conductances g syn,i and g noise,i are set to zero. Synaptic weights are distributed according to a bimodal distribution with w i→ j (t = 0) ∈ {0, 1} such that a given mean weight w initial := w i→ j (t = 0) is obtained. In order to wait for the neuronal activity to equilibrate STDP is switched off for the first 20 s. Simulations were performed using the explicit Euler method with time step h = 0.1 ms. Stationary networks were obtained by simulating until the mean synaptic weight equilibrates.

B. Coexistence of synchronized and desynchronized state
We performed a detailed analysis of the network's multistability. For the parameters given in Table I, we found coexisting states of synchronized and desynchronized activity for depression-dominated STDP functions, Eq. (1), with asymmetric decay times τ R > 1. This is in accordance with theoretical predictions for weakly connected networks presented in Ref. [72]. In the main text we use τ R = 4, see Fig. 8(d).
Raster plots of neuronal spiking activity in desynchronized and synchronized states are shown in Figs. 8(a) and 8(b), respectively. Accordingly, corresponding distributions of synaptic weights and individual neurons' firing rates are depicted in Figs. 8(c) and 8(e), respectively. In order to distinguish between desynchronized and synchronized activity, we calculate the time-averaged Kuramoto order parameter from Eq. (15) as a function of the initial mean synaptic weight for stationary networks. Results are shown in Fig. 8(f). We find coexisting desynchronized and synchronized states for depressiondominated STDP (where β > 1). For the parameter set used in the main text, i.e., τ R = 4 and β = 1.4, networks with w initial 0.3 approach a stationary state with synchronized neuronal activity and mean synaptic weight w i→ j ≈ 0.38. Networks with w initial 0.25 approach a stationary state of desynchronized activity. We find that network realizations with mean initial weights 0.25 w initial 0.3 can approach either state.
To study the effect of stimulation on the network dynamics, we use the final network states of our study on multistability as initial networks. For simulation results presented in Figs. 5 and 6, we used the networks with initial mean synaptic weight w initial = 0.5. These networks were simulated for 2000 s of biological time in order to wait for the mean synaptic weight to become stationary. Longer simulations showed that their mean weights fluctuated around an average value of w i→ j ≈ 0.38.

C. Stimulus waveform
During stimulation, sequences of charge-balanced pulses were delivered to neurons in respective subpopulations. This is modeled by applying the stimulation current I stim,i (t ) = AX (t ).
A is the stimulation amplitude in units mS/cm 2 . Following earlier studies [69], X (t ) in units of mV is given by the sum of stimuli administered to neuron i. Each stimulus is characterizes by a stimulus wave form consisting of two rectangular pulses. The first excitatory one of amplitude 1 mV and duration of 0.4 ms is followed by an inhibitory one of amplitude −4/30 mV duration 3 ms. Both, rectangular pulses are separated by a gap of 0.2 ms. We set the minimum time between subsequent stimulus deliveries to T λ ≈ 7.69 ms. This corresponds to a maximal stimulation frequency of 130 Hz. The stimulus wave form and a representative response of a neuron is shown in Fig. 7.

D. Generation of spike trains according to stimulation-controlled spiking approximation
Spike trains according to the stimulation-controlled spiking approximation, Eq. (3), were generated as follows: First, a stimulation sequence S is generated according to one of the stimulation protocols, i.e., RR or CR stimulation. Then, at the stimulation times s k one spike time t i is generated for each stimulated neuron, i ∈ P k . These times are distributed according to λ(t i − s k ). In Figs. 3(a) and 3(b), we set λ(t i − s k ) ∝ exp((t i − s k ) 2 /2σ 2 λ ). To study stimulation-induced network topology [Figs. 3(d), 3(e), 3(i), and 3(j)], we generated all-to-all coupled networks with 1000 nodes. The initial mean weight is set to w initial = 0.38, which is realized by distributing initial weights according to a bimodal distribution w i→ j ∈ {0, 1} with this mean value. Then, spike trains for each node were generated according to the stimulation-controlled spiking approximation. Given these spike trains, we calculated the weight dynamics of individual synaptic weights w i→ j by applying the nearest neighbor STDP scheme from Eq. (1). The connectivity matrices presented in Figs. 3(d), 3(e), 3(i), and 3(j) show snapshots of individual synaptic weights w i→ j .

ACKNOWLEDGMENTS
We gratefully acknowledge support of this study by Boston Scientific Neuromodulation (Grant No. IC2019-0221). Some of the computing for this project was performed on Stanford's Sherlock High-Performance 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.

APPENDIX A: RESHAPING OF NETWORK CONNECTIVITY DURING CR STIMULATION
We calculate p CR i j (S| ) for P CR min synapses. Without loss of generality, we assume that the presynaptic neuron i is part of subpopulation P q i and the postsynaptic neuron j is part of subpopulation P q j . The mean number of occurrences for interstimulus intervals between subsequent stimulations of these neurons can be obtained as follows. We consider a stimulus delivered to the presynaptic neuron at time s k . It triggers a spike of the presynaptic neuron. As the postsynaptic neuron does not receive a stimulus at time s k , time lags involving this presynaptic spike can only result from pairings with postsynaptic spikes triggered by earlier stimuli, i.e., at time s k−ξ , or by future stimulation of the postsynaptic neuron, i.e., at time s k+ξ . Here, ξ is a natural number. For evaluation of time lags resulting for the former case, we consider all stimulation sequences with a stimulus delivered to P q i at time s k and a stimulation of P q j at time s k−ξ but not at times s k−ξ +1 , s k−ξ +2 , . . . , s k−1 . For the latter case, we consider all stimulation sequences with a stimulus delivered to P q i at time s k and stimulation of P q j at time s k+ξ but not at times s k+ξ −1 , s k+ξ −2 , . . . , s k+1 . The mean number of occurrences per presynaptic spike are listed in Table II for all possible ξ .  for P q i = P q j . (A1) Using this in Eq. (6) yields the result for N CR i j (t ) which is depicted in Fig. 3(g).

APPENDIX B: PREDICTED PERFORMANCE OF CR AND RR STIMULATION FOR DIFFERENT STDP FUNCTIONS
We calculate J ∞ i j using Eqs. (5) and (6) for different shapes of STDP functions. In the brain, different types of Hebbian and anti-Hebbian STDP have been observed, see for instance Ref. [106]. For both types, qualitatively different profiles have been reported [106]. Here, we consider symmetric and asymmetric shapes. An asymmetric Hebbian STDP function was considered in the main text. In the following, we consider an symmetric Hebbian one, with weight updates given by a symmetric anti-Hebbian one and an asymmetric anti-Hebbian one Parameters are chosen as in the main text. β scales the ratio between long-term depression and long-term potentiation such that β > 1 yields depression-dominated, β = 1 balanced, and β < 1 potentiation-dominated STDP.
Results for the symmetric Hebbian STDP function W SH (t ), Eq. (B1), are shown in Fig. 9. We find fast strengthening of P RR max and P CR max synapses, Figs. 9(d)-9(f), and weakening of P RR min and P CR min synapses for balanced and depressiondominated STDP functions. For asymmetric anti-Hebbian STDP, Fig. 10, we find mainly strengthening of P RR max and P CR max synapses, Figs. 10(d)-10(f), and weakening of P RR min synapses for balanced and depression-dominated STDP functions. Weakening P CR min synapses by CR stimulation is only observed for a very narrow range of small interpulse intervals, Fig. 10(g)-10(i). For symmetric anti-Hebbian STDP functions, Fig. 11, fast weakening occurs for P RR max and P CR max synapses during both types of stimulation. Here, we only find weakening of P RR min and P CR min synapses for depressiondominated STDP. For the latter, weakening of P RR min synapses during RR stimulation is more pronounced that weakening of P CR min synapses during CR stimulation. For all types of depression-dominated STDP functions, we find that weakening of P RR min synapses during RR stimulation is more robust with respect to changes of interstimulus intervals, than weakening of P CR min synapses during CR stimulation.