Excessive neuronal synchrony is a hallmark of several neurological disorders, e.g., Parkinson’s disease. An established treatment for medically refractory Parkinson’s disease is high-frequency deep brain stimulation. However, it provides only acute relief, and symptoms return shortly after cessation of stimulation. A theory-based approach called coordinated reset (CR) has shown great promise in achieving long-lasting effects. During CR stimulation, phase-shifted stimuli are delivered to multiple stimulation sites to counteract neuronal synchrony. Computational studies in plastic neuronal networks reported that synaptic weights reduce during stimulation, which may cause sustained structural changes leading to stabilized desynchronized activity even after stimulation ceases. Corresponding long-lasting effects were found in recent preclinical and clinical studies. We study long-lasting desynchronization by CR stimulation in excitatory recurrent neuronal networks of integrate-and-fire neurons with spike-timing-dependent plasticity (STDP). We focus on the impact of the stimulation frequency and the number of stimulation sites on long-lasting effects. We compare theoretical predictions to simulations of plastic neuronal networks. Our results are important regarding CR calibration for two reasons. We reveal that long-lasting effects become most pronounced when stimulation parameters are adjusted to the characteristics of STDP—rather than to neuronal frequency characteristics. This is in contrast to previous studies where the CR frequency was adjusted to the dominant neuronal rhythm. In addition, we reveal a nonlinear dependence of long-lasting effects on the number of stimulation sites and the CR frequency. Intriguingly, optimal long-lasting desynchronization does not require larger numbers of stimulation sites.

Abnormally strong neuronal synchrony is associated with several neurological disorders, e.g., essential tremor1 and Parkinson’s disease.2 An established treatment is high-frequency deep brain stimulation (HF DBS); however, symptoms typically return shortly after stimulation ceases.3 Modern, theory-based approaches, such as coordinated reset (CR) stimulation,4 specifically counteract abnormal synchrony. During CR, multiple stimulation sites are used to deliver stimuli to different neuronal subpopulations. Previous studies analyzed CR stimulation in neuronal networks with spike-timing-dependent plasticity (STDP). In such networks, synchronized states with strong synaptic connections may coexist with desynchronized states with weak connections. CR stimulation can weaken synapses and drive the network into the basin of attraction of a stable desynchronized state. This may explain long-lasting symptom relief observed in preclinical and clinical studies.5–7 We study how parameters such as the number of stimulation sites and the CR frequency impact long-lasting effects. We find that long-lasting desynchronization is most pronounced when CR parameters are adjusted to the STDP and not to the dominant neuronal rhythm, as suggested by prior studies.4,6,8 Furthermore, for certain numbers of stimulation sites, CR has pronounced long-lasting effects, while it stabilizes abnormal synchrony for others. Our results indicate that tuning the spatial resolution of lead electrodes may help to exploit neuronal plasticity for long-lasting therapeutic effects.

Synchronization phenomena are observed in various fields of the natural sciences.9–17 Often, synchronization is critical for functional performance. In some systems, however, excessive synchrony impairs functionality. In the brain, for instance, abnormally strong neuronal synchrony has been related to several neurological disorders, including essential tremor,1 Parkinson’s disease,2 and epilepsy.18 

An established treatment for advanced Parkinson’s disease is high-frequency (>100Hz) deep brain stimulation (HF DBS). Sufficiently strong HF DBS of the subthalamic nucleus (STN) is able to suppress Parkinson’s symptoms such as dyskinesia, tremor, rigidity, and bradykinesia.19,20 However, its mechanism of action remains enigmatic, and symptoms return shortly after cessation of stimulation;3 thus, the clinical standard is permanently delivered DBS. This increases the risk of unwanted side effects such as dyskinesia, depression, cognitive decline, speech difficulty, instability, and gait disorders.20 Moreover, the limited longevity of batteries in internal pulse generators of DBS devices requires the patient to undergo battery-replacement surgeries.21 

Early experimental and theoretical desynchronization studies in cardiology and neuroscience focused on the complex response characteristics of synchronized oscillators to single stimulus pulses.22–26 While weak stimuli slightly advance or delay the phase of the collective rhythm,24 a single stimulus pulse of moderate amplitude delivered at a vulnerable phase may temporarily desynchronize networks of coupled oscillators.16 Corresponding desynchronization by repetitive or demand-controlled pulse delivery has been demonstrated computationally16,27 and in experiments on coupled electrochemical oscillators.28 In contrast, a sufficiently strong stimulus can cause a phase reset, where phases of simultaneously stimulated oscillators align regardless of their individual phases before stimulus administration.16,23,24 Thus, phase-resetting stimuli provide control over the collective dynamics of oscillators.

Other stimulation approaches intend to desynchronize oscillators by ongoing pulse stimulation. A weak periodic forcing with suitable stimulation frequency leads to phase synchronization;29 nonetheless, tuning the stimulation frequency out of these parameter regions results in chaotic desynchronization.30 The latter might contribute to the therapeutic effects of deep brain stimulation.30 One can calculate the desynchronizing stimulation frequencies from estimated phase response curves of the synchronous rhythm.30 Based on phase response curves of individual oscillators, stimuli can be timed such that they cause growing phase differences. A corresponding closed-loop setup for invasive deep brain stimulation was suggested in Ref. 31. Other stimulation techniques deliver linear or nonlinear delayed feedback either permanently or on-demand.32–38 

Limitations of demand-controlled single pulse stimulation and complex stimuli combining a strong phase reset followed by a desynchronizing single pulse28,39 led to the development of CR stimulation, a multisite stimulation technique that does not require precisely timed stimulus delivery (e.g., targeting a vulnerable phase of the synchronized oscillation).4 As predicted by computational studies in neuronal networks with STDP,40,41 acute effects of CR stimulation may also entail long-lasting desynchronization and therapeutic after-effects.5–7Acute effects denote stimulation-induced effects observed during stimulus administration, whereas long-lasting aftereffects refer to sustained effects emerging after cessation of stimulus delivery for t.

In general, a desynchronizing stimulation may cause acute desynchronization. In plastic neuronal networks, this may entail a long-lasting effect if stimulation causes adequate plastic changes in network structure. However, acute desynchronization does not necessarily result in long-lasting desynchronization.8 Below, we briefly review current research on acute and long-lasting desynchronization by CR stimulation and present novel results on the impact of the number of stimulation sites on long-lasting effects.

Originally, CR stimulation was introduced in a network of N all-to-all coupled Kuramoto oscillators.4 A phase-resetting stimulus delivered to oscillator i at phase ψi(t) resets its phase to a constant value. CR stimulation is typically delivered periodically at CR frequency fCR. Each oscillator is stimulated exactly once during a CR cycle of length 1/fCR. For illustration, we consider a set of synchronized phase oscillators with collective frequency ω, where oscillators receive phase-resetting stimuli at times ti=2π(i/N+kcycle)/fCR (kcycle counts CR cycles). The resulting phase distribution can be tuned by varying fCR. Acute desynchronization occurs for a uniform distribution of phases. It is typically obtained when the CR frequency matches that of the dominant synchronous rhythm fCRω.

In the Kuramoto model, the order parameter (reflecting in-phase synchronization, see below) acts on a slow time scale and interacts with higher-order modes (reflecting different types of cluster states) that operate on a fast time scale.16 

The initial development of CR stimulation aimed at explicitly using the center manifold theorem42 (a dynamical systems theorem related to the slaving principle43 of statistical physics) by delivering CR stimulation to a set of M neuronal subpopulations using, e.g., multipolar deep brain electrodes.6 The underlying goal was to quench the order parameter (i.e., synchrony) by temporarily exciting higher-order modes, which, in turn, get rapidly quenched by the vanishing order parameter based on the center manifold theorem (for review, see Refs. 4 and 44). Considering the setup of synchronized oscillators, we group the N oscillators into M subpopulations; and apply stimuli at times ti=2π(ki/M+kcycle)/fCR. Here, ki=0,1,2,3,,M1 refers to the subpopulation of oscillator i. Sufficiently fast stimulation leads to an M-cluster state with synchronized activity within each cluster. The degree of synchrony can be quantified using the Kuramoto order parameter11 

ρ¯Δ(t)=1ΔtΔ2t+Δ2dt|1Nk=0N1eIψj(t)|.
(1)

ρ¯Δ quantifies the degree of synchronization during the time interval Δ. ρ¯Δ1 refers to maximum in-phase synchronization and ρ¯Δ0 to the absence of in-phase synchronization. Upon cessation of stimulation, the cluster state becomes unstable and breaks. Then, the system slowly resynchronizes via a transient. Figure 1 illustrates this concept for M=4 and the neuronal model used throughout the paper. Stimulation is turned on and a four-cluster state appears. As stimulation ceases, the cluster state becomes unstable and the system resynchronizes via a long transient.

FIG. 1.

Acute desynchronization by CR stimulation. (a) Raster plot of spiking activity for the integrate-and-fire model without STDP (η=0) (see Sec. V). Stimuli (yellow) are applied to M=4 neuronal subpopulations with fCR=3.8Hz and stimulation amplitude Astim=0.1. y-axis shows neuron indices sorted according to their coordinates. (b) Time-averaged Kuramoto order parameter [Eq. (1)] before, during (yellow), and after cessation of stimulation (Δ=0.5s). Markers indicate evaluation times t of Eq. (1). Stimulation induces acute desynchronization (b). After cessation of stimulation, the network slowly resynchronizes. During the transition, the network remains partially desynchronized for several periods of the original synchronous rhythm.

FIG. 1.

Acute desynchronization by CR stimulation. (a) Raster plot of spiking activity for the integrate-and-fire model without STDP (η=0) (see Sec. V). Stimuli (yellow) are applied to M=4 neuronal subpopulations with fCR=3.8Hz and stimulation amplitude Astim=0.1. y-axis shows neuron indices sorted according to their coordinates. (b) Time-averaged Kuramoto order parameter [Eq. (1)] before, during (yellow), and after cessation of stimulation (Δ=0.5s). Markers indicate evaluation times t of Eq. (1). Stimulation induces acute desynchronization (b). After cessation of stimulation, the network slowly resynchronizes. During the transition, the network remains partially desynchronized for several periods of the original synchronous rhythm.

Close modal

Originally, CR stimulation was delivered on demand, such that short stimulation episodes, as shown in Fig. 1, were applied whenever the system resynchronized.4 In an alternative embodiment, CR stimuli were delivered periodically with stimulus amplitude adapted to the amount of neuronal synchrony.4 Most studies, however, consider intermittent CR, where short stimulation “ON” periods of, say, three CR cycles are interrupted by stimulation “OFF” periods, e.g., by two cycles without stimulus delivery (see, e.g., Ref. 45).

Studies of neuronal networks with STDP showed that acute desynchronization by CR stimulation might entail long-lasting effects, outlasting stimulation.8,40,46 In plastic networks, connections dynamically change according to the activity of connected oscillators. This allows the network structure to adapt and stabilize different dynamical patterns. In particular, states with distinct dynamics, such as synchronized and desynchronized states or cluster states, may coexist.40,47–53 A common plasticity mechanism is STDP, which regulates the weights of synaptic connections according to the time lags between post- and presynaptic spikes.54,55

CR stimulation may alter the statistics of time lags between post- and presynaptic spikes in a way that STDP causes a weakening of synaptic weights. This destabilizes the synchronized state and may drive the network into the attractor of a stable desynchronized state.40 This observation has led to the hypothesis that long-lasting therapeutic effects of CR stimulation result from a stimulation-induced transition between a pathological and a physiological state;40 these states correspond to a synchronized state with strong synaptic connections and a desynchronized state with weak connections, respectively. This is illustrated in Fig. 2.

FIG. 2.

CR stimulation induces a transition to the physiological state. In plastic neuronal networks, stable states with synchronized and desynchronized activity may coexist. These are associated with pathological and physiological neuronal activity. CR stimulation may drive the network into the basin of attraction of a stable desynchronized state, leading to long-lasting desynchronization that persists after cessation of stimulation. This may explain long-lasting therapeutic effects that outlast the duration of CR stimulation.5–7 

FIG. 2.

CR stimulation induces a transition to the physiological state. In plastic neuronal networks, stable states with synchronized and desynchronized activity may coexist. These are associated with pathological and physiological neuronal activity. CR stimulation may drive the network into the basin of attraction of a stable desynchronized state, leading to long-lasting desynchronization that persists after cessation of stimulation. This may explain long-lasting therapeutic effects that outlast the duration of CR stimulation.5–7 

Close modal

Studies on plasticity in the STN are limited; however, it was shown that STN synapses reshape during HF DBS.56–58 This creates evidence that appropriate stimulation may cause structural changes that may contribute to long-lasting therapeutic effects.

Of particular interest is how the number of stimulation sites M impacts desynchronization effects. In clinical DBS, M is closely related to electrode design and placement accuracy.59 The impact of electrode location on therapeutic effects is extensively discussed in the context of HF DBS.60–64 Modern segmented multisite electrodes possess dozens of stimulation contacts.65 So far, CR stimulation was delivered through electrodes with four stimulation contacts.5–7 However, understanding how the number of contacts influences stimulation outcome, might help to exploit modern electrodes for improved therapeutic effects.66 

So far, computational studies analyzed the impact of the number of stimulation sites M on acute desynchronization in all-to-all coupled networks.67 In Ref. 4, larger M was associated with faster decay of the M-cluster state during OFF periods. Further, for sufficiently narrow spatial stimulation profiles (with a sufficiently small spatial decay rate of the stimulation current), acute desynchronization improves with an increasing number of stimulation sites M, saturating at values around M=10 (see Ref. 67). M=4 is typically used in theoretical and computational studies on CR stimulation.4,8,45,46 It is not known how these results translate to complex networks and how the number of stimulation sites correlates with long-lasting effects in plastic networks.

In the present paper, for the first time, we study how the number of stimulation sites affects long-lasting desynchronization by CR stimulation in a neural network with STDP. We perform extensive simulations of excitatory neuronal networks of leaky integrate-and-fire (LIF) neurons with STDP. We consider a one-dimensional setup with a linear alignment of M stimulation sites. In addition, we perform a detailed theoretical analysis of the impact of M on synaptic reshaping. Our theoretical predictions are compared to simulation results. Surprisingly, long-lasting effects of CR stimulation are most pronounced when stimulation parameters are adjusted to the STDP and not to the dominant synchronous rhythm as suggested earlier.4 For stimulation frequencies adapted to STDP parameters, M greatly impacts the dynamics of synaptic weights. While high M may prevent long-lasting desynchronization effects, rendering CR stimulation inefficient, an appropriate number of stimulation sites results in pronounced long-lasting effects, even for short stimulation times.

We consider N=1000 LIF neurons that are randomly placed along the x-axes. Neuron centers xi are uniformly distributed in the interval [2.5,2.5]mm between neurons.68 Excitatory synapses between neurons i and j are randomly added. The connection probability exp(dij/0.5mm) depends on the distance dij=|xjxi|. In total, 7% of all possible non-autaptic connections are added.68 

The subthreshold membrane potential Vi(t) of neuron i obeys

CidVidt=gleak(VrestVi)+gsyn,i(t)(VsynVi)+Istim(t)+Inoise,i(t).
(2)

Ci is the membrane capacitance. Terms on the right-hand side are: the leak current with leakage conductance gleak and resting potential Vrest, the excitatory synaptic input with synaptic conductance gsyn,i(t) and reversal potential Vsyn, the applied stimulation current Istim(t), and the noisy input current Inoise,i(t). Spikes are generated whenever Vi(t) crosses the dynamic threshold Vth,i(t), given by

τthdVth,idt=(Vth,restVth,i).
(3)

We implement a rectangular spike shape by setting Vi(t) to Vspike for a duration of τspike. Afterward, Vth,i(t) and Vi(t) are reset: Vth,i(t)Vth,spike and Vi(t)Vreset.

Dynamic synaptic conductances gsyn,i(t) obey

τsyndgsyn,idt=gsyn,i+κτsynNjGiwji(t)ljδ(ttljjtd).
(4)

τsyn is the synaptic timescale, tljj the timing of the ljth spike of the presynaptic neuron j, and td the synaptic transmission delay. The outer sum runs over all presynaptic neurons Gi of neuron i. κ is the maximal coupling strength. Strengths of individual synapses between presynaptic neuron j and postsynaptic neuron i are scaled by the time-dependent weights wji(t).

Neurons are subject to noisy input Inoise,i(t), modeled by independent Poisson input with firing rate fnoise that is fed into excitatory synapses on neuron i. The resulting current is given by

Inoise,i(t)=gnoise,i(t)(VsynVi),
(5)

with noise conductances gnoise,i(t) resulting from

τsyndgnoise,idt=gnoise,i+κnoiseτsynkiδ(tkit).
(6)

Here, κnoise scales the noise intensity. tki is the kith spike time of the presynaptic Poisson spike train fed into neuron i.

Parameters are set according to Ref. 69: gleak=0.02mS/cm2, Vrest=38mV, Vreset=67mV, Vth,spike=0mV, Vth,rest=40mV, τth=5ms, Vsyn=0mV, τsyn=1ms, td=3ms, κ=8mS/cm2, κnoise=0.026mS/cm2, and fnoise=20Hz. Membrane capacitances Ci are Gaussian distributed [N(μC,σC), μC=3μF/cm2, and σC=0.05μC].

Initially, membrane potentials are uniformly distributed, Vi[Vreset,Vth,rest], and dynamic conductances set to zero. Equations (2)–(6) are integrated using an explicit Euler scheme with time step dt=0.1 ms.

We implement a nearest neighbor STDP scheme by updating the weights wij whenever a postsynaptic neuron spikes, t=tpost, or a presynaptic spike, with spike time tpre, has arrived at the postsynaptic neuron, t=tpre+td. Weight updates are given by wijwij+W(Δt) with time lag Δt=tposttpretd and STDP function

W(t)=η{exp(Δtτ+),Δt>0,0,Δt=0,βτRexp(|Δt|τ),Δt<0.
(7)

η1 scales the weight update per spike, τ+ and τ=τRτ+ are the STDP decay times for long-term potentiation (W>0) and long-term depression (W<0), respectively, and β is the ratio of overall depression 0dtW(t) to potentiation 0dtW(t). In addition, we apply hard bounds wji(t)[0,1].

Throughout the paper, we set β=1.4, τ+=10ms, and τR=4 for which desynchronized and synchronized states coexist.69 We consider slow STDP by setting η=0.02.

In neuroscience and medicine, electrical stimuli are typically charge-balanced to avoid tissue damage.70 Hence, we deliver charge-balanced stimuli to the neurons. Individual stimuli consist of two rectangular current pulses: an excitatory and an inhibitory one. The excitatory one has pulse duration νe=0.4ms and amplitude of Ae=AstimΔVμC/νe and the inhibitory one νi=3ms and amplitude Ai=Aeνe/νi. The inhibitory pulse is separated by 0.2 ms from the excitatory one. Here, ΔV=Vth,spikeVreset. Astim is the dimensionless stimulation strength.

We deliver stimuli according to the CR protocol with rapidly varying sequence (RVS CR) introduced in Refs. 45 and 71 for M=4 stimulation sites and used in earlier preclinical and clinical studies.5–7 Here, we consider a generalized version with M2. The CR frequency fCR is realized by delivering stimuli at times t=k/MfCR, k=0,1, after stimulation onset. Each stimulus is delivered to one out of M neuronal subpopulations. The mth subpopulation contains neurons with xi[2.5+5m/M,2.5+5(m+1)/M[mm. RVS CR is delivered in cycles of M stimuli. The sequence of stimulated subpopulations is randomly drawn, such that each subpopulation receives exactly one stimulus per cycle.45,71

We perform long simulations for different initial distributions of synaptic weights wij(t). In particular, we randomly set wij(t) either to one or zero such that a given mean weight w(t)=ijwij(t) is realized. We evaluate w(t) and the degree of synchrony as quantified by the time-averaged Kuramoto order parameter ρ¯Δ(t) [Eq. (1)].

Figure 3 depicts the results for ρ¯Δ(t) and w(t). We find that initial networks with high mean synaptic weight approach a state with synchronized activity and strong synaptic connections. In contrast, initial networks with low mean synaptic weight approach a desynchronized state with weak connections.

FIG. 3.

Coexistence of synchronized and desynchronized state. Simulation results for different mean initial weights (gray scale). (a) Time-averaged Kuramoto order parameter ρ¯Δ [Eq. (1)], with Δ=2s. (b) Mean synaptic weight w(t). Networks with low initial mean weights (light) approach a desynchronized state with weak connections. In contrast, networks with high initial mean weights (dark) approach a synchronized state with strong connections. Parameters: w(t=0)=0.1,0.2,0.3, and 0.5 (from bottom to top).

FIG. 3.

Coexistence of synchronized and desynchronized state. Simulation results for different mean initial weights (gray scale). (a) Time-averaged Kuramoto order parameter ρ¯Δ [Eq. (1)], with Δ=2s. (b) Mean synaptic weight w(t). Networks with low initial mean weights (light) approach a desynchronized state with weak connections. In contrast, networks with high initial mean weights (dark) approach a synchronized state with strong connections. Parameters: w(t=0)=0.1,0.2,0.3, and 0.5 (from bottom to top).

Close modal

In the following, we prepare the network in the synchronized state, by selecting trajectories with w(t=0)=0.5 and simulating for 500 s until the mean weight becomes stationary (see Fig. 3). Then, we apply RVS CR stimulation for 500 s. After cessation of stimulation, we simulate the network for 1000 s in order to test for long-lasting effects.

We evaluate traces of the time-averaged Kuramoto order parameter ρ¯Δ(t) and the mean synaptic weight w(t), before, during, and after cessation of stimulation.

Figure 4 shows results for ρ¯Δ(t) and w(t) for different numbers of stimulation sites M and strong stimulation Astim=1. We find that the stimulation-induced weight dynamic strongly depends on the number of stimulation sites M. While RVS CR with small M=2 and rather large M=24 results in a weakening of synaptic weights, stimulation with M=12 strengthens synapses. Accordingly, while long-lasting desynchronization is observed for the former values of M, the network returns to the synchronized state and stimulation effects disappear shortly after stimulation ceases for the latter (M=12).

FIG. 4.

Long-lasting stimulation effects depend on number of stimulation sites. (a) Kuramoto order parameter ρ¯Δ [Eq. (1)], time-averaged over time windows of Δ=2s. (b) Mean synaptic weight, w for three numbers of stimulation sites and no-stimulation control. The stimulation period is marked yellow. Panels (a) and (b) show parts of (a) and (b) right after stimulation onset. Black dotted lines mark theoretical predictions for the mean weight w(t) based on Eq. (10) and the assumption that weights are close to the hardbounds at stimulation onset. Parameters: Astim=1, fCR=17.5Hz.

FIG. 4.

Long-lasting stimulation effects depend on number of stimulation sites. (a) Kuramoto order parameter ρ¯Δ [Eq. (1)], time-averaged over time windows of Δ=2s. (b) Mean synaptic weight, w for three numbers of stimulation sites and no-stimulation control. The stimulation period is marked yellow. Panels (a) and (b) show parts of (a) and (b) right after stimulation onset. Black dotted lines mark theoretical predictions for the mean weight w(t) based on Eq. (10) and the assumption that weights are close to the hardbounds at stimulation onset. Parameters: Astim=1, fCR=17.5Hz.

Close modal

In the following, we study the influence of RVS CR stimulation on the synaptic weight dynamics in more detail. To this end, we expand theoretical results on the weight dynamics during RVS CR with four stimulation sites from Ref. 69 to arbitrary M.

Following Ref. 69, we consider the dynamics of a single synaptic weight wij(t) with presynaptic neuron i and postsynaptic neuron j. The dynamics of wij(t) results from the STDP scheme [Eq. (7)] and is determined by the statistics of time lags between post- and presynaptic spikes.72–74 The average rate of weight change Jij(t,T) in a time interval [t,t+T] is given by

Jij(t,T)=1Tti,tjpairsW(tj(ti+td)).
(8)

The sum runs over all pairs of spike times ti and tj that contribute to weight changes in the interval.

We assume that stimulation is sufficiently strong and fast so that spiking occurs mainly in response to stimuli. This applies if the stimulation frequency is larger than that of the dominant synchronous rhythm, fdom3.5Hz, and the stimulation strength, Astim, is close to one. Then, neuronal spike trains can be approximated by

xl(t)=k,skSlδ(t(sk+ϵk,l)).
(9)

Here, δ(t) is the Dirac delta distribution, Sl are the stimuli sk delivered to neuron l, and ϵk,ls are time lags between spikes and the stimuli that trigger them.

For CR stimulation, Jij(t,T) saturates and becomes independent of t when T becomes long compared to 1/fstim, Jij(t,T)Jij. We consider ϵk,l to be distributed according to a distribution λ(ϵk,l). Using Eq. (9) in Eq. (8) and averaging over realizations of the stimulation sequence and spike times yields

JijfCRdtGij(t)W(ttd).
(10)

Gij(t) is the distribution of time lags between post- and presynaptic spikes that contribute to weight updates. We derive analytical results for Gij(t) in the following and analyze its dependence on stimulation parameters.

Gij(t) results from the statistics of interstimulus intervals between stimuli skj and ski applied to the post- and presynaptic neuron, respectively. We introduce the probability ppairij(s,ϵ) that the post- and presynaptic neurons receive stimuli with interstimulus interval s=skjski and the resulting spike pair contributes to weight updates. ϵ=ϵkj,jϵki,i is the difference in response times to stimuli skj and ski [see Eq. (9)].

RVS CR stimulation is delivered to M neuronal subpopulations. Based on the subpopulations of post- and presynaptic neurons, we distinguish between two groups of synapses. Synapses of the first group connect neurons that are part of the same subpopulation, while the second group contains synapses for which post- and presynaptic neurons belong to different subpopulations. We will refer to the first group as intrapopulation synapses and to the second group as interpopulation synapses in the following. In order to distinguish between respective quantities, we use the suffixes “intra” and “inter” instead of those referring to neuron indices.

In Ref. 69, the case M=4 was studied. However, their results apply to general M as long as interstimulus intervals are long compared to the synaptic transmission delay td. Their results are given in the supplementary material. As 1/MfCR determines interstimulus intervals for RVS CR, their assumptions may not hold for RVS CR with large numbers of stimulation sites M.

Here, we extend the result of Ref. 69 to the case of 1/MfCR<td<2/MfCR, i.e., spikes of a stimulated subpopulation do not arrive at postsynaptic neurons before the next subpopulation is stimulated. This rearranges the order of postsynaptic spikes and presynaptic spike arrival times and leads to different post-pre-pairings. We describe the change of ppairA(s,ϵ) by introducing a correction term δppairA(s,ϵ) such that ppairA(s,ϵ)=ppair,0A(s,ϵ)+δppairA(s,ϵ). Suffix “A” refers to “intra” or “inter” synapses. ppair,0A(s,ϵ) is the result of Ref. 69 and given in the supplementary material. Note that ppairA(s,ϵ) and ppair,0A(s,ϵ) are normalized to two, i.e., on average two pairings per spike result in weight updates. In contrast, δppairA(s,ϵ) is normalized to zero and can attain negative values.

In Ref. 69, ppair,0A(s,ϵ) was derived by considering stimulation times for the postsynaptic neuron for each of the M possible stimulation times of the presynaptic neuron within a CR cycle and evaluating the resulting time lags. δppairA(s,ϵ) can be derived by considering stimulation times of the postsynaptic neuron for each CR sequence in which the postsynaptic neuron receives a stimulus right after the presynaptic one. Then, resulting spike pairings for 1/MfCR<td<2/MfCR are compared with the case td<1/MfCR. The difference of the distributions of interstimulus times s for these two cases yields the correction term δppairA(s,ϵ). For intrapopulation synapses, we find

δppairintra(s,ϵ)=1M2δ(s)+ξ=0M1δ(sM+1+ξfCRM)M3,1MfCR+ϵ<td.
(11)

For interpopulation synapses, we find

δppairinter(s,ϵ)=u=0M2v=0M1δ(sBuvfCRM)δ(s+CuvfCRM)M2(M1)+v=0M1δ(sDuvfCRM)M3v=1M1δ(s+vfCRM)M2(M1),1MfCR+ϵ<td,
(12)

with Buv=Mu+v, Cuv=1+u+v, and Duv=M+u+v. Note that δppairA(s,ϵ)=0 for 1/MfCR+ϵtd.

The resulting distribution of time lags GA(t) is given by

GA(t)=dsΛ(ts)ppairA(s,ts).
(13)

Here, Λ(t)=dtλ(t)λ(t+t). Ginter(t) and Gintra(t) for different values of M are shown in Fig. 5. As M increases, more peaks appear at multiples of 1/MfCR. Note that peaks appear at 1/MfCR<τd for large M (see panels for M=24).

FIG. 5.

Distribution of time lags between post- and presynaptic spikes. (a) Shape of STDP function shifted by delay time td. (b) and (c) log[Ginter(t)] (b) and log[Gintra(t)] (c) for different M obtained from Eqs. (13), (11), and (12). Red region marks time lags that lead to long-term potentiation [see Eq. (10)]. Parameters: fCR=17.5Hz and λ(t) is approximated by a Gaussian distribution with zero mean and standard deviation 0.1ms.

FIG. 5.

Distribution of time lags between post- and presynaptic spikes. (a) Shape of STDP function shifted by delay time td. (b) and (c) log[Ginter(t)] (b) and log[Gintra(t)] (c) for different M obtained from Eqs. (13), (11), and (12). Red region marks time lags that lead to long-term potentiation [see Eq. (10)]. Parameters: fCR=17.5Hz and λ(t) is approximated by a Gaussian distribution with zero mean and standard deviation 0.1ms.

Close modal

The effect of increasing M and fCR on Jinter and Jintra can be studied by using GA(t) in Eq. (10). Results are shown in Figs. 6(a) and 6(b), for inter- and intrapopulation synapses, respectively. While Jintra decays monotonically with increasing CR frequency fCR, Jinter possesses a complex nonlinear dependence on both fCR and M [see Fig. 6(a)]. Therefore, we will focus our analysis on the dynamics of interpopulation synapses. As illustrated in Fig. 5(b), Ginter(t) consists of peaks at integer multiples of 1/MfCR. For 1/MfCRτ+, corresponding time lags hardly contributed to weight updates. As M increases, new peaks appear and corresponding time lags lead to synaptic potentiation. However, for the first M with 1/MfCR<τd, a peak appears at time lags that contribute to synaptic depression. As M increases further, this behavior repeats with the kth peaks leading to synaptic depression if k/MfCR<τd. This leads to the complex shape of Jinter [Fig. 6(a)], with “tongue”-like regions of synaptic depression and potentiation, respectively.

FIG. 6.

Stimulation-induced synaptic weight dynamics and long-lasting desynchronization. (a) and (b) Theoretical predictions for Jinter (a) and Jintra (b). Gray curve marks Jinter=0. Dotted vertical line marks frequency of dominant rhythm fdom in synchronized state. Stars mark parameter values used in Fig. 4. Dashed gray lines in (a) mark parameter combinations with k/MfCR=τd, k=1,2,3,4 from left bottom to right top. (c) and (d) Theoretical predictions for Jinter and Jintra compared to numerical simulations (symbols) for different CR frequencies and Astim=1. (e) Acute effect on mean weight wAC when stimulation is turned off. Panels show Astim=1, 0.5, and 0.1 from left to right. Gray contour curve from panel (a) is added. (f) Time-averaged Kuramoto order parameter ρ¯Δ evaluated 1000 s after cessation of stimulation for a time window of 2 s. Values for Astim as in (e). Parameters: Stimulation is delivered for 500 s after 500 s of equilibration time. Simulation results in panels (c)–(f) are averaged over three network realizations and three realizations of CR sequences per network realization.

FIG. 6.

Stimulation-induced synaptic weight dynamics and long-lasting desynchronization. (a) and (b) Theoretical predictions for Jinter (a) and Jintra (b). Gray curve marks Jinter=0. Dotted vertical line marks frequency of dominant rhythm fdom in synchronized state. Stars mark parameter values used in Fig. 4. Dashed gray lines in (a) mark parameter combinations with k/MfCR=τd, k=1,2,3,4 from left bottom to right top. (c) and (d) Theoretical predictions for Jinter and Jintra compared to numerical simulations (symbols) for different CR frequencies and Astim=1. (e) Acute effect on mean weight wAC when stimulation is turned off. Panels show Astim=1, 0.5, and 0.1 from left to right. Gray contour curve from panel (a) is added. (f) Time-averaged Kuramoto order parameter ρ¯Δ evaluated 1000 s after cessation of stimulation for a time window of 2 s. Values for Astim as in (e). Parameters: Stimulation is delivered for 500 s after 500 s of equilibration time. Simulation results in panels (c)–(f) are averaged over three network realizations and three realizations of CR sequences per network realization.

Close modal

We compare theoretical predictions for Jinter and Jintra with simulation results in Figs. 6(c) and 6(d). We find good quantitative agreement for strong stimulation Astim=1.

Next, we perform extensive simulations for different fCR and M, with stimulation time as in Fig. 4. Results for the mean weight shortly before cessation of stimulation wAC are depicted in Fig. 6(e) for Astim=1, 0.5, and 0.1 corresponding to strong, moderate, and weak stimulation, respectively. Note that for very small stimulation strengths (Astim0), stimulation hardly affects neuronal spiking and the network remains in the synchronized state (data not shown). We find that the characteristic shape of Jinter translates into regions of weak synaptic weights (for Jinter<0) and strong synaptic weights (for Jinter>0) at the end of the stimulation period. For moderate and weak stimulation, parts of this pattern remain; however, we find that high fCR and M lead to synaptic weakening rather than strengthening.

We evaluate the time-averaged Kuramoto order parameter [Eq. (1)], 1000 s after cessation of stimulation. Results are shown in Fig. 6(f). The complex shape of Jinter translates into distinct regions of long-lasting desynchronization (for Jinter<0) and synchronization (for Jinter>0), respectively. Remarkably, moderate or weak RVS CR stimulation with high stimulation frequencies and large numbers of stimulation sites causes long-lasting desynchronization. In contrast, strong stimulation with the same parameters is inefficient [see Fig. 6(f)].

We studied CR stimulation of excitatory recurrent neuronal networks with STDP. We focused on the impact of the CR frequency and the number of stimulation sites on long-lasting desynchronization effects, caused by synaptic depression during stimulation. Our theoretical analysis and simulations revealed pronounced synaptic reshaping for CR frequencies in the range of inverse STDP decay times, thus, higher than that of the dominant synchronous rhythm. Synapses are reshaped differently, depending on whether they connect neurons close to the same or to different stimulation sites. Remarkably, for synapses of the latter type, regions of depression and potentiation form an alternating pattern in parameter space as MfCR increases. Therefore, specific parameter combinations are needed to achieve long-lasting effects. For moderate and weak stimulation, parts of the pattern are still observable; however, we find well-pronounced long-lasting desynchronization for high fCR and M. Our results may assist parameter selection for CR stimulation in preclinical and clinical studies.

We consider CR with rapidly varying sequence (RVS CR), which was used in preclinical and clinical studies.5–7 An earlier computational study by Lysyansky et al. on the effect of the number of stimulation sites M on acute desynchronization,67 considered spatially selective stimulation profiles and a fixed CR sequence with and without ON:OFF intermittent pattern. Stimulation-induced desynchronization by CR with narrow profile and without intermittent pattern improved up to M=10 and then saturated. We studied RVS CR without ON:OFF intermittent pattern. We find that acute desynchronization effects of strong stimulation are independent of the number of stimulation sites (see Fig. 4). This is because, in subsequent cycles, the same subpopulation receives resetting stimuli at different phases due to shuffling of the stimulation sequence. Hence, their oscillation periods vary among CR cycles. This results in oscillations of the Kuramoto order parameter that result in a non-zero time average (see Fig. 4). This differs from the finding in Ref. 67. The difference might result from the fact that Lysyansky et al. considered a fixed CR sequence and a spatially selective stimulation profile. Note that we observed acute desynchronization for fixed CR sequence (see Fig. 1).

We focus on long-lasting desynchronization effects that result from stimulation-induced synaptic depression. We reveal that long-lasting effects on interpopulation synapses possess a complex nonlinear dependence on the CR frequency and number of stimulation sites [see Fig. 6(a)]. In parameter space, this results in an alternating pattern of regions with pronounced long-lasting desynchronization and regions where stimulation effects decay shortly after cessation of stimulation. In contrast, synaptic depression of intrapopulation synapses shows only a weak dependence on the number of stimulation sites and becomes more pronounced at high stimulation frequencies [Fig. 6(b)]. Parameter regions with long-lasting desynchronization and synchronization coincide with regions of stimulation-induced synaptic depression and potentiation of interpopulation synapses, respectively (see Fig. 6). For strong stimulation, regions of synaptic depression and potentiation can be predicted by our theory [see Figs. 6(c) and 6(d)]. Note that theoretical predictions for low numbers of stimulation sites still hold for weaker stimulation [Fig. 6(f)]. For very low frequencies, neuronal spiking is dominated by synaptic input rather than stimulation, which leads to deviations. We find pronounced long-lasting desynchronization for weak stimulation with high CR frequency and large number of stimulation sites. For such stimulation, stimuli arrive at short time lags; however, neurons cannot respond to weak stimuli when they are far from the spiking threshold. We hypothesize that weak stimulation acts like noise and induces random spiking, which in turn results in weakening of synaptic weights and long-lasting desynchronization [see Figs. 6(e) and 6(f)].

The alternating pattern (see Fig. 6) results from the interplay of two competing effects: first, as the product MfCR increases, small time lags between separately stimulated subpopulations occur more frequently (see Fig. 5). As potentiation dominates over depression for short time lags, this increases the contribution of synaptic potentiation to the stimulation-induced dynamics of interpopulation synapses. Second, a delay-induced effect causes synaptic depression for strongly synchronized neuronal subpopulations.75,76 If time lags between spikes become shorter than the delay time, presynaptic spikes arrive after postsynaptic ones. This results in pronounced synaptic depression.75,76 For RVS CR stimulation, this affects interpopulation weights if time lags between spiking of separately stimulated subpopulations are small. In particular, whenever integer multiples of 1/MfCR become smaller than the delay time td, the contribution of synaptic depression to the dynamics of interpopulation synapses increases. The interplay of both potentiation for large MfCR and the increase of delay-induced depression leads to the alternating pattern of parameter regions with pronounced synaptic depression and potentiation, respectively.

Our results may have important consequences for parameter adjustment in clinical CR stimulation for the treatment of Parkinson’s disease. Here, symptoms such as tremor, bradykinesia, and rigidity are associated with pronounced synchronization in certain frequency bands, such as theta (3--10Hz) (tremor)77–79 and beta band (13--30Hz) (bradykinesia and rigidity).80,81 The CR frequency is typically adjusted to the dominant peak in the power spectrum.5–7 In contrast, our results on long-lasting effects due to stimulation-induced synaptic depression suggest that adjusting the CR frequency to the inverse STDP time scale may yield better long-lasting aftereffects (see Figs. 5 and 6). This time scale is typically in the range of 10–100 ms.82 Consequently, regions with strongly pronounced synaptic depression, Jintra/inter<0, are mainly located at CR frequencies that are significantly higher than that of the dominant synchronous rhythm fdom (see Fig. 6). In a previous computational study for M=4, it was found that fCRfdom yields most pronounced long-lasting desynchronization.8 The authors considered a network of excitatory and inhibitory neurons, with STDP affecting all synapses. Thus, desynchronization results from the interplay of excitatory and inhibitory neurons. While our results on stimulation-induced weight dynamics can be applied to their model, our mechanism for desynchronization (weak excitatory connections) differs significantly from theirs. This may be the reason for the different results.

Furthermore, the complex nonlinear dependence of long-lasting effects and synaptic depression on stimulation parameters (see Fig. 6) raises the question whether aftereffects of CR stimulation are robust with respect to further aspects of the treatment such as electrode placement and additional medication. Electrode placement determines which neurons receive stimuli. Furthermore, it affects whether stimuli are delivered directly to the somata or as sensory input via fibers. An earlier computational study found pronounced long-lasting effects of CR stimulation for both direct electrical and sensory stimulation. They concluded that the effects of CR are robust with respect to electrode placement.46 While Ref. 46 considered the same STDP function for both types of stimulation, its shape may strongly depend on the type of neuron and how neurons receive stimuli.82 For instance, even the location of the synapse receiving input can affect the shape of the STDP function.83 We showed that the adjustment of the CR frequency to the STDP decay times is critical for synaptic depression and resulting long-lasting effects. Consequently, long-lasting effects of sensory CR stimulation may differ from those of direct stimulation.

Stimulation-induced reshaping of excitatory synapses in the STN also depends on the concentration of dopamine.57 Medication might, therefore, change the STDP function governing the evolution of synaptic weights and, consequently, the parameter region for which synaptic depression is observed in Fig. 6. Our theoretical and computational results highlight the importance of the number of stimulation sites, the CR frequency, and the properties of synaptic plasticity for long-lasting aftereffects of CR stimulation. In a clinical setup, however, information about plasticity in target brain regions is often unavailable. Therefore, modern multisite stimulation electrodes with directional steering might by advantageous to exploit STDP for long-lasting effects. In these electrodes, independent stimulus delivery to a large number of stimulation contacts is used to generate spatial current profiles. This could be used to scan for suitable combinations of stimulation parameters to obtain the desired therapeutic effects.65,66,84 So far, CR stimulation was tested in rat hippocampal slices,85 in the 1-methyl-4-phenyl-1,2,3,6-tetrahydropyridine Parkinsonian monkey model5,7 , and in Parkinson’s disease patients,6 where up to four stimulation contacts were used. Tass et al.85 used an array of four different unipolar stimulation electrodes. As yet, the preclinical and clinical CR Deep Brain Stimulation (CR-DBS) studies5–7 were performed with electrodes with four cylindrical contacts where stimulation was delivered through the three lower contacts, whereas the upper contact served as common reference. Recently, directional DBS electrodes with radially segmented contacts were developed to shape and direct the electrical field in the horizontal plane.65,86 When delivering conventional HF DBS, segmented electrodes may improve the therapeutic window, by enhancing efficacy and reducing side effects.65 Furthermore, using moderate stimulation amplitudes, segmented electrodes may enable to direct the electrical field to different targets, in this way providing more options for multisite stimulation protocols, hence, potentially improving the outcome of CR-DBS.86 Optimal steering of the electrical field for improved multichannel stimulation protocols, such as CR-DBS, will benefit from biophysically realistic neuronal network models. However, as a first step, to set expectations and to develop adequate multisite protocols, for the first time, we here considered the impact of the number of stimulation sites on the long-lasting effects of CR stimulation. Intriguingly, our study indicates that in the presence of STDP, the relation between long-lasting desynchronization effects and the number of stimulation sites is actually complex (see Fig. 6). This is somewhat counterintuitive as one would expect that stimulation patterns that separate a neuronal population into a large number of subpopulations are better suited for desynchronization.67 

Regarding standard HF DBS, so far, the mechanism of action of convention HF DBS is a matter of debate (see, e.g., Refs. 87–89). The blockage of nerve conduction induced by HF DBS is just one hypothesis (see references above and Ref. 90). As yet, the computational studies of CR stimulation led to non-trivial hypotheses, such as long-lasting and cumulative effects,40,41 that were verified experimentally.5–7 Permanent HF electrical bursts delivered through electrodes may have different effects, depending on the activated target structure. However, brief HF bursts of sufficient strength, delivered through the soma and/or through excitatory and/inhibitory synapses may reliably cause a phase reset, so that the corresponding CR stimulus patterns induce long-lasting desynchronization.46 This fundamental dynamical property of CR stimulation argues in favor of the robustness of CR stimulation. However, in principle, alternative mechanisms of CR stimulation will have to be taken into account and explored in future studies.

In future studies, we anticipate using detailed models of Parkinson’s disease-related brain regions and spatial stimulation profiles to study multisite CR stimulation. These studies may yield insight into how robust the observed long-lasting effects are with respect to more complex intrinsic neuronal dynamics and input from other brain regions. We hope that our results ultimately motivate preclinical and clinical studies on multisite stimulation techniques that exploit neuronal plasticity for long-lasting therapeutic effects.

See the supplementary material for results on ppair,0intra(s,ϵ) and ppair,0inter(s,ϵ) published in Ref. 69.

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

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1.
A.
Schnitzler
,
C.
Münks
,
M.
Butz
,
L.
Timmermann
, and
J.
Gross
, “
Synchronized brain network associated with essential tremor as revealed by magnetoencephalography
,”
Mov. Disord.
24
,
1629
1635
(
2009
).
2.
C.
Hammond
,
H.
Bergman
, and
P.
Brown
, “
Pathological synchronization in Parkinson’s disease: Networks, models and treatments
,”
Trends Neurosci.
30
,
357
364
(
2007
).
3.
P.
Temperli
,
J.
Ghika
,
J.-G.
Villemure
,
P.
Burkhard
,
J.
Bogousslavsky
, and
F.
Vingerhoets
, “
How do parkinsonian signs return after discontinuation of subthalamic DBS?
,”
Neurology
60
,
78
81
(
2003
).
4.
P. A.
Tass
, “
A model of desynchronizing deep brain stimulation with a demand-controlled coordinated reset of neural subpopulations
,”
Biol. Cybern.
89
,
81
88
(
2003
).
5.
P. A.
Tass
,
L.
Qin
,
C.
Hauptmann
,
S.
Dovero
,
E.
Bezard
,
T.
Boraud
, and
W. G.
Meissner
, “
Coordinated reset has sustained aftereffects in parkinsonian monkeys
,”
Ann. Neurol.
72
,
816
820
(
2012
).
6.
I.
Adamchic
et al., “
Coordinated reset neuromodulation for Parkinson’s disease: Proof-of-concept study
,”
Mov. Disord.
29
,
1679
1684
(
2014
).
7.
J.
Wang
,
S.
Nebeck
,
A.
Muralidhara
,
M. D.
Johnson
,
J. L.
Vitek
, and
K. B.
Baker
, “
Coordinated reset deep brain stimulation of subthalamic nucleus produces long-lasting, dose-dependent motor improvements in the 1-methyl-4-phenyl-1, 2, 3, 6-tetrahydropyridine non-human primate model of parkinsonism
,”
Brain Stimul.
9
,
609
617
(
2016
).
8.
T.
Manos
,
M.
Zeitler
, and
P. A.
Tass
,
PLoS Comput. Biol.
14
,
e1006113
(
2018
).
9.
A. T.
Winfree
, “
Biological rhythms and the behavior of populations of coupled oscillators
,”
J. Theor. Biol.
16
,
15
42
(
1967
).
10.
Y.
Kuramoto
, “Self-entrainment of a population of coupled non-linear oscillators,” in International Symposium on Mathematical Problems in Theoretical Physics, edited by H. Araki (Springer, New York, 1975), pp. 420–422.
11.
Y.
Kuramoto
,
Chemical Oscillations, Waves, and Turbulence
(
Springer
,
Berlin
,
1984
).
12.
N.
Kopell
, “Toward a theory of modelling central pattern generators,” in Neural Control of Rhythmic Movements in Vertebrates, edited by A. H. Cohen, S. Rossignol, and S. Grillner (Wiley, New York, 1988), pp. 369–413.
13.
R. E.
Mirollo
and
S. H.
Strogatz
, “
Synchronization of pulse-coupled biological oscillators
,”
SIAM J. Appl. Math.
50
,
1645
1662
(
1990
).
14.
G. B.
Ermentrout
and
N.
Kopell
, “
Multiple pulse interactions and averaging in systems of coupled neural oscillators
,”
J. Math. Biol.
29
,
195
217
(
1991
).
15.
D.
Hansel
,
G.
Mato
, and
C.
Meunier
, “
Phase dynamics for weakly coupled Hodgkin–Huxley neurons
,”
Europhys. Lett.
23
,
367
(
1993
).
16.
P. A.
Tass
,
Phase Resetting in Medicine and Biology: Stochastic Modelling and Data Analysis
(
Springer
,
Berlin
,
1999
).
17.
A.
Pikovsky
,
M.
Rosenblum
, and
J.
Kurths
,
Synchronization: A Universal Concept in Nonlinear Sciences
(
Cambridge University Press
,
Cambridge
,
2001
).
18.
F.
Mormann
,
K.
Lehnertz
,
P.
David
, and
C. E.
Elger
, “
Mean phase coherence as a measure for phase synchronization and its application to the EEG of epilepsy patients
,”
Physica D
144
,
358
369
(
2000
).
19.
P.
Krack
,
A.
Batir
,
N.
Van Blercom
,
S.
Chabardes
,
V.
Fraix
,
C.
Ardouin
,
A.
Koudsie
,
P. D.
Limousin
,
A.
Benazzouz
,
J. F.
LeBas
et al., “
Five-year follow-up of bilateral stimulation of the subthalamic nucleus in advanced Parkinson’s disease
,”
N. Engl. J. Med.
349
,
1925
1934
(
2003
).
20.
M. C.
Rodriguez-Oroz
,
J.
Obeso
,
A.
Lang
,
J.-L.
Houeto
,
P.
Pollak
,
S.
Rehncrona
,
J.
Kulisevsky
,
A.
Albanese
,
J.
Volkmann
,
M.
Hariz
et al., “
Bilateral deep brain stimulation in Parkinson’s disease: A multicentre study with 4 years follow-up
,”
Brain
128
,
2240
2249
(
2005
).
21.
M.
Bin-Mahfoodh
,
C.
Hamani
,
E.
Sime
, and
A. M.
Lozano
, “
Longevity of batteries in internal pulse generators used for deep brain stimulation
,”
Stereotact. Funct. Neurosurg.
80
,
56
60
(
2003
).
22.
G. R.
Mines
, “
On circulating excitations in heart muscle and their possible relation to tachycardia and fibrillation
,”
Trans. R. Soc. Can.
8
,
43
52
(
1914
).
23.
A. T.
Winfree
, “
Phase control of neural pacemakers
,”
Science
197
,
761
763
(
1977
).
24.
A. T.
Winfree
,
The Geometry of Biological Time
(
Springer
,
New York
,
1980
).
25.
E.
Warman
and
D.
Durand
, “Desynchronization of epileptiform activity by phase resetting,” in Images of the Twenty-First Century. Proceedings of the Annual International Engineering in Medicine and Biology Society (IEEE, 1989), pp. 1286–1287.
26.
P.
Tass
, “
Resetting biological oscillators—A stochastic approach
,”
J. Biol. Phys.
22
,
27
64
(
1996
).
27.
P. A.
Tass
, “
Stochastic phase resetting: A theory for deep brain stimulation
,”
Prog. Theor. Phys. Suppl.
139
,
301
313
(
2000
).
28.
Y.
Zhai
,
I. Z.
Kiss
,
P. A.
Tass
, and
J. L.
Hudson
, “
Desynchronization of coupled electrochemical oscillators with pulse stimulations
,”
Phys. Rev. E
71
,
065202(R)
(
2005
).
29.
V.
Arnold
, “
Cardiac arrhythmias and circle mappings
,”
Chaos
1
,
20
24
(
1991
).
30.
C. J.
Wilson
,
B.
Beverlin
, and
T.
Netoff
, “
Chaotic desynchronization as the therapeutic mechanism of deep brain stimulation
,”
Front. Syst. Neurosci.
5
,
50
(
2011
).
31.
A. B.
Holt
,
D.
Wilson
,
M.
Shinn
,
J.
Moehlis
, and
T. I.
Netoff
, “
Phasic burst stimulation: A closed-loop approach to tuning deep brain stimulation parameters for Parkinson’s disease
,”
PLoS Comput. Biol.
12
,
e1005011
(
2016
).
32.
M.
Rosenblum
and
A.
Pikovsky
, “
Delayed feedback control of collective synchrony: An approach to suppression of pathological brain rhythms
,”
Phys. Rev. E
70
,
041904
(
2004
).
33.
M. G.
Rosenblum
and
A. S.
Pikovsky
, “
Controlling synchronization in an ensemble of globally coupled oscillators
,”
Phys. Rev. Lett.
92
,
114102
(
2004
).
34.
O. V.
Popovych
,
C.
Hauptmann
, and
P. A.
Tass
, “
Effective desynchronization by nonlinear delayed feedback
,”
Phys. Rev. Lett.
94
,
164102
(
2005
).
35.
Y.
Zhai
,
I. Z.
Kiss
, and
J. L.
Hudson
, “
Control of complex dynamics with time-delayed feedback in populations of chemical oscillators: Desynchronization and clustering
,”
Ind. Eng. Chem. Res.
47
,
3502
3514
(
2008
).
36.
O. V.
Popovych
,
B.
Lysyansky
, and
P. A.
Tass
, “
Closed-loop deep brain stimulation by pulsatile delayed feedback with increased gap between pulse phases
,”
Sci. Rep.
7
,
1033
(
2017
).
37.
O. V.
Popovych
,
B.
Lysyansky
,
M.
Rosenblum
,
A.
Pikovsky
, and
P. A.
Tass
, “
Pulsatile desynchronizing delayed feedback for closed-loop deep brain stimulation
,”
PLoS One
12
,
e0173363
(
2017
).
38.
O. V.
Popovych
and
P. A.
Tass
, “
Adaptive delivery of continuous and delayed feedback deep brain stimulation—A computational study
,”
Sci. Rep.
9
,
10585
(
2019
).
39.
P. A.
Tass
, “
Effective desynchronization by means of double-pulse phase resetting
,”
Europhys. Lett.
53
,
15
(
2001
).
40.
P. A.
Tass
and
M.
Majtanik
, “
Long-term anti-kindling effects of desynchronizing brain stimulation: A theoretical study
,”
Biol. Cybern.
94
,
58
66
(
2006
).
41.
C.
Hauptmann
and
P.
Tass
, “
Cumulative and after-effects of short and weak coordinated reset stimulation: A modeling study
,”
J. Neural Eng.
6
,
016004
(
2009
).
42.
V. A.
Pliss
, “
A reduction principle in the theory of stability of motion
,”
Izv. Akad. Nauk SSSR Ser. Mat.
28
,
1297
1324
(
1964
).
43.
H.
Haken
,
Synergetics: An Introduction
(
Springer
,
Berlin
,
1977
).
44.
P. A.
Tass
, “
Desynchronization by means of a coordinated reset of neural sub-populations—A novel technique for demand-controlled deep brain stimulation
,”
Prog. Theor. Phys. Suppl.
150
,
281
296
(
2003
).
45.
B.
Lysyansky
,
O. V.
Popovych
, and
P. A.
Tass
, “
Desynchronizing anti-resonance effect of m:n ON–OFF coordinated reset stimulation
,”
J. Neural Eng.
8
,
036019
(
2011
).
46.
O. V.
Popovych
and
P. A.
Tass
, “
Desynchronizing electrical and sensory coordinated reset neuromodulation
,”
Front. Hum. Neurosci.
6
,
58
(
2012
).
47.
P.
Seliger
,
S. C.
Young
, and
L. S.
Tsimring
, “
Plasticity and learning in a network of coupled phase oscillators
,”
Phys. Rev. E
65
,
041906
(
2002
).
48.
D. H.
Zanette
and
A. S.
Mikhailov
, “
Dynamical clustering in oscillator ensembles with time-dependent interactions
,”
Europhys. Lett.
65
,
465
471
(
2004
).
49.
M. G.
Zimmermann
,
V. M.
Eguíluz
, and
M.
San Miguel
, “
Coevolution of dynamical states and interactions in dynamic networks
,”
Phys. Rev. E
69
,
065102(R)
(
2004
).
50.
J. M.
Pacheco
,
A.
Traulsen
, and
M. A.
Nowak
, “
Coevolution of strategy and structure in complex networks with dynamical linking
,”
Phys. Rev. Lett.
97
,
258103
(
2006
).
51.
Y. L.
Maistrenko
,
B.
Lysyansky
,
C.
Hauptmann
,
O.
Burylko
, and
P. A.
Tass
, “
Multistability in the Kuramoto model with synaptic plasticity
,”
Phys. Rev. E
75
,
066207
(
2007
).
52.
T.
Aoki
and
T.
Aoyagi
, “
Co-evolution of phases and connection strengths in a network of phase oscillators
,”
Phys. Rev. Lett.
102
,
034101
(
2009
).
53.
R.
Berner
,
J.
Sawicki
, and
E.
Schöll
, “
Birth and stabilization of phase clusters by multiplexing of adaptive networks
,”
Phys. Rev. Lett.
124
,
088301
(
2020
).
54.
H.
Markram
,
J.
Lübke
,
M.
Frotscher
, and
B.
Sakmann
, “
Regulation of synaptic efficacy by coincidence of postsynaptic APs and EPSPs
,”
Science
275
,
213
215
(
1997
).
55.
G.-Q.
Bi
and
M.-M.
Poo
, “
Synaptic modifications in cultured hippocampal neurons: Dependence on spike timing, synaptic strength, and postsynaptic cell type
,”
J. Neurosci.
18
,
10464
10472
(
1998
).
56.
K.-Z.
Shen
,
Z.-T.
Zhu
,
A.
Munhall
, and
S. W.
Johnson
, “
Synaptic plasticity in rat subthalamic nucleus induced by high-frequency stimulation
,”
Synapse
50
,
314
319
(
2003
).
57.
N.
Yamawaki
,
P.
Magill
,
G.
Woodhall
,
S.
Hall
, and
I.
Stanford
, “
Frequency selectivity and dopamine-dependence of plasticity at glutamatergic synapses in the subthalamic nucleus
,”
Neuroscience
203
,
1
11
(
2012
).
58.
L.
Gorodetski
,
R.
Zeira
,
H.
Lavian
, and
A.
Korngreen
, “
Long-term plasticity of glutamatergic input from the subthalamic nucleus to the entopeduncular nucleus
,”
Eur. J. Neurosci.
48
,
2139
2151
(
2018
).
59.
Z.
Li
,
J.-G.
Zhang
,
Y.
Ye
, and
X.
Li
, “
Review on factors affecting targeting accuracy of deep brain stimulation electrode implantation between 2001 and 2015
,”
Stereotact. Funct. Neurosurg.
94
,
351
362
(
2016
).
60.
J.
Voges
,
J.
Volkmann
,
N.
Allert
,
R.
Lehrke
,
A.
Koulousakis
,
H.-J.
Freund
, and
V.
Sturm
, “
Bilateral high-frequency stimulation in the subthalamic nucleus for the treatment of Parkinson disease: Correlation of therapeutic effect with anatomical electrode position
,”
J. Neurosurg.
96
,
269
279
(
2002
).
61.
S. H.
Paek
,
J.-Y.
Lee
,
H.-J.
Kim
,
D.
Kang
,
Y. H.
Lim
,
M. R.
Kim
,
C.
Kim
,
B. S.
Jeon
, and
D. G.
Kim
, “
Electrode position and the clinical outcome after bilateral subthalamic nucleus stimulation
,”
J. Korean Med. Sci.
26
,
1344
1355
(
2011
).
62.
F.
Caire
,
D.
Ranoux
,
D.
Guehl
,
P.
Burbaud
, and
E.
Cuny
, “
A systematic review of studies on anatomical position of electrode contacts used for chronic subthalamic stimulation in Parkinson’s disease
,”
Acta Neurochir.
155
,
1647
1654
(
2013
).
63.
N.
Nakano
,
M.
Taneda
,
A.
Watanabe
, and
A.
Kato
, “
Computed three-dimensional atlas of subthalamic nucleus and its adjacent structures for deep brain stimulation in Parkinson’s disease
,”
ISRN Neurol.
2012
,
592678
(
2012
).
64.
S.
Guo
,
P.
Zhuang
,
M.
Hallett
,
Z.
Zheng
,
Y.
Zhang
,
J.
Li
, and
Y.
Li
, “
Subthalamic deep brain stimulation for Parkinson’s disease: Correlation between locations of oscillatory activity and optimal site of stimulation
,”
Parkinsonism Relat. D
19
,
109
114
(
2013
).
65.
F.
Steigerwald
,
C.
Matthies
, and
J.
Volkmann
, “
Directional deep brain stimulation
,”
Neurotherapeutics
16
,
100
104
(
2019
).
66.
M. F.
Contarino
,
L. J.
Bour
,
R.
Verhagen
,
M. A.
Lourens
,
R. M.
de Bie
,
P.
van den Munckhof
, and
P.
Schuurman
, “
Directional steering: A novel approach to deep brain stimulation
,”
Neurology
83
,
1163
1169
(
2014
).
67.
B.
Lysyansky
,
O. V.
Popovych
, and
P. A.
Tass
, “
Optimal number of stimulation contacts for coordinated reset neuromodulation
,”
Front. Neuroeng.
6
,
5
(
2013
).
68.
M.
Ebert
,
C.
Hauptmann
, and
P. A.
Tass
, “
Coordinated reset stimulation in a large-scale model of the STN-GPE circuit
,”
Front. Comput. Neurosci.
8
,
154
(
2014
).
69.
J. A.
Kromer
and
P. A.
Tass
, “
Long-lasting desynchronization by decoupling stimulation
,”
Phys. Rev. Res.
2
,
033101
(
2020
).
70.
D.
Harnack
,
C.
Winter
,
W.
Meissner
,
T.
Reum
,
A.
Kupsch
, and
R.
Morgenstern
, “
The effects of electrode material, charge density and stimulation duration on the safety of high-frequency stimulation of the subthalamic nucleus in rats
,”
J. Neurosci. Methods
138
,
207
216
(
2004
).
71.
M.
Zeitler
and
P. A.
Tass
, “
Augmented brain function by coordinated reset stimulation with slowly varying sequences
,”
Front. Syst. Neurosci.
9
,
49
(
2015
).
72.
R.
Kempter
,
W.
Gerstner
, and
J.
van Hemmen
, “
Hebbian learning and spiking neurons
,”
Phys. Rev. E
59
,
4498
4514
(
1999
).
73.
A. N.
Burkitt
,
H.
Meffin
, and
D. B.
Grayden
, “
Spike-timing-dependent plasticity: The relationship to rate-based learning for models with weight dynamics determined by a stable fixed point
,”
Neural Comput.
16
,
885
940
(
2004
).
74.
M.
Gilson
,
A.
Burkitt
, and
L. J.
van Hemmen
, “
STDP in recurrent neuronal networks
,”
Front. Comput. Neurosci.
4
,
23
(
2010
).
75.
E. V.
Lubenov
and
A. G.
Siapas
, “
Decoupling through synchrony in neuronal circuits with propagation delays
,”
Neuron
58
,
118
131
(
2008
).
76.
A.
Knoblauch
,
F.
Hauser
,
M.-O.
Gewaltig
,
E.
Körner
, and
G.
Palm
, “
Does spike-timing-dependent synaptic plasticity couple or decouple neurons firing in synchrony?
,”
Front. Comput. Neurosci.
6
,
55
(
2012
).
77.
F.
Steigerwald
,
M.
Potter
,
J.
Herzog
,
M.
Pinsker
,
F.
Kopper
,
H.
Mehdorn
,
G.
Deuschl
, and
J.
Volkmann
, “
Neuronal activity of the human subthalamic nucleus in the parkinsonian and nonparkinsonian state
,”
J. Neurophysiol.
100
,
2515
2524
(
2008
).
78.
P.
Tass
,
D.
Smirnov
,
A.
Karavaev
,
U.
Barnikol
,
T.
Barnikol
,
I.
Adamchic
,
C.
Hauptmann
,
N.
Pawelcyzk
,
M.
Maarouf
,
V.
Sturm
et al., “
The causal relationship between subcortical local field potential oscillations and parkinsonian resting tremor
,”
J. Neural Eng.
7
,
016009
(
2010
).
79.
M. F.
Contarino
,
L. J.
Bour
,
M.
Bot
,
P.
van den Munckhof
,
J. D.
Speelman
,
P. R.
Schuurman
, and
R. M.
de Bie
, “
Tremor-specific neuronal oscillation pattern in dorsal subthalamic nucleus of parkinsonian patients
,”
Brain Stimul.
5
,
305
314
(
2012
).
80.
A. A.
Kühn
,
A.
Kupsch
,
G.-H.
Schneider
, and
P.
Brown
, “
Reduction in subthalamic 8-35 Hz oscillatory activity correlates with clinical improvement in Parkinson’s disease
,”
Eur. J. Neurosci.
23
,
1956
1960
(
2006
).
81.
M.
Weinberger
,
N.
Mahant
,
W. D.
Hutchison
,
A. M.
Lozano
,
E.
Moro
,
M.
Hodaie
,
A. E.
Lang
, and
J. O.
Dostrovsky
, “
Beta oscillatory activity in the subthalamic nucleus and its relation to dopaminergic response in Parkinson’s disease
,”
J. Neurophysiol.
96
,
3248
3256
(
2006
).
82.
D. E.
Feldman
, “
The spike-timing dependence of plasticity
,”
Neuron
75
,
556
571
(
2012
).
83.
R. C.
Froemke
,
M.-M.
Poo
, and
Y.
Dan
, “
Spike-timing-dependent synaptic plasticity depends on dendritic location
,”
Nature
434
,
221
(
2005
).
84.
C.
Pollo
,
A.
Kaelin-Lang
,
M. F.
Oertel
,
L.
Stieglitz
,
E.
Taub
,
P.
Fuhr
,
A. M.
Lozano
,
A.
Raabe
, and
M.
Schüpbach
, “
Directional deep brain stimulation: An intraoperative double-blind pilot study
,”
Brain
137
,
2015
2026
(
2014
).
85.
P. A.
Tass
,
A. N.
Silchenko
,
C.
Hauptmann
,
U. B.
Barnikol
, and
E.-J.
Speckmann
, “
Long-lasting desynchronization in rat hippocampal slice induced by coordinated reset stimulation
,”
Phys. Rev. E
80
,
011902
(
2009
).
86.
J.
Buhlmann
,
L.
Hofmann
,
P. A.
Tass
, and
C.
Hauptmann
, “
Modeling of a segmented electrode for desynchronizing deep brain stimulation
,”
Front. Neuroeng.
4
,
15
(
2011
).
87.
M. D.
Johnson
,
S.
Miocinovic
,
C. C.
McIntyre
, and
J. L.
Vitek
, “
Mechanisms and targets of deep brain stimulation in movement disorders
,”
Neurotherapeutics
5
,
294
308
(
2008
).
88.
V.
Gradinaru
,
M.
Mogri
,
K. R.
Thompson
,
J. M.
Henderson
, and
K.
Deisseroth
, “
Optical deconstruction of parkinsonian neural circuitry
,”
Science
324
,
354
359
(
2009
).
89.
J.-M.
Deniau
,
B.
Degos
,
C.
Bosch
, and
N.
Maurice
, “
Deep brain stimulation mechanisms: Beyond the concept of local functional inhibition
,”
Eur. J. Neurosci.
32
,
1080
1091
(
2010
).
90.
K.
Pyragas
,
V.
Novičenko
, and
P. A.
Tass
, “
Mechanism of suppression of sustained neuronal spiking under high-frequency stimulation
,”
Biol. Cybern.
107
,
669
684
(
2013
).

Supplementary Material