Impact of number of stimulation sites on long-lasting desynchronization effects of coordinated reset stimulation

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,¹ Parkinson’s disease,² and epilepsy.¹⁸ 

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;³ 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.²⁰ Moreover, the limited longevity of batteries in internal pulse generators of DBS devices requires the patient to undergo battery-replacement surgeries.²¹ 

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,²⁴ a single stimulus pulse of moderate amplitude delivered at a vulnerable phase may temporarily desynchronize networks of coupled oscillators.¹⁶ Corresponding desynchronization by repetitive or demand-controlled pulse delivery has been demonstrated computationally^16,27 and in experiments on coupled electrochemical oscillators.²⁸ 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;²⁹ nonetheless, tuning the stimulation frequency out of these parameter regions results in chaotic desynchronization.³⁰ The latter might contribute to the therapeutic effects of deep brain stimulation.³⁰ One can calculate the desynchronizing stimulation frequencies from estimated phase response curves of the synchronous rhythm.³⁰ 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 pulse^28,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).⁴ 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–7 Acute 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.⁸ 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.⁴ 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.¹⁶ 

The initial development of CR stimulation aimed at explicitly using the center manifold theorem⁴² (a dynamical systems theorem related to the slaving principle⁴³ of statistical physics) by delivering CR stimulation to a set of M neuronal subpopulations using, e.g., multipolar deep brain electrodes.⁶ 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,…,M−1$ 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 parameter¹¹ 

$ρ¯Δ$ 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.

Originally, CR stimulation was delivered on demand, such that short stimulation episodes, as shown in Fig. 1, were applied whenever the system resynchronized.⁴ In an alternative embodiment, CR stimuli were delivered periodically with stimulus amplitude adapted to the amount of neuronal synchrony.⁴ 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.⁴⁰ 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;⁴⁰ 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.

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.⁵⁹ 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.⁶⁵ 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.⁶⁶ 

So far, computational studies analyzed the impact of the number of stimulation sites $M$ on acute desynchronization in all-to-all coupled networks.⁶⁷ 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.⁴ 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.⁶⁸ Excitatory synapses between neurons $i$ and $j$ are randomly added. The connection probability $∝exp⁡(dij/0.5mm)$ depends on the distance $dij=|xj−xi|$⁠. In total, 7% of all possible non-autaptic connections are added.⁶⁸ 

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

$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

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

$τsyn$ is the synaptic timescale, $tljj$ the timing of the $lj$th 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 $wj→i(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

with noise conductances $gnoise,i(t)$ resulting from

Here, $κnoise$ scales the noise intensity. $tki$ is the $ki$th 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 $wi→j$ 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 $wi→j→wi→j+W(Δt)$ with time lag $Δt=tpost−tpre−td$ and STDP function

$η≪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 $∫0∞dtW(t)$⁠. In addition, we apply hard bounds $wj→i(t)∈[0,1]$⁠.

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

In neuroscience and medicine, electrical stimuli are typically charge-balanced to avoid tissue damage.⁷⁰ 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,spike−Vreset$⁠. $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 $M≥2$⁠. 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 $m$th 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 $wi→j(t)$⁠. In particular, we randomly set $wi→j(t)$ either to one or zero such that a given mean weight $⟨w(t)⟩=∑ijwi→j(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.

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$⁠).

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 $wi→j(t)$ with presynaptic neuron $i$ and postsynaptic neuron $j$⁠. The dynamics of $wi→j(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

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, $fdom≈3.5Hz$⁠, and the stimulation strength, $Astim$⁠, is close to one. Then, neuronal spike trains can be approximated by

Here, $δ(t)$ is the Dirac delta distribution, $Sl$ are the stimuli $sk$ delivered to neuron $l$⁠, and $ϵk,l$s 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

$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=skj−ski$ 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

For interpopulation synapses, we find

with $Buv=M−u+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

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$⁠).

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 $k$th 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.

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 (⁠$Astim→0$⁠), 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,⁶⁷ 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.⁸² 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 $fCR≈fdom$ yields most pronounced long-lasting desynchronization.⁸ 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.⁴⁶ 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.⁸² For instance, even the location of the synapse receiving input can affect the shape of the STDP function.⁸³ 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.⁵⁷ 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,⁸⁵ in the 1-methyl-4-phenyl-1,2,3,6-tetrahydropyridine Parkinsonian monkey model^5,7 , and in Parkinson’s disease patients,⁶ where up to four stimulation contacts were used. Tass et al.⁸⁵ used an array of four different unipolar stimulation electrodes. As yet, the preclinical and clinical CR Deep Brain Stimulation (CR-DBS) studies^5–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.⁶⁵ 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.⁸⁶ 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.⁶⁷ 

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