Coordinated reset stimulation of plastic neural networks with spatially dependent synaptic connections

Background: Abnormal neuronal synchrony is associated with several neurological disorders, including Parkinson’s disease (PD), essential tremor, dystonia, and epilepsy. Coordinated reset (CR) stimulation was developed computationally to counteract abnormal neuronal synchrony. During CR stimulation, phase-shifted stimuli are delivered to multiple stimulation sites. Computational studies in plastic neural networks reported that CR stimulation drove the networks into an attractor of a stable desynchronized state by down-regulating synaptic connections, which led to long-lasting desynchronization effects that outlasted stimulation. Later, corresponding long-lasting desynchronization and therapeutic effects were found in animal models of PD and PD patients. To date, it is unclear how spatially dependent synaptic connections, as typically observed in the brain, shape CR-induced synaptic downregulation and long-lasting effects.

Methods: We performed numerical simulations of networks of leaky integrate-and-fire neurons with spike-timing-dependent plasticity and spatially dependent synaptic connections to study and further improve acute and long-term responses to CR stimulation.

Results: The characteristic length scale of synaptic connections relative to the distance between stimulation sites plays a key role in CR parameter adjustment. In networks with short synaptic length scales, a substantial synaptic downregulation can be achieved by selecting appropriate stimulus-related parameters, such as the stimulus amplitude and shape, regardless of the employed spatiotemporal pattern of stimulus deliveries. Complex stimulus shapes can induce local connectivity patterns in the vicinity of the stimulation sites. In contrast, in networks with longer synaptic length scales, the spatiotemporal sequence of stimulus deliveries is of major importance for synaptic downregulation. In particular, rapid shuffling of the stimulus sequence is advantageous for synaptic downregulation.

Conclusion: Our results suggest that CR stimulation parameters can be adjusted to synaptic connectivity to further improve the long-lasting effects. Furthermore, shuffling of CR sequences is advantageous for long-lasting desynchronization effects. Our work provides important hypotheses on CR parameter selection for future preclinical and clinical studies.

1 Introduction

The human body consists of various interacting physiological systems forming a complex network in which stimulation of one system can have long-ranging effects on others (Bashan et al., 2012). To develop novel treatments, researchers study how the complex interactions between physiological systems and the various signaling pathways lead to the emergence of pathological and physiological states (Bartsch et al., 2015; Ivanov et al., 2016) and how transitions between these states can be induced, e.g., by stimulation.

Parkinson’s disease is accompanied by synaptic reorganization (Fan et al., 2012; Miguelez et al., 2012; Chu et al., 2017; Pamukcu et al., 2020; Madadi Asl et al., 2022) and changes of neuronal activity (Hammond et al., 2007; Neumann et al., 2016; Vinding et al., 2024) in several brain areas. The severity of motor symptoms, such as bradykinesia and rigidity, correlates with excessive neuronal synchrony (see, e.g., Hammond et al., 2007; Neumann et al., 2016). Brain stimulation has been found to alleviate symptoms, which goes along with a reduction of neuronal synchrony (Hammond et al., 2007). These observations motivated a model of a multistable complex network of PD-related brain areas (Tass and Majtanik, 2006; Tass and Hauptmann, 2007). In this model, pathological states characterized by excessive neuronal synchrony and the presence of symptoms (Hammond et al., 2007) and with reorganized synaptic connectivity (Mallet et al., 2019; Madadi Asl et al., 2022) coexist with healthy states characterized by reduced neuronal synchrony and less severe symptoms or even the absence thereof (Tass and Hauptmann, 2007).

High-frequency (HF) deep brain stimulation (DBS) is an established treatment for medically refractory PD (Krack et al., 2003). However, symptoms return shortly after cessation of stimulation (Temperli et al., 2003). Applying the interpretation of coexisting physiological and pathological states suggests that while HF DBS induces acute therapeutic effects during stimulation, it may not drive the network into the attractor of a stable, healthy state where the network remains after stimulation ceases.

CR stimulation is a multisite stimulation technique computationally developed to specifically counteract pathological neuronal synchrony (Tass, 2003a; Tass, 2003b). During CR stimulation, a spatiotemporal pattern of phase-shifted stimuli is delivered to multiple stimulation sites to target segregated neuronal subpopulations. Originally, CR stimulation was developed to desynchronize neuronal networks with fixed synaptic connections (Tass, 2003a; Tass, 2003b). Later, computational studies of CR stimulation in plastic networks revealed that stimulation not only reduced synchronization during stimulation but also caused downregulated synaptic connections, which drove the network into the attractor of a stable desynchronized state such that desynchronization effects persisted even after cessation of stimulation (Tass and Majtanik, 2006). Later preclinical studies in animal models for PD (Tass et al., 2012; Wang et al., 2016; Bore et al., 2022; Wang et al., 2022) and PD patients (Adamchic et al., 2014; Pfeifer et al., 2021) reported corresponding therapeutic effects that outlasted stimulation for several days to weeks. In contrast to HF DBS, this suggests that CR stimulation indeed drove the network into the attractor of a stable, healthy state, characterized by desynchronized neuronal activity, such that desynchronization effects persisted after cessation of stimulation.

Based on the computational results, the downregulation of (pathological, see Hauptmann and Tass, 2010) synaptic connectivity was key for inducing long-lasting desynchronization effects (Tass and Majtanik, 2006). Since then, several computational studies analyzed the role of different CR parameters for inducing long-lasting effects in plastic networks of phase oscillators (Tass and Majtanik, 2006), conductance-based neuron models (Popovych and Tass, 2012; Ebert et al., 2014; Lourens et al., 2015; Manos et al., 2018), and leaky integrate-and-fire (LIF) neurons (Kromer et al., 2020). These studies pointed out that the spatiotemporal pattern of stimulus deliveries, characterized by the stimulation frequency and the timing and the locations of stimulation site activations, has a strong impact on the CR-induced synaptic dynamics (Kromer et al., 2020; Khaledi-Nasab et al., 2021; Khaledi-Nasab et al., 2022). For CR stimulation, stimuli are delivered in cycles such that each stimulation site receives exactly one stimulus per cycle. The sequence at which sites receive stimuli is called CR sequence. In preclinical and clinical studies on CR, the CR sequence is typically shuffled after each cycle, which will be referred to as shuffled CR below (Tass et al., 2012; Adamchic et al., 2014; Wang et al., 2022). A recent computational and theoretical study also showed that the stimulus shape, characterized by the number of pulses per stimulus burst and the intraburst frequency, has a strong impact on the CR-induced synaptic dynamics (Kromer and Tass, 2022). Beyond the scope of CR stimulation other computational studies analyzed how the pattern of stimulus deliveries can shape plastic synaptic connections (see for instance, Madadi Asl et al., 2023; Anil et al., 2023), how synaptic reshaping might interact with heterogeneous intrinsic neuronal dynamics (Pariz et al., 2023), how spatially dependent connections impact synchronization (Berner and Yanchuk, 2021), and how stimulation-induced synaptic reshaping can lead to lasting effects on neuronal synchrony (Schmalz and Kumar 2019). However, none of these studies analyzed how different spatial patterns of synaptic connectivity might affect the CR-induced synaptic weight dynamics and how such spatial patterns can be harnessed for improved long-lasting effects of CR stimulation. The latter is an important question, as typical target brain regions, for instance, the subthalamic nucleus (STN) for CR DBS in PD, present complex patterns of synaptic connectivity (Emmi et al., 2020). Recently, a brief computational study reported that long-lasting effects of CR in spatially inhomogeneous networks might depend on the selected CR sequence and suggested that rapid shuffling of the CR sequence might be advantageous for long-lasting effects (Kromer and Tass, 2024), indicating that CR might have different effects depending on the synaptic network structure in the target brain area.

The present paper studies CR stimulation in plastic networks with spatially dependent synaptic connections. The network dynamics is modeled by LIF neurons with spike-timing-dependent plasticity (STDP) (Kromer and Tass, 2020). The network parameters are adjusted such that a strongly connected synchronized state and a weakly connected desynchronized state coexist. This was motivated by the presence of pathological states with excessive neuronal synchrony (Hammond et al., 2007) and healthy states characterized by the absence of excessive neuronal synchrony in PD (Tass and Majtanik, 2006). In contrast to previous studies, we vary the spatial structure of synaptic connections and derive stimulation strategies for different network structures. Specifically, we introduce a characteristic length scale of synaptic connections, s, that controls the distance dependence of the connection probability between two neurons. We assumed a decay of the connection probability with the distance between neurons, based on the analysis of synaptic connectivity in several brain regions, including the cortex (Hellwig, 2000), the striatum (Humphries et al., 2010), and as previously assumed in a detailed computational model of the STN (Ebert et al., 2014).

Based on our computational results, we derive the hypothesis that local connections dominate networks with short synaptic length scales and that selecting an appropriate stimulus shape has a significant impact on synaptic downregulation regardless of the spatiotemporal pattern of stimulus deliveries. In contrast, networks with long synaptic connections possess mainly connections between neurons close to different stimulation sites, and CR-induced synaptic downregulation strongly depends on the spatiotemporal pattern of stimulus deliveries.

Our paper is organized as follows. First, we present the neuronal network model and show that different stable states coexist for a wide range of synaptic length scales. In the results section, we present our simulation results on CR stimulation of networks of LIF neurons with STDP and with different synaptic length scales and analyze the stimulation-induced network structure. Then, we derive efficient stimulation strategies for such networks and show how the spatiotemporal stimulus pattern affects the dynamics of synaptic connections. Finally, we discuss our results.

2 Model and methods

2.1 Network of leaky integrate-and-fire neurons

To simulate the effect of CR stimulation on plastic neural networks, we performed simulations of networks of excitatory LIF neurons with STDP. The model equations were taken from previous studies (Kromer et al., 2020; Kromer and Tass, 2020). Here, we extend the model by considering networks with different synaptic length scales as described below.

Individual neurons were modeled as LIF neurons. The ith neuron’s membrane potential, V_i(t), obeyed the dynamics.

$C_i\frac{d{V}_i}{dt}=g_{\text{leak}}\left(V_{\text{rest}}-V_i\right)+g_{\text{syn},i}\left(t\right)\left(V_{\text{syn}}-V_i\right)+I_{\text{stim},i}\left(t\right)+I_{\text{noise},i}\left(t\right).(1)$

Here, C_i is the membrane capacitance and the terms on the right-hand side model the leak current, with leakage conductance g_leak and resting potential V_rest; the excitatory synaptic input current, with synaptic conductance g_syn,i(t) and reversal potential V_syn; the stimulation current I_stim,i(t); and the noisy input current I_noise,i(t), modelling input from other brain regions. Action potentials (spikes) were not modelled directly. Instead neuron i was assumed to exhibit a spike whenever its membrane potential crossed the threshold potential V_th,i(t) from below. V_th,i(t) obeyed the dynamics

${\tau }_{\text{th}}\frac{d{V}_{\text{th},i}}{dt}=-\left(V_{\text{th},i}-V_{\text{th,rest}}\right).(2)$

After threshold crossing, we implemented a rectangular voltage spike by setting V_i(t) to V_spike for a duration of τ_spike. Then, we performed an instantaneous reset of the threshold potential and the membrane potential: V_th,i(t) → V_th,spike and V_i(t) → V_reset.

The synaptic conductances g_syn,i(t) were increased instantaneously whenever a presynaptic spike arrived, and decayed exponentially in between incoming presynaptic spikes:

${\tau }_{\text{syn}}\frac{d{g}_{\text{syn},i}}{dt}=-g_{\text{syn},i}+\kappa \frac{{\tau }_{\text{syn}}}{N}{\displaystyle \sum _{j\in G_i}}{w}_{j\to i}\left(t\right){\displaystyle \sum _{l\in Y_j}}\delta \left(t-t_l^j-t_{\text{d}}\right).(3)$

τ_syn is the synaptic timescale, $t_l^j$ the timing of the lth spike of neuron j, Y_j the set of spikes of neuron j, and t_d the synaptic (axonal) transmission delay. The first sum runs over the set of all presynaptic neurons G_i of neuron i. κ scales the maximum synaptic coupling strength. Additionally, the strength of individual synapses was scaled by the time-dependent synaptic weights w_j→i(t), with j being the index of the presynaptic neuron and i the index of the postsynaptic neuron.

To model noisy input from other brain areas, we delivered independent Poisson spike trains with mean firing rate f_noise to each neuron. These spike trains were fed into the neurons through excitatory synapses. The resulting input current to neuron i is given by

$I_{\text{noise},i}\left(t\right)=g_{\text{noise,i}}\left(t\right)\left(V_{\text{syn}}-V_i\right).(4)$

g_noise,i(t) is the noise conductance given by

${\tau }_{\text{syn}}\frac{d{g}_{\text{noise},i}}{dt}=-g_{\text{noise,i}}+\frac{{\kappa }_{\text{noise}}{\tau }_{\text{syn}}}{N}{\displaystyle \sum _{k\in U_i}}\delta \left(t_k-t\right).(5)$

κ_noise sets the noise intensity and t_k is the timing of the kth spike of the Poisson spike train that was fed into neuron i, and U_i is the set of all spikes in this spike train.

All parameters were taken from Kromer and Tass (2020): g_leak = 0.02 mS/cm², V_rest = −38 mV, V_reset = −67 mV, V_th,spike = 0 mV, V_th,rest = −40 mV, τ_th = 5 ms, V_syn = 0 mV, τ_syn = 1 ms, t_d = 3 ms, κ = 8 mS/cm², κ_noise = 0.026 mS/cm², and f_noise = 20 Hz. The membrane capacitances C_i were Gaussian distributed ($\mathcal{N}({\mu }_{\text{C}},{\sigma }_{\text{C}})$, with mean value μ_C = 3 μF/cm² and standard deviation σ_C = 0.05μ_C). This led to heterogeneity of the intrinsic dynamics of individual neurons. This parameter set was chosen such that the frequency and the range of membrane potential oscillations matched recordings of periodically spiking STN neurons in brain slices of healthy rats (Bevan and Wilson, 1999). The STN is a major target region for HF DBS in PD (Krack et al., 2003) and target for CR DBS in Parkinsonian monkeys (Tass et al., 2012) and PD patients (Adamchic et al., 2014).

2.2 Networks with spatially dependent synaptic connections

To study CR stimulation in networks with spatially dependent synaptic connectivity, we considered networks of N = 1000 LIF neurons. The neurons were placed along the x-axis such that their center coordinates, x_i, were uniformly distributed in the interval 0 ≤ x_i < L. This is a similar distribution as in previous computational studies (Kromer et al., 2020) and was motivated by the distribution of neurons along a DBS electrode. Synaptic connections were introduced randomly between the neurons such that 7% of all N² possible connections are implemented. This was motivated by earlier detailed studies on the effect of CR stimulation on neuronal activity in a detailed model of the STN (Ebert et al., 2014). We did not allow for autaptic connections.

The spatial dependence of the synaptic connections was implemented by considering a distance-dependent connection probability p (d_ij) ∝ exp (−d_ij/s) (Ebert et al., 2014). Here, d_ij = |x_j − x_i| is the Euclidian distance between the presynaptic neuron i and the postsynaptic neuron j. Periodic boundary conditions were not considered. s is the synaptic length scale and will be varied throughout the present paper. Network realizations for the three values of s that were considered are visualized in Figure 1.

Figure 1

Figure 1. Connection probability and connectivity diagrams for network realizations with different synaptic length scales. (A–C): Plots of the connection probability as a function of the distance, d_ij, between two neurons. We normalized p (d_ij) such that ${\int }_0^{\infty }dup(u)=1$. Labels indicate synaptic length scales. (A’–C’) Corresponding connectivity diagrams. Black dots indicate a synaptic connection between a presynaptic neuron at location x_pre and a postsynaptic neuron at location x_post. Panels show results for synaptic length scales s =0.08L (A, A′), 0.4L (B, B′), and 2.0L (C, C′). Here, L is the system’s length scale. In units of the distance between adjacent stimulation sites, d, the synaptic lengths scales are s =0.32 d (A), 1.6 d (B), and 8.0 d (C) (see below).

For s ≪ L the network of synaptic connections is dominated by local connections whereas s ≫ L leads to rather homogeneous synaptic connectivity that becomes independent of the distance between neurons in the limit of long synaptic length scales s/L → ∞. In the present paper, we studied how CR stimulation reshapes the synaptic connectivity for the three synaptic length scales in Figure 1 which are representative of cases where the network is dominated by local connections (s = 0.08 L, Figures 1A, A’), the network shows close-to-homogeneous synaptic connectivity (s = 2 L, Figures 1C, C’), and a case in between these two extreme cases (s = 0.4 L, Figures 1B, B’).

For comparison, we discuss values of L and s used in previous studies modelling the STN. In Kromer et al. (2020), L = 5.0 mm was used, which was motivated by the length of a short axis of an ellipsoidal volume approximation of the STN based on experimental measurements of STN tissue volume (Lévesque and Parent, 2005; Ebert et al., 2014). Magnetic resonance imaging (MRI) studies reported variations of STN dimensions among PD patients with measured mean lengths of 5.9, 3.7, and 5.0 mm along the anteroposterior, mediolateral, and dorsoventral axes, respectively; however, it was noted that MRI yields smaller estimates of STN dimensions than other techniques (Richter et al., 2004). In a detailed computational study, s = 0.5 mm was used to reproduce experimentally measured mean connection lengths (Afsharpour, 1985; Hellwig, 2000; Ebert et al., 2014). Together, these values for L and s would correspond to a synaptic length scale of s ≈ 0.1L, close to the values used in Figures 1A, A’. However, note that the STN is part of the basal ganglia, and its activity is shaped by input from other nuclei (Emmi et al., 2020), making a comparison difficult.

2.3 Spike-timing-dependent plasticity

To model synaptic plasticity, we employed a nearest neighbor STDP scheme (Morrison et al., 2008). The synaptic weights w_i→j were updated whenever a presynaptic or a postsynaptic spike arrived at the synapse. Respective arrival times are given by t_pre + t_d for presynaptic and t_post for postsynaptic spikes. t_pre and t_post are the spike times of the presynaptic and postsynaptic neurons, respectively. At the arrival times, instantaneous synaptic weight updates were performed according to w_i→j → w_i→j + W (Δt) with the time lag Δt = t_post − t_pre − t_d and the STDP function (Song et al., 2000)

$W\left(\mathrm{\Delta }t\right)=\eta \left\{\begin{array}{cc}\exp \left(-\frac{\mathrm{\Delta }t}{{\tau }_{+}}\right),\hfill & \mathrm{\Delta }t>0,\hfill \\ 0,\hfill & \mathrm{\Delta }t=0,\hfill \\ -\frac{\beta }{{\tau }_{\text{R}}}\exp \left(-\frac{|\mathrm{\Delta }t|}{{\tau }_{-}}\right),\hfill & \mathrm{\Delta }t<0.\hfill \end{array}\right.(6)$

η ≪ 1 scales the weight update per spike arrival. 1/η sets the time scale on which STDP affects the synaptic weights. The parameters τ₊ and τ₋ = τ_Rτ₊ determine the shape of the STDP function and are typically in the range of tens of milliseconds (Bi and Poo, 1998). β scales the ratio of overall long-term depression to long-term potentiation, i.e., the integrals over the time lags for which W(Δt) < 0 $(|{\int }_{-\infty }^{0}dtW(t)|)$ and for which W (Δt) > 0 $(|{\int }_0^{\infty }dtW(t)|)$, respectively. STDP is typically considered to be either balanced (β ≈ 1), depression dominated (β > 1), or potentiation dominated (β < 1). In the simulations, we used hard bounds by clipping the synaptic weights after each update to ensure that w_i→j(t) ∈ [0, 1].

Throughout the present paper, we used the parameters β = 1.4, τ₊ = 10 ms, τ_R = 4 (Kromer et al., 2020; Kromer and Tass, 2020). These parameters ensured that desynchronized and synchronized states coexisted for all considered synaptic length scales (Figure 2). We set η = 0.01 which yielded slow changes of the synaptic weights relative to the interspike intervals.

Figure 2

Figure 2. Coexistence of strongly connected synchronized and weakly connected desynchronized states for different synaptic length scales. Simulated traces of the Kuramoto order parameter ρ, Eq. 7, (A–C) for networks with different synaptic length scales, s (columns). Panels in the bottom row show corresponding traces of the mean synaptic weight, ⟨w⟩ obtained by averaging the weights of all synapses in the network. Gray tones mark trajectories for different initial mean synaptic weights. At t = 0 individual synaptic weights were randomly set to either zero or one such that a given mean synaptic weight was realized. Synaptic length scales were the same as in Figure 1, i.e., s = 0.08L (A,A′), s = 0.4L (B,B′), s =2L (C,C′).

2.4 Degree of in-phase synchronization

The degree of in-phase synchronization was quantified by evaluating the Kuramoto order parameter (Kuramoto, 1984)

$\rho \left(t\right)≔\left|\frac{1}{T}{\displaystyle \sum _{i=1}^N}{e}^{2\pi I{\varphi }_i\left(t\right)}\right|.(7)$

ϕ_i(t) is a phase function that attains subsequent integer values at the timings of subsequent spikes of neuron i and increases linearly in between subsequent spikes (Rosenblum et al., 2001). I is the imaginary unit. ρ(t) ≈ 1 indicates in-phase synchronization and ρ(t) ≈ 0 the absence of in-phase synchronization.

2.5 Coordinated reset stimulation

CR stimulation was delivered to four stimulation sites (Tass, 2003b). The sites were placed at locations x_K = L/8, 3L/8, 5L/8, and 7L/8, respectively, such that the distance between two adjacent sites was d = L/4.

Stimulation was delivered in cycles of four stimuli T_CR/4, such that each site received exactly one stimulus per cycle. The cycle period T_CR is given by the inverse CR frequency T_CR = 1/f_CR. Hence, individual sites received stimuli at an average inter-stimulus interval of T_CR. Following, we refer to the sequence of stimulus deliveries during a single CR cycle as CR sequence and to the overall spatio-temporal pattern of stimulus deliveries as CR pattern (Figure 3). We considered two qualitatively different CR patterns: non-shuffled CR and shuffled CR. During non-shuffled CR, one of the 4! possible CR sequences is selected and delivered for the entire stimulation period. In contrast, during shuffled CR, a new CR sequence is randomly selected at the beginning of each CR cycle, i.e., every T_CR (Figure 3). Shuffled CR has been used in preclinical (Tass et al., 2012) and clinical studies (Adamchic et al., 2014) on PD.

Figure 3

Figure 3. Illustration of stimuli and CR pattern-related parameters. We distinguish between individual stimuli and the spatio-temporal pattern of stimulus deliveries (CR pattern). Individual stimuli are characterized by the stimulus shape (left). Two different stimulus shapes are considered throughout the present paper: burst stimuli consisting of three charge-balanced pulses separated by T_intra = 1/f_intra (top left) and single-pulse stimuli (bottom left). Stimulus-related parameters include the stimulus amplitude, the number of pulses per stimulus, the intraburst frequency, f_intra, and the waveform during individual charge-balanced pulses, here characterized by the widths of the excitatory and the inhibitory rectangular pulses (left). Several CR patterns are used throughout the present paper (right). In the right panels, red rectangles mark individual stimuli. We distinguish between non-shuffled CR patterns in which the same CR sequence is repeated for the entire stimulation period (top right) and shuffled CR patterns in which a new CR sequence of stimulus deliveries to the stimulation sites (Roman numerals) is generated after each cycle period, T_CR, by shuffling the M = 4 Roman numerals, referring to individual sites, and randomly drawing M Roman numerals without replacement (bottom right). The shuffled CR pattern shown in the bottom right consists of the CR sequences “I,III,II,IV”, and before the first shuffling, and “II,I,III,IV”, after the first shuffling. After the second shuffling, only the first half of a CR sequence is shown with stimulus deliveries to sites IV and II.

Individual stimuli were modelled as a sequence of charge-balanced pulses. In neuroscience and medicine, electrical stimuli are typically charge-balanced to avoid tissue damage (Harnack et al., 2004). Each charge-balanced pulse consisted of two rectangular current pulses: an excitatory one and an inhibitory one. The excitatory one had a pulse duration of ν_e = 0.4 ms and an amplitude of ${\mathcal{A}}_{\text{e}}=A\mathrm{\Delta }V{\mu }_{\text{C}}/{\nu }_{\text{e}}$. It was immediately followed by the inhibitory one, which had a duration of ν_i = 0.8 ms and an amplitude of ${\mathcal{A}}_{\text{i}}=-{\mathcal{A}}_{\text{e}}{\nu }_{\text{e}}/{\nu }_{\text{i}}$ (Figure 3). ΔV = V_th,spike − V_reset is the maximum voltage range in the subthreshold voltage region of a single LIF neuron in the absence of stimulation and A the dimensionless stimulation strength. A = 1 marks stimuli for which the excitatory pulse resulted in an approximate elevation of the membrane potential of ΔV. Thus, individual stimuli typically elicited a spiking responses.

We employed a spatial stimulus profile to account for the fact that neurons that are further away from a stimulation site, experience weaker stimuli than those close to the site (Richardson et al., 2003). This effectively reduced the stimulus amplitude for neurons at a finite distance |x_i − x_K| from a site located at x_K. Specifically, an LIF neuron with center coordinated x_i experienced stimuli with a reduced amplitude A(x_i, x_K) from a site at x_K (Lysyansky et al., 2013)

$A\left(x_i,x_K\right)=A_{\mathrm{s}\mathrm{t}\mathrm{i}\mathrm{m}}{\left(1+{\left(\frac{x_i-x_K}{\sigma }\right)}^2\right)}^{-1}.(8)$

A(x_i, x_K) is the dimensionless stimulus amplitude experienced by a neuron at location x_i for a stimulus delivered to the stimulation site at location x_K. σ scales the width of the stimulus profile and was set to σ = d/4π throughout the present paper. Thus, during a stimulus delivery to a site at x_K a neuron at location x_i received an excitatory stimulation current, I_stim,i(t), (Eq. 1), with amplitude A(x_i, x_K)ΔVμ_C/ν_e, during the excitatory rectangular pulses and a current with amplitude −A(x_i, x_K)ΔVμ_C/ν_i during the inhibitory pulses.

Below, we study the impact of the stimulus shape, referring to the time course of I_stim,i(t) after stimulus onset, on the outcome of CR stimulation. To this end, we deliver single-pulse stimuli, consisting of one charge-balanced pulse, i.e., one excitatory rectangular pulse followed by an inhibitory one as described above, and burst stimuli consisting of multiple such charge-balanced pulses (Figure 3). Within a burst stimulus, the time lag between subsequent charge-balanced pulses was controlled by the intraburst frequency, f_intra, specifying the inverse time between the onset of two subsequent charge-balanced pulses within a burst.

2.6 Analytical approximation of relative number of synaptic connections between neuronal subpopulations affected by two stimulation sites

In the results section, we present approximations for the mean synaptic weight based on the relative numbers of synaptic connections b_xy≔B_xy/∑_xyB_xy between a presynaptic neuronal subpopulation x and a postsynaptic neuronal subpopulation y. In particular, we consider neuronal subpoplations that are the closest to one of the stimulation sites. Here, B_xy is the total number of connections between the two neuronal subpopulations. The neuronal subpopulation that is the closest to stimulation site K at location x_K = (2K − 1)L/8, K = 1, 2, 3, 4 incorporates all neurons with center coordinates x_i ∈ [(K − 1)L/4, KL/4).

Next, we consider the distribution p_ΔK(d_ij, d) of distances d_ij between postsynaptic neurons j that are the closest to site K_post and presynaptic neurons i that are the closest to stimulation site K_pre with ΔK = K_post − K_pre. For a uniform distribution of neuron center coordinates in each subpopulation, we find.

$p_{\mathrm{\Delta }K}\left(d_{ij},d\right)=\frac{1}{d}\left\{\begin{array}{cc}\left(\frac{d_{ij}}{d}-\left(\mathrm{\Delta }K-1\right)\right),\hfill & d\left(\mathrm{\Delta }K-1\right)\le d_{ij}<d\mathrm{\Delta }K\hfill \\ \left(-\frac{d_{ij}}{d}+\left(\mathrm{\Delta }K+1\right)\right),\hfill & d\mathrm{\Delta }K\le d_{ij}<d\left(\mathrm{\Delta }K+1\right)\hfill \\ 0,\hfill & \text{otherwise}\hfill \end{array}\right..(9)$

For a large number of connections between individual subpopulations, the relative numbers of connections between subpopulations can by approximated by,

$b_{xy}\approx \frac{\underset{-\infty }{\overset{\infty }{\int }}dx{p}_{y-x}\left(dx,d\right)e^{-\frac{|dx|}{s}}}{{\displaystyle \sum _{x,y=1}^M}\underset{-\infty }{\overset{\infty }{\int }}dx{p}_{y-x}\left(dx,d\right)e^{-\frac{|dx|}{s}}}.(10)$

Here, M = 4 is the total number of stimulation sites. This integral can be solved analytically, and the solution depends only on the ratio l = d/s of the distance between adjacent stimulation sites d, characterizing the characteristic length scale of the spatiotemporal stimulus pattern, and the synaptic length scale s. In particular, for l → 0 the result becomes independent of y and x and yield results for homogeneous networks, which have been studied in Kromer and Tass (2022).

2.7 Approximation of mean synaptic weight after CR stimulation in networks with spatially dependent synaptic connections

Below, we explore the effect of the stimulus shape and the CR pattern on the mean synaptic weight after CR stimulation in more detail and develop effective stimulation strategies for synaptic weight reduction in networks with spatially dependent synaptic connections. Recently, we derived theoretical approximations of the stimulation-induced dynamics of the mean synaptic weight for non-shuffled CR stimulation with rectangular stimulation profile in networks with spatially homogeneous synaptic connectivity and STDP (Kromer and Tass, 2022). Given an estimate of the average rate of synaptic weight change during stimulation, J_x→y, of synapses between the presynaptic neuronal population x and postsynaptic neuronal population y and the mean synaptic weight of these synapses at a reference time t₀, ⟨w_x→y (t₀)⟩, their mean synaptic weight at a later time t > t₀ can be approximated by Kromer and Tass (2022)

$⟨w_{x\to y}\left(t\right)⟩\approx {\left[⟨w_{x\to y}\left(t_0\right)⟩+S\left(J_{x\to y},⟨w_{x\to y}\left(t_0\right)⟩\right)\mathrm{t}-{\mathrm{t}}_0\right]}_{\text{clip},\mathrm{0,1}}.(11)$

${\left[a\right]}_{\text{clip},\mathrm{0,1}}=a$ for a ∈ [0, 1], 0 for a < 0 and 1 for a > 1 accounts for the hard bounds for individual synaptic weights. The function

$S\left(\mathcal{J},w\right)\mathrm{t}-{\mathrm{t}}_0≔\left\{\begin{array}{cc}w\mathcal{J}\left(t-t_0\right),\hfill & \mathcal{J}\le 0\hfill \\ \left(1-w\right)\mathcal{J}\left(t-t_0\right),\hfill & \mathcal{J}>0\hfill \end{array}\right.,(12)$

accounts for a linear increase of the weights of initially downregulated synapses and a linear decrease of weights of initially upregulated ones. Note that during the derivation of Eq. 11, it was assumed that a mean synaptic weight of w_x→y(t₀) corresponds to a fraction of w_x→y(t₀) synapses with synaptic weights one and a fraction of 1 − w_x→y(t₀) synapses with synaptic weights zero. This assumption was motivated by the observation that individual synaptic weights are close to the hard bounds in the stationary states (Song et al., 2000; Rubin et al., 2001).

From Eq. 11, an approximation of the mean synaptic weight was derived in Kromer and Tass (2022). To this end, the matrix B was introduced with its entries B_xy being the total numbers of synaptic connections with presynaptic neurons in neuronal subpopulation x and postsynaptic neurons in subpopulaton y. For a given matrix B, the overall mean synaptic weight ⟨w(t)⟩ can be approximated by

$⟨w\left(t\right)⟩\approx \frac{1}{{\sum }_{xy}{B}_{xy}}{\displaystyle \sum _{xy}}{B}_{xy}⟨w_{xy}\left(t\right)⟩={\displaystyle \sum _{xy}}{b}_{xy}⟨w_{xy}\left(t\right)⟩.(13)$

The sums run over all possible combinations of presynaptic neuronal subpopulations x and postsynaptic subpopulations y. b_xy≔B_xy/∑_xyB_xy denotes the relative number of connections with presynaptic neuron in subpopulation x and postsynaptic neuron in subpopulation y.

In the previous section, we derived an analytical estimate for the relative number of connections b_xy for the networks of spatially dependent synaptic connections considered throughout the present paper for the case where a subpopulation incorporates the neurons that are closest to one of the stimulation sites (Eq. 10). Using Eq. 10 in Eq. 13, we calculated the approximate dynamics of the overall mean synaptic weight during non-shuffled CR stimulation with a given CR sequence and a given ratio l of the distance between adjacent stimulation sites d and the synaptic length scale s. In Figure 7, results from Eq. 10 are compared to estimates from networks generated by the algorithm, used in our simulations. Note that Eq. 10 also yields the relative portion of intra-population synapses as b_intra = ∑_xB_xx and inter-population synapses b_inter = 1 − b_intra, i.e., synaptic connections between presynaptic and postsynaptic neurons that are closest to the same and to different sites, respectively.

2.8 Details of numerical simulations and data analyses

Numerical integration of the network model of LIF neurons with STDP was performed using an Euler scheme with integration time step 0.1 ms. The python code used for the simulations and the generation of the figures is available on github.

For simulations of CR stimulation, we used the networks for an initial mean synaptic weight of 0.45 from Figure 2 at t = 5000 s and started stimulation. Thus, the time point at t = 5000 s from Figure 2 corresponds to the onset of stimulation at t = 0 s in Figures 4–6, 8, 9.

Figure 4

Figure 4. CR stimulation-induced dynamics depend on network structure. Simulated traces of the Kuramoto order parameter ρ, Eq. 7, for networks with different synaptic length scales. s = 0.08 L (A,D,G), 0.4 L (B,E,H), and 2.0 L (C,F,I), and for different stimulation frequencies, f_CR = 4 Hz (A–C), 10 Hz (D–F), and 21 Hz (G–I). Individual stimuli consisted of bursts of three stimulus pulses and an intraburst frequency of 130 Hz, the corresponding time between subsequent pulses within a burst, T_intra, is 1/130 s (see schematic on the top right). Black traces show results for shuffled CR, often referred to as CR with rapidly varying sequence (Tass and Majtanik, 2006; Zeitler and Tass, 2015)). The stimulation period of 1000 s is marked in light red. Colored traces show results for non-shuffled CR with CR sequences “I,II,III,IV”,“I,II,IV,III”,“I,III,II,IV”,“I,III,IV,II”,“I,IV,II,III”, and “I,IV,III,II”. The corresponding sequences of stimulation site activations within a CR cycle are illustrated on the right-hand side. Parameters: A_stim = 2.5 and σ = d_s/4 Π.

Figure 5

Figure 5. Synaptic pathways induced by CR stimulation with burst stimuli. Results of simulations of the LIF network model prior to stimulation (A–D), and during stimulation with non-shuffled CR with CR sequences “I,II,III,IV” (E–H) and “I,IV,II,III” (I–L) and shuffled CR (M–P) are shown. Left panels show raster plots of neuronal spiking activity. Black horizontal bars mark a time interval of 100 ms. The other panels show connectivity diagrams with dark colors marking strong connections and light gray marking weak connections between presynaptic neurons at locations x_pre and postsynaptic neurons at locations x_post. Different columns show results for networks with different synaptic length scales, s (see labels at the top of respective columns). In the brackets, we give synaptic length scales in units of the distance between adjacent stimulation sites, d. The labels inside individual panels show the mean synaptic weight ⟨w⟩. Parameters: A = 2.5, f_CR = 10 Hz and σ = d/4 Π. Raster plots in panels (A,E,I, M) are shown for s = 0.08L. Connectivity diagrams (F–H,J–L, N–P) were recorded at the end of the 1000 s stimulation period. Raster plots (E,I,M) show the last second of this period.

Figure 6

Figure 6. Synaptic pathways induced by CR stimulation with single-pulse stimuli. Results of stimulation of the LIF network model prior to stimulation (A–D), and during simulation with non-shuffled CR with CR sequences “I,II,III,IV” (E–H) and “I,IV,II,III” (I–L) and shuffled CR (M–P) are shown. In contrast to Figure 5, we used single-pulse stimuli. The raster plots on the left show neuronal spiking activity. Black horizontal bars mark a time interval of 100 ms. The other panels show connectivity diagrams with dark colors marking strong connections and light gray marking weak connections between presynaptic neurons at locations x_pre and postsynaptic neurons at locations x_post. Different columns show results for networks with different synaptic length scales, s (see labels at the top of respective columns). Labels inside the panels show the value of the mean synaptic weight. Parameters: A = 2.5 and σ = d/4 Π. Raster plots in panels (A,E,I, M) are shown for s = 0.08L. Connectivity diagrams were recorded at the end of the 1000 s stimulation period. The raster plots on the left show the last second before stimulation onset (A) and the last second of the 1000 s of CR stimulation period (E,I,M).

The relative numbers of connections b_xy in Figure 7 were obtained by averaging over five network realizations. The same networks were used in the simulations for Figure 8. Analytical approximations in Figure 7 and for the approximations in Figure 9 were obtained from Eq. 10.

Figure 7

Figure 7. Relative numbers of connections between neuronal subpopulations closest to different combinations of stimulation sites for different synaptic length scales (A–C). Black horizontal lines represents analytical estimates (Eq. 10) and crosses show data from five network realizations. Presynaptic and postsynaptic neuronal subpopulations are denoted by the Roman numeral of the closest stimulation site, i.e., “I-II” refers to the presynaptic subpopulation of all neurons within distance d/2 of stimulation site “I” and postsynaptic subpopulation of all neurons within distance d/2 of stimulation site “II”.

Figure 8

Figure 8. Stimulation-induced change of mean synaptic weight depends on phase lag between stimuli delivered to postsynaptic and presynaptic neuronal subpopulation. (A): Raster plot of neuronal spiking activity during two-site stimulation with single-pulse stimuli for ϕ_x→y = 0.3. Phase lags were normalized to one, i.e., ϕ_x→y ∈ [0, 1), such that ϕ_x→y/f is the time lag between the onsets of stimuli delivered to population y and x, respectively. (A)’: Change of mean synaptic weight of synapses between postsynaptic neurons closest to a site located at 5L/8 and presynaptic neurons closest to a site located at 3L/8 in the first 20 s of stimulation (enclosed by red dashed lines in (A) as function of ϕ_x→y and stimulation frequency, f. (B) and (B)’: Same as (A) and (A′) but for stimulation with burst stimuli with three pulses per burst and an intraburst frequency of 130 Hz. (C, D): Examples of CR sequences (top). Roman numerals denote stimulation sites. Bottom panels show corresponding matrices of resulting phase lags between stimuli delivered to postsynaptic and presynaptic neurons affected by respective stimulation sites.

Figure 9

Figure 9. Mean synaptic weights of intra- (A–C) and (D–F) and inter-population synapses (A’–C’) and (D’–F’) after long stimulation with non-shuffled CR with different CR sequences and shuffled CR. Approximations from Eq. 14 (lines) for t = 3000 s are compared to simulation results after 3000 s of stimulation (symbols) for networks with different synaptic length scales, s, (columns) and for single-pulse (A–C) and (A’–C’) and burst stimuli (D–F) and (D'-F'). Colors mark different CR patterns (see legend in panel (A) and Figure 4). Parameters: A_stim = 2.5, σ = d/4 Π.

Estimates of the mean weight change during stimulation, Δw, in Figures 8A’, 8B’, were obtained by averaging the weights of all synapses between presynaptic neurons within a distance of less than d/2 from the stimulation site at x_II = 3L/8 and postsynaptic neurons within a distance of less than d/2 from the stimulation site at x_III = 5L/8, at the onset of stimulation and after 20 s of stimulation with a phase lag ϕ_x→y and a stimulation frequency f. Then the corresponding mean synaptic weights were averaged over five network realizations.

3 Results

We performed simulations of non-shuffled and shuffled CR stimulation of networks of LIF neurons with STDP. Before stimulation, the networks were in a strongly connected, synchronized state. Specifically, we used the networks obtained for an initial mean synaptic weight of 0.45 from Figure 2 after a simulation time of 5000 s. We delivered CR stimulation with the primary goal to downregulate synaptic weights and drive the network into the attractor of a stable, weakly connected desynchronized state (see Figure 2) such that stimulation entailed long-lasting desynchronization effects (see Figure 4). We recorded traces of the mean synaptic weight ⟨w(t)⟩ and the Kuramoto order parameter, ρ(t), quantifying the degree of neuronal synchrony (Eq. 7).

3.1 CR stimulation-induced dynamics depends on network structure

Representative trajectories of the Kuramoto order parameter, ρ(t), and the mean synaptic weight, ⟨w(t)⟩, before, during, and after CR stimulation with burst stimuli are shown in Figure 4. Simulations were performed for networks with different synaptic length scales. CR stimulation induced complex dynamics of the Kuramoto order parameter and the mean synaptic weight (Figure 4). While various CR patterns led to acute desynchronization (reflected by low values of Kuramoto order parameter during stimulation), the degree of long-term synchronization varied across networks with different synaptic length scales and stimulation parameters. In networks with short synaptic length scales, s, non-shuffled CR with different CR sequences led to different values of the mean synaptic weight during CR stimulation. In contrast, non-shuffled CR stimulation of networks with long synaptic length scales led to similar mean synaptic weights during stimulation (compare first and last column in Figure 4).

After cessation of stimulation, the networks freely evolved and approached either the synchronized, a desynchronized, or a partially synchronized state. In the latter, neurons in the immediate vicinity of stimulation sites remained synchronized whereas neurons further away from stimulation sites showed desynchronized activity (see raster plots of spiking activity and snapshots of synaptic connectivity for simulations of shuffled CR with high CR frequencies in Supplementary Figure S1). In networks with short synaptic length scales, it took a long time for the networks to approach a stationary state with transient dynamics occurring on time scales that were even longer than the stimulation period (Figure 4, left column). In some cases, non-shuffled CR with different sequences also led to different long-term effects of stimulation (Figure 4F). Surprisingly, in these cases, different stationary states were approached even though similar mean synaptic weights were obtained at the end of the stimulation period, indicating a complex attractor landscape of the stationary states.

For most parameter sets used in Figure 4 and sufficiently long stimulation periods, shuffled CR led to a more pronounced reduction of the mean synaptic weight during stimulation and lower values of the Kuramoto order parameter after cessation of stimulation than non-shuffled CR. For large CR frequencies, only shuffled CR kept the system from returning to the synchronized state. Here, the system approached a partially synchronized state that persisted even for long simulation times (we simulated up to 10000 s after cessation of stimulation) (Figures 4G–I). However, at the beginning of the stimulation period, non-shuffled CR led to a faster reduction of the mean synaptic weight than shuffled CR for a low stimulation frequency (see Figures 4A–C). It is currently unclear whether this effect is strong enough to enable long-lasting desynchronization effects after non-shuffled CR for shorter stimulation periods than required for long-lasting desynchronization after shuffled CR.

3.2 CR pattern determined spatial orientation of stimulation-induced synaptic pathways

Next, we analyzed the effect of non-shuffled CR on the synaptic weight dynamics. Snapshots of connectivity diagrams during CR stimulation are shown in Figure 5. Non-shuffled CR upregulated certain populations of synaptic connections while down-regulating others, thereby inducing strongly connected clusters of neurons. To analyze which synapses were affected, we drew red dashed lines in the connectivity diagrams in Figure 5 at coordinates L/4, L/2, and 3L/4, i.e., marking halfway points between adjacent stimulation sites. Vertical lines separate presynaptic neuron coordinates into four groups, each containing presynaptic neurons that were the closest to the same stimulation site. Accordingly, horizontal lines separate postsynaptic neuron coordinates into four groups containing postsynaptic neurons that were the closest to the same stimulation site. As can be seen in Figure 5, these lines approximately separate synaptic connections into different populations with presynaptic and postsynaptic neurons closest to the same stimulation sites.

Remarkably, non-shuffled CR with different CR sequences up- and downregulated different synaptic populations, thereby shaping the orientation of groups of strong synapses marking pathways of strong synaptic connections in the network. For example, in Figure 5H, the CR sequence I,II,III, IV induced the following pathways of strong synaptic connections I → II, II → III, III → IV, IV → I. Here, X → Y, denotes strong synaptic connections between presynaptic neurons with closest site X and postsynaptic neurons with closest site Y. In contrast, the CR sequence I, IV, II, III induced the pathways I → IV, II → III, III → I, and VI → II of strong synapses (Figure 5L).

In networks with short synaptic length scales, some of these pathways only contained a small number of synapses (compare Figures 5F–H; Figures 5J–L). This resulted in a different structure of CR stimulation-induced synaptic pathways due to the lag of synaptic connections between certain neuronal subpopulations. Consequently, different CR sequences led to different mean synaptic weights after non-shuffled CR stimulation of networks with short synaptic length scales if up- and downregulated blocks contained different numbers of synaptic connections.

Comparing the results for non-shuffled CR using burst stimuli (Figure 5) as opposed to single-pulse stimuli (Figure 6), we found that both induced similar patterns of strongly connected synaptic pathways. Noteworthy, for both types of stimuli, no such pathways were induced by shuffled CR variants (Figures 5, 6, bottom row). This was robust with respect to different realizations of random CR sequences during shuffled CR. These results suggests that the structure of CR stimulation-induced synaptic pathways is mainly determined by the CR pattern and less affected by the employed type of stimulus.

3.3 Stimulus shape determines network structure in the vicinity of stimulation sites

Next, we studied how stimuli affected the local network structure in the vicinity of the stimulation sites. In the connectivity diagrams in Figures 5, 6, the local network structure consists of the populations of synapses that are located along the diagonals. We will refer to these synapses as intra-population synapses as they connect presynaptic and postsynaptic neurons that are close to the same stimulation site, as opposed to inter-population synapses where presynaptic and postsynaptic neurons are close to different sites. We thereby follow the terminology of previous studies in which a rectangular spatial stimulus profile was implemented (i.e., all neurons within a certain distance of a stimulation site received the same stimulation current (Kromer et al., 2020; Kromer and Tass, 2022)). Of note, in the present study, the stimulation current decays with increasing distance from the stimulation site (Equation 8).

Stimulation-induced reshaping of intra-population synapses was strongly affected by the employed stimulus shape. Delivering CR stimulation with burst stimuli and a spatial stimulus profile (Eq. 8) caused a reshaping of local synaptic connectivity. Specifically, synaptic connections with postsynaptic neurons that were close to the stimulation site were upregulated. In contrast, connections with postsynaptic neurons that were further away from their closest stimulation site were downregulated. This resulted in the horizontal black regions along the diagonals in the connectivity diagrams during CR stimulation with burst stimuli in Figure 5. This complex local network structure was in marked contrast to CR stimulation with single-pulse stimuli where all intra-population synapses were downregulated (Figure 6). Consequently, in networks with short synaptic length scales, a substantial reduction of the mean synaptic weight can be achieved if the “right” stimulus is used, i.e., single pulse instead of burst stimuli, in our model. Note that inter-population synapses were affected in a similar way for both burst and single-pulse stimuli (Figures 5, 6), as discussed in the previous paragraph.

3.4 Effect of stimulus shape and CR pattern on synaptic weights

The results in Figure 7 show that networks with short synaptic length scales (Figure 7A) contain mostly intra-population synapses; therefore, for effective synaptic downregulation, the focus should be on reducing the strength of intra-population synapses. In contrast, in networks with long synaptic length scales (Figure 7C), the distribution of synapses across the different combinations of presynaptic and postsynaptic neuronal subpopulation is close to uniform. Thus, the major portion of synapses are inter-population synapses, and the focus should be on reducing the strength of such synapses when aiming for synaptic downregulation.

Combining the pronounced difference between the distributions of synapses in networks with different synaptic length scales and the results on the different effects of the CR pattern and the stimulus shape on the connectivity diagrams in Figures 5, 6 from the previous sections, we formulate the hypothesis, that selecting an effective stimulus shape can substantially improve stimulation-induced synaptic weight reduction in networks with short synaptic length scales; whereas, in networks with long synaptic length scales the selection of the CR pattern strongly affects synaptic weight reduction. Here, we use the term CR pattern to refer to the spatiotemporal characteristics of stimulus deliveries, characterized by the CR sequence, the presence or absence of shuffling as well as the CR frequency. In contrast, the stimulus shape is characterized by the stimulation amplitude, the number of stimulus pulses per burst, and the intra-burst frequency (Figure 3).

Accordingly, comparing the results for networks with short synaptic length scales in Figures 5, 6, we see that the difference between mean synaptic weights after CR stimulation with single-pulse (Figure 6) and burst stimuli (Figure 5) is substantial for all three CR patterns, i.e., for non-shuffled CR with sequence I,II,III, IV compare Figures 5F, 6F and with sequence I,IV,II, III compare Figures 5J, 6J, and for shuffled CR compare Figures 5N, 6N. In contrast, for long synaptic length scales the difference between connectivity diagrams and mean synaptic weights for the same CR pattern for burst and single pulse stimuli is smaller (compare Figures 5H, 6H, 5L, 6L; Figures 5P, 6P, respectively).

3.5 CR pattern selection for STDP-induced synaptic weight reduction

Finally, we studied how synaptic weight reduction depended on the selected CR pattern. We focused on the case where neuronal responses to stimuli were short compared to the inter-stimulus intervals.

The dynamics of the mean synaptic weight for a single synaptic population ⟨w_x→y⟩ during CR stimulation strongly depends on the mean rate of weight change (Eq. 11 for non-shuffled CR). Note that for shuffled CR, the dynamics is more complex and was studied for single-pulse stimuli in Kromer et al. (2020).

For non-shuffled CR, an estimate of the function S (Eq. 11) for given CR parameters can be obtained by stimulating two subpopulations with stimulation frequency f for a time period T (Kromer and Tass, 2022) and evaluating the change of the mean synaptic weight of synapses interconnecting the two subpopulations Δw = ⟨w_x→y(t₀ + T)⟩ − ⟨w_x→y(t₀)⟩. This setup is illustrated in Figures 8A,B for single-pulse and burst stimuli, respectively.

In Figures 8A’, 8B’, we show Δw for single-pulse stimuli (Figures 8A, A’) and burst stimuli (Figures 8B, B’) with three pulses per burst for different stimulation frequencies f and phase lags, ϕ_x→y, between stimuli delivered to the two subpopulations. We used short stimulation periods of T = 20 s to minimize the effect of the hard bounds on individual synaptic weights. Note that for long stimulation times T → ∞, synaptic weights approach the hard bounds, and Δw either approaches −⟨w(t₀)⟩ indicating a reduction of individual weights towards the lower hard bound, or 1 − ⟨w(t₀)⟩ indicating an increase of individual weights towards the upper hard bound.

In a wide range of stimulation frequencies, the sign of Δw depended on the phase lag ϕ_x→y, indicating that the adjustment of the phase lag between stimuli determines whether synapses between these neuronal subpopulations strengthen (Δw > 0) or weaken (Δw < 0).

In the schematics in Figures 8C, C’, and Figures 8D, D’, we show the phase lags between stimuli delivered to different stimulation sites for two CR sequences. For M = 4 stimulation sites, each CR sequence induces M times the phase lags 0, 1/M, 2/M, 3/M, … , (M − 1)/M, between neuronal subpopulations that are closest to the M stimulation sites; however, the pairs of neuronal subpopulations that receive stimuli at the phase lags vary across CR sequences (Figures 8C’, D’). As a consequence, the resulting change of the mean synaptic weight (Eq. 13) may depend on the CR sequence in inhomogeneous networks $b_{x_1{y}_1}\ne b_{x_2{y}_2}$. For a given CR sequence, the mean synaptic weight at time t > t₀ during stimulation can be approximated by

${\left.⟨w\left(t\right)⟩\right|}_{\mathit{\varphi }}\approx {\displaystyle \sum _{xy}}{b}_{xy}{\left[⟨w_{x\to y}\left(t_0\right)⟩+S\left(J\left({\varphi }_{x\to y},f_{\text{CR}}\right),⟨w_{x\to y}\left(t_0\right)⟩\right)\mathrm{t}-{\mathrm{t}}_0\right]}_{\text{clip},\mathrm{0,1}}.(14)$

ϕ denotes the matrix of the M² phase lags ϕ_x→y between the different combinations of postsynaptic subpopulations, y, and presynaptic subpopulations, x (see examples in Figures 8C’, D’).

Approximating $\mathcal{S}(J({\varphi }_{x\to y},f_{\text{CR}}),⟨w_{x\to y}(t_0)⟩)\mathrm{t}-{\mathrm{t}}_0\approx \mathrm{\Delta }w({\varphi }_{x\to y},f=f_{\text{CR}})(\mathrm{t}-{\mathrm{t}}_0)/\mathrm{T}$ and b_xy using Eq. 10, we evaluated ${\left.⟨w(t)⟩\right|}_{\mathit{\varphi }}$ for the matrices ϕ for different CR sequences and networks with different synaptic length scales. Note that for long stimulation times t − t₀ → ∞ the sign of Δw is sufficient as the mean weight of synapses between subpopulations will either approach zero or one. The results are compared to simulations of the LIF network in Figure 9.

As can be seen in Figure 9, the approximations based on Eq. 14 well approximated the mean weight of inter-population synapses for single-pulse and burst stimuli (Figures 9A’–F’) and for intra-population synapses when single-pulse stimuli were used (Figures 9A–C). For burst stimuli, large deviations between approximations and simulation results occurred. These deviations were caused by the local network structure (see Figure 5), which was not considered by the approximations. Instead, approximations were solely derived based on the mean weight change for the synaptic populations and assumed that weights of all synapses within the population approach either one or zero. Furthermore, for low CR frequencies, neurons have enough time between stimuli to respond to input from other subpopulations, which causes deviations from the two-site stimulation results (Figure 8) used to calculate the approximations.

Overall, the stimulation-induced mean synaptic weight of intra-population synapses was more robust with respect to changes in the CR frequency than the one of inter-population synapses.

4 Discussion

CR DBS induced long-lasting therapeutic effects in Parkinsonian monkeys (Tass et al., 2012; Wang et al., 2016; Bore et al., 2022; Wang et al., 2022) and PD patients (Adamchic et al., 2014). Therapeutic long-lasting effects were also observed in PD patients who received non-invasive vibrotactile CR finger-tip stimulation (Pfeifer et al., 2021). Computational studies in plastic neural networks and networks of phase oscillators predicted that CR stimulation may have long-lasting effects resulting from a downregulation of plastic synaptic connections, in this way driving the network from a pathological state, prior to stimulation, into the basin of attraction of a stable physiological state in which the network remains after cessation of stimulation (Tass and Majtanik, 2006). In previous computational studies, CR stimulation was studied in different network models with plasticity, including networks of phase oscillators (Tass and Majtanik, 2006) and spiking neuronal networks (Popovych and Tass, 2012; Manos et al., 2018; Kromer et al., 2020). There, the focus was on the influence of stimulation parameters such as the CR frequency and the stimulus amplitude on synaptic downregulation and long-lasting desynchronization effects. Here, we focused on networks with spatially dependent synaptic connections and derived stimulation strategies for CR stimulation that harness spatially inhomogeneous network structure for efficient synaptic downregulation during stimulation.

Multiple stable states can coexist in plastic networks (Seliger et al., 2002; Zanette and Mikhailov, 2004; Maistrenko et al., 2007; Madadi Asl et al., 2016; Madadi Asl et al., 2018a; Madadi Asl et al., 2018b). The networks studied here showed coexisting synchronized and desynchronized stationary states (Figure 2) and partially synchronized states in which the system remained for several hours after cessation of shuffled CR stimulation with burst stimuli (Figure 4 and Supplementary Figure S1). During CR stimulation, we observed rapid changes of the synaptic connectivity that were followed by sometimes long and complex transient dynamics after cessation of stimulation before the network eventually approached either stationary state (Figure 4). These transient dynamics included short periods of rebound synchronization (Figures 4D–I) (Tass and Hauptmann, 2006) and long transient desynchronization that was sometimes followed by a relaxation back to a stationary synchronized state (Figure 4F). These results suggest that in a clinical setting therapeutic CR stimulation may be followed by epochs of transient dynamics, that are potentially longer than the actual stimulation period.

Figure 4 also tells us that pronounced acute desynchronization by CR stimulation does not necessarily imply long-term desynchronization: For instance, while shuffled CR as well as all types of non-shuffled CR induced a pronounced acute desynchronization for f_CR = 21 Hz, shuffled CR entailed partially synchronized network dynamics and the system resynchronized after non-shuffled CR stimulation (see Figures 4G–I). Hence, in this computational model, acute desynchronization, measured during stimulus delivery, does not provide sufficient information for calibration and/or prediction of long-term outcome of CR stimulation. In clinical applications, this may make it difficult to predict long-term therapeutic outcome based on solely measuring acute effects during stimulus delivery.

Targeting synaptic plasticity for long-lasting therapeutic effects remains challenging as stimulation solely provides control over the synaptic weight dynamics during stimulation. However, even networks with similar mean synaptic weights at the end of the stimulation period may approach different stationary states after cessation of stimulation (see, e.g., Figure 4F). Furthermore, it is currently not possible to assess synaptic weights directly during stimulation. In the present study, we implemented a procedure suggested by Kromer and Tass (2022) during which an estimate of the synaptic weight dynamics induced by the employed stimulation is obtained during a two-site stimulation setup (Figure 8), and then this estimate is used to predict which CR patterns lead to the desired synaptic weight reshaping (Figure 9). However, it is unclear how direct estimates of synaptic weight changes can be obtained in patients as classic STDP experiments, e.g., as performed in Bi and Poo (1998), are not feasible. For such a procedure, one may rely on results from corresponding animal experiments or indirect measures, such as changes of evoked responses following stimulation as commonly used in the context of plasticity induced by transcranial magnetic stimulation (Kozyrev et al., 2018; Schwab et al., 2019).

Our results on CR stimulation of neuronal networks with spatially dependent synaptic connections and STDP revealed that fundamental characteristics of such networks can be harnessed to adjust CR stimulation for efficient synaptic weight downregulation. Considering CR stimulation of networks with different synaptic length scales, our results suggest that the ratio between the characteristic distance between stimulation sites, d, and the characteristic synaptic length scale, s, determines whether a network is dominated by intra-population connections (s ≪ d) or inter-population synapses (s ≫ d). In the former case, the shape of the employed stimulus, characterized by parameters such as the number of pulses per stimulus burst, the intra-burst frequency, and the stimulus amplitude, is critical for efficient synaptic weight reduction. Our computational results suggested that such a stimulus could be combined with the shuffled CR pattern (often referred to as CR with rapidly varying sequence, Tass and Majtanik (2006)) for effective synaptic downregulation. We find that shuffled CR is capable of down-regulating synaptic connections between separately stimulated neuronal subpopulation, whereas some of these connections remain strong after non-shuffled CR (Figures 5, 6). The results in Figure 4 further suggest that long-lasting desynchronization effects of shuffled CR are more robust with respect to changes of the CR frequency and across networks with different synaptic length scales than those of non-shuffled CR. This is in accordance with previous computational (Kromer and Tass, 2024) and experimental results in Parkinsonian monkeys that reported that shuffled CR outperformed non-shuffled CR in terms of long-lasting therapeutic effects (Wang et al., 2022). We speculate that in some target brain regions, it may not be necessary to stimulate according to a CR pattern to induce long-lasting therapeutic effects as long as appropriate stimuli are employed. In this context, we note that periodic delivery of burst stimuli to the external segment of the globus pallidus, which is strongly interconnected with the STN, of dopamine-depleted mice induced long-lasting therapeutic effects that depended on the stimulus triggering certain neuronal responses (Spix et al., 2021). In contrast, if networks consist of mostly inter-population synapses, our results suggest that the employed CR pattern is of great importance, and that the choice of parameters characterizing the CR pattern, such as the CR frequency and the CR sequences and their shuffling, has a strong impact on synaptic weight reduction, as suggested by our computational results (Figures 5, 6). So far, clinical studies in PD patients used CR frequencies of 3–20 Hz that were adjusted to the individual patients’ dominant peaks in the LFP power spectrum between 2 and 35 Hz (Adamchic et al., 2014). An adjustment of the CR frequency to the frequency of the dominant synchronous rhythm was suggested by the original computational studies on CR in networks of phase oscillators (Tass, 2003b). Later, computational studies in plastic networks found a non-trivial dependence of long-lasting effects on the CR frequency and observed substantially weaker long-lasting effects for certain unfavorable frequencies (Manos et al., 2018; Kromer and Tass, 2020). Another study found that the interplay of synaptic transmission delays, plasticity parameters, and stimulation determines ranges of unfavorable CR frequencies. Robust long-lasting desynchronization was observed for low CR frequencies (approximately $\le 20$ Hz for four stimulation sites) (Kromer et al., 2020). As details on individual plasticity parameters in the STN are currently not known (see, however, Shen et al. (2003)), it remains to be seen to which extend the CR frequency impacts long-lasting effects in STN CR DBS.

We further found that carefully selected stimuli have the potential to induce non-trivial local network structure in the vicinity of the stimulation site (Figure 5). In particular, CR stimulation employing burst stimuli resulted in the strengthening of local synaptic connections with postsynaptic neurons that were very close to the site and a weakening of the connections with postsynaptic neurons that were further away from the site (Figure 5). This is a consequence of a complex interplay of synaptic plasticity and the distance dependence of stimulus-induced spiking responses due to the spatial stimulus profile (Eq. 8). Specifically, neurons close to the center of the stimulus profile experienced stronger stimulation (Eq. 8) that was able to induce one neuronal spike per stimulus pulse. In contrast, neurons further away from the center responded only reliably to the first stimulus pulse in each stimulus burst as they experienced a weaker stimulation that was not able to cause a spiking response shortly after a neuron’s spike. As a consequence, several postsynaptic spikes followed one presynaptic spike, leading to the upregulation of synapses with postsynaptic neurons close to the center of the stimulus profile (Figure 5). In marked contrast, single-pulse stimuli caused a single synchronized spiking response of high fidelity that is known to lead to synaptic downregulation in networks with axonal transmission delays that are longer than the width of the synchronized spiking response (Lubenov and Siapas, 2008; Knoblauch et al., 2012). Here, neurons were distributed along a line and the stimulus profile depended solely on the distance to the stimulation site. It would be interesting to study the stimulus-induced local network structures in higher-dimensional distributions of neurons with stimulus profiles without spherical symmetry, which could be realized by directional steering in DBS (Contarino et al., 2014) or as a consequence of anisotropic tissue in the target brain area (McIntyre et al., 2004).

It remains to be shown whether CR stimulation with single pulses (Figure 6) as opposed to bursts (Figure 5) is actually feasible in pre-clinical and clinical applications. So far, in pre-clinical and clinical CR studies bursts were used instead of single pulses (Tass et al., 2012; Adamchic et al., 2014; Wang et al., 2016; Bore et al., 2022; Wang et al., 2022). Based on the results obtained above, one may try to deliver single pulse CR as it might be more effective in inducing synaptic downregulation. The decoupling properties of single-pulse stimulation were also observed in hippocampus (Lubenov and Siapas, 2008); however, it is unclear whether this also applies to target regions for DBS in PD, e.g., the STN. We also want to note that qualitatively different model neurons as well as biological neurons may require stronger stimuli. In this case, bursts may provide stimuli of sufficient strength obeying tissue safety limits, as, e.g., demonstrated for tremor entrainment computationally as well as in a patient with spinocerebellar ataxia type 2 (SCA2) with periodic stimuli delivered through a subthalamic-thalamic DBS electrode (Barnikol et al., 2008). For the computational part of that study a network of FitzHugh-Nagumo neurons (FitzHugh, 1961; Nagumo et al., 1962) was used.

Our results suggest that the shuffling of CR sequences substantially improves the induced downregulation of synaptic weights (Figures 4, 9). While the outcomes of non-shuffled CR with different CR sequences differed in some networks, shuffled CR consistently led to either similar or better long-lasting desynchronization effects then non-shuffled CR with the best-performing sequence. It is important to note that this effect may be due to the considered spatial dependence of synaptic connections. In particular, the probability for synaptic connections between two neurons was decaying exponentially in our networks, which was motivated by previous works (Hellwig, 2000; Ebert et al., 2014). However, the resulting connectivity diagrams are approximately symmetric (Figure 1), and the networks do not realize a preferred direction of synaptic connections. The differences between CR sequences in networks with preferred direction of synaptic connections may be larger (see Kromer and Tass, 2024). For instance, based on Figures 8A’, D’, we would expect that the CR stimulation-induced mean synaptic weight for the CR sequence shown in Figure 8D would approach zero for long stimulation times if a network had only connections between neurons close to the following sites: I → III, II → IV, III → II, and IV → I, the first Roman numeral indicating the site closest to presynaptic neurons and the second one the site closest to postsynaptic neurons.

It will be important to test our results and predictions in detailed models of the basal ganglia region, in particular, the STN and globus pallidus. STN neurons are subject to several inputs from other nuclei, e.g., gamma-aminobutyric acidergic inputs from neurons in the external segment of the globus pallidus. In our model, we only considered excitatory neurons. Furthermore, synaptic inputs to the STN and its projections to other nuclei are somatotopically organized (Nambu, 2011). Studies in animal models of PD suggest that PD is accompanied by substantial synaptic reorganization (Fan et al., 2012; Miguelez et al., 2012; Chu et al., 2017; Pamukcu et al., 2020; Madadi Asl et al., 2022), and that the somatotopic organization of synaptic inputs is perturbed (Filion et al., 1988; Boraud et al., 2000; Cho et al., 2002). In our LIF model, stimulation-induced downregulation of synaptic connections leads to a stabilization of desynchronized neuronal activity. Evidence from other computational studies suggests that CR is not only able to downregulate synaptic connections, but also to restore more complex physiological connectivity (Hauptmann and Tass, 2010). This is important as more complex synaptic reshaping might accompany long-lasting therapeutic effects of CR stimulation in animal models and PD patients (Tass et al., 2012; Adamchic et al., 2014; Wang et al., 2016; Bore et al., 2022; Wang et al., 2022). The qualitatively different results for networks dominated by local as opposed to by non-local connections in our LIF model suggest that one might observed different stimulation-induced dynamics of connections between neurons representing similar versus different body parts. In the future, we plan on studying CR stimulation in more detailed computational models to understand how more complex synaptic connectivity affects long-lasting effects.

Our computational study reveals that knowledge of fundamental characteristics of target neuronal networks with spatially dependent synaptic connections and STDP can be harnessed to adjust CR stimulation parameters for improved synaptic downregulation. Downregulation of pathological synaptic connectivity (Tass and Majtanik, 2006) is currently assumed to be the mechanism by which CR stimulation induces long-lasting therapeutic effects in PD patients and related animal models (Tass et al., 2012; Adamchic et al., 2014; Wang et al., 2016; Pfeifer et al., 2021; Wang et al., 2022). Thereby, our results provide important hypotheses for future clinical studies and computational studies on CR stimulation of detailed computational models of target brain regions for deep brain stimulation, e.g., the STN, or vibrotactile stimulation, e.g., sensorimotor cortical areas.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary material, further inquiries can be directed to the corresponding author. The code that was used to generate the figures in this article is available on github

Author contributions

JK: Conceptualization, Data curation, Formal Analysis, Investigation, Methodology, Software, Validation, Visualization, Writing–original draft, Writing–review and editing. PT: Conceptualization, Funding acquisition, Project administration, Resources, Supervision, Validation, Visualization, Writing–review and editing, Investigation, Methodology.

Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. PT gratefully acknowledges funding support by the Vibrotactile Therapy Research Fund, by the John A. Blume Foundation, by the Alda Parkinson’s Research Fund, and by the Foundation for OCD Research (New Venture Fund 011665-2020-08-01, url: https://www.ffor.org/). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Some of the computing for this project was performed on Stanford’s Sherlock Computing cluster.

Acknowledgments

We are grateful to Stanford University and Stanford’s Sherlock Computing cluster for computational resources and support that contributed to these research results.

Conflict of interest

JK and PT filed a Stanford-owned provisional patent related to the presented results.

The author(s) declared that they were an editorial board member of Frontiers, at the time of submission. This had no impact on the peer review process and the final decision.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnetp.2024.1351815/full#supplementary-material

References

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Afsharpour, S. (1985). Light microscopic analysis of golgi-impregnated rat subthalamic neurons. J. Comp. Neurol. 236, 1–13. doi:10.1002/cne.902360102

PubMed Abstract | CrossRef Full Text | Google Scholar

Anil, S., Lu, H., Rotter, S., and Vlachos, A. (2023). Repetitive transcranial magnetic stimulation (rTMS) triggers dose-dependent homeostatic rewiring in recurrent neuronal networks. PLoS Comput. Biol. 19, e1011027. doi:10.1371/journal.pcbi.1011027

PubMed Abstract | CrossRef Full Text | Google Scholar

Barnikol, U. B., Popovych, O. V., Hauptmann, C., Sturm, V., Freund, H.-J., and Tass, P. A. (2008). Tremor entrainment by patterned low-frequency stimulation. Phil. Trans. R. Soc. A 366, 3545–3573. doi:10.1098/rsta.2008.0104

PubMed Abstract | CrossRef Full Text | Google Scholar

Bartsch, R. P., Liu, K. K. L., Bashan, A., and Ivanov, P. C. (2015). Network physiology: how organ systems dynamically interact. PLoS ONE 10, e0142143. doi:10.1371/journal.pone.0142143

PubMed Abstract | CrossRef Full Text | Google Scholar

Bashan, A., Bartsch, R. P., Kantelhardt, J. W., Havlin, S., and Ivanov, P. C. (2012). Network physiology reveals relations between network topology and physiological function. Nat. Commun. 3, 702. doi:10.1038/ncomms1705

PubMed Abstract | CrossRef Full Text | Google Scholar

Berner, R., Fialkowski, J., Kasatkin, D., Nekorkin, V., Yanchuk, S., and Schöll, E. (2019). Hierarchical frequency clusters in adaptive networks of phase oscillators. Chaos 29, 103134. doi:10.1063/1.5097835

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Berner, R., and Yanchuk, S. (2021). Synchronization in networks with heterogeneous adaptation rules and applications to distance-dependent synaptic plasticity. Front. Appl. Math. Stat. 7, 714978. doi:10.3389/fams.2021.714978

CrossRef Full Text | Google Scholar

Bevan, M. D., and Wilson, C. J. (1999). Mechanisms underlying spontaneous oscillation and rhythmic firing in rat subthalamic neurons. J. Neurosci. 19, 7617–7628. doi:10.1523/JNEUROSCI.19-17-07617.1999

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Boraud, T., Bezard, E., Bioulac, B., and Gross, C. (2000). Ratio of inhibited-to-activated pallidal neurons decreases dramatically during passive limb movement in the mptp-treated monkey. J. Neurophysiol. 83, 1760–1763. doi:10.1152/jn.2000.83.3.1760

PubMed Abstract | CrossRef Full Text | Google Scholar

Bore, J. C., Campbell, B. A., Cho, H., Pucci, F., Gopalakrishnan, R., Machado, A. G., et al. (2022). Long-lasting effects of subthalamic nucleus coordinated reset deep brain stimulation in the non-human primate model of parkinsonism: a case report. Brain Stimul. 15, P598–P600. doi:10.1016/j.brs.2022.04.005

CrossRef Full Text | Google Scholar

Cho, J., Duke, D., Manzino, L., Sonsalla, P. K., and West, M. O. (2002). Dopamine depletion causes fragmented clustering of neurons in the sensorimotor striatum: evidence of lasting reorganization of corticostriatal input. J. Comp. Neurol. 452, 24–37. doi:10.1002/cne.10349

PubMed Abstract | CrossRef Full Text | Google Scholar

Chu, H.-Y., McIver, E. L., Kovaleski, R. F., Atherton, J. F., and Bevan, M. D. (2017). Loss of hyperdirect pathway cortico-subthalamic inputs following degeneration of midbrain dopamine neurons. Neuron 95, P1306–P1318. doi:10.1016/j.neuron.2017.08.038

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Emmi, A., Antonini, A., Macchi, V., Porzionato, A., and De Caro, R. (2020). Anatomy and connectivity of the subthalamic nucleus in humans and non-human primates. Front. Neuroanat. 14, 13. doi:10.3389/fnana.2020.00013

PubMed Abstract | CrossRef Full Text | Google Scholar

Fan, K. Y., Baufreton, J., Surmeier, D. J., Chan, C. S., and Bevan, M. D. (2012). Proliferation of external globus pallidus-subthalamic nucleus synapses following degeneration of midbrain dopamine neurons. J. Neurosci. 32, 13718–13728. doi:10.1523/JNEUROSCI.5750-11.2012

PubMed Abstract | CrossRef Full Text | Google Scholar

Filion, M., Tremblay, L., and Bédard, P. J. (1988). Abnormal influences of passive limb movement on the activity of globus pallidus neurons in parkinsonian monkeys. Brain. Res. 444, 165–176. doi:10.1016/0006-8993(88)90924-9

PubMed Abstract | CrossRef Full Text | Google Scholar

FitzHugh, R. (1961). Impulses and physiological states in theoretical models of nerve membrane. Biophys. J. 1, 445–466. doi:10.1016/S0006-3495(61)86902-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Hammond, C., Bergman, H., and Brown, P. (2007). Pathological synchronization in Parkinson’s disease: networks, models and treatments. Trends Neurosci. 30, 357–364. doi:10.1016/j.tins.2007.05.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Harnack, D., Winter, C., Meissner, W., Reum, T., Kupsch, A., and Morgenstern, R. (2004). 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. doi:10.1016/j.jneumeth.2004.04.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Hauptmann, C., and Tass, P. A. (2010). Restoration of segregated, physiological neuronal connectivity by desynchronizing stimulation. J. Neural Eng. 7, 056008. doi:10.1088/1741-2560/7/5/056008

PubMed Abstract | CrossRef Full Text | Google Scholar

Hellwig, B. (2000). A quantitative analysis of the local connectivity between pyramidal neurons in layers 2/3 of the rat visual cortex. Biol. Cybern. 82, 111–121. doi:10.1007/PL00007964

PubMed Abstract | CrossRef Full Text | Google Scholar

Humphries, M. D., Wood, R., and Gurney, K. (2010). Reconstructing the three-dimensional gabaergic microcircuit of the striatum. PLoS Comput. Biol. 6, e1001011. doi:10.1371/journal.pcbi.1001011

PubMed Abstract | CrossRef Full Text | Google Scholar

Ivanov, P. C., Liu, K. K., and Bartsch, R. P. (2016). Focus on the emerging new fields of network physiology and network medicine. New J. Phys. 18, 100201. doi:10.1088/1367-2630/18/10/100201

PubMed Abstract | CrossRef Full Text | Google Scholar

Khaledi-Nasab, A., Kromer, J. A., and Tass, P. A. (2021). Long-lasting desynchronization effects of coordinated reset stimulation improved by random jitters. Front. Physiol. 12, 719680. doi:10.3389/fphys.2021.719680

PubMed Abstract | CrossRef Full Text | Google Scholar

Khaledi-Nasab, A., Kromer, J. A., and Tass, P. A. (2022). Long-lasting desynchronization of plastic neuronal networks by double-random coordinated reset stimulation. Front. Netw. Physiol. 2, 864859. doi:10.3389/fnetp.2022.864859

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Kozyrev, V., Staadt, R., Eysel, U. T., and Jancke, D. (2018). TMS-induced neuronal plasticity enables targeted remodeling of visual cortical maps. Proc. Natl. Acad. Sci. U.S.A. 115, 6476–6481. doi:10.1073/pnas.1802798115

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

Kromer, J. A., and Tass, P. A. (2022). Synaptic reshaping of plastic neuronal networks by periodic multichannel stimulation with single-pulse and burst stimuli. PLoS Comput. Biol. 18, e1010568. doi:10.1371/journal.pcbi.1010568

PubMed Abstract | CrossRef Full Text | Google Scholar

Kromer, J. A., and Tass, P. A. (2024). Sequences and their shuffling may crucially impact coordinated reset stimulation – a theoretical study. Brain Stimul. 17, P194–P196. doi:10.1016/j.brs.2024.02.004

CrossRef Full Text | Google Scholar

Kuramoto, Y. (1984). Chemical oscillations, waves, and turbulence. Berlin: Springer.

Google Scholar

Lévesque, J.-C., and Parent, A. (2005). GABAergic interneurons in human subthalamic nucleus. Mov. Disord. 20, 574–584. doi:10.1002/mds.20374

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Lysyansky, B., Popovych, O. V., and Tass, P. A. (2013). Optimal number of stimulation contacts for coordinated reset neuromodulation. Front. Neuroeng. 6, 5. doi:10.3389/fneng.2013.00005

PubMed Abstract | CrossRef Full Text | Google Scholar

Madadi Asl, M., Vahabie, A. H., Valizadeh, A., and Tass, P. A. (2022). Spike-timing-dependent plasticity mediated by dopamine and its role in Parkinson’s disease pathophysiology. Front. Netw. 2, 817524. doi:10.3389/fnetp.2022.817524

CrossRef Full Text | Google Scholar

Madadi Asl, M., Valizadeh, A., and Tass, P. A. (2016). Dendritic and axonal propagation delays determine emergent structures of neuronal networks with plastic synapses. Sci. Rep. 6, 39682. doi:10.3389/fphys.2018.01849

CrossRef Full Text | Google Scholar

Madadi Asl, M., Valizadeh, A., and Tass, P. A. (2018a). Delay-induced multistability and loop formation in neuronal networks with spike-timing-dependent plasticity. Sci. Rep. 8, 12068. doi:10.1038/s41598-018-30565-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Madadi Asl, M., Valizadeh, A., and Tass, P. A. (2018b). Propagation delays determine neuronal activity and synaptic connectivity patterns emerging in plastic neuronal networks. Chaos 28, 106308. doi:10.1063/1.5037309

PubMed Abstract | CrossRef Full Text | Google Scholar

Madadi Asl, M., Valizadeh, A., and Tass, P. A. (2023). Decoupling of interacting neuronal populations by time-shifted stimulation through spike-timing-dependent plasticity. PLoS Comput. Biol. 19, e1010853. doi:10.1371/journal.pcbi.1010853

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Mallet, N., Delgado, L., Chazalon, M., Miguelez, C., and Baufreton, J. (2019). Cellular and synaptic dysfunctions in Parkinson’s disease: stepping out of the striatum. Cells 8, 1005. doi:10.3390/cells8091005

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

McIntyre, C. C., Mori, S., Sherman, D. L., Thakor, N. V., and Vitek, J. L. (2004). Electric field and stimulating influence generated by deep brain stimulation of the subthalamic nucleus. Clin. Neurophysiol. 115, 589–595. doi:10.1016/j.clinph.2003.10.033

PubMed Abstract | CrossRef Full Text | Google Scholar

Miguelez, C., Morin, S., Martinez, A., Goillandeau, M., Bezard, E., Bioulac, B., et al. (2012). Altered pallido-pallidal synaptic transmission leads to aberrant firing of globus pallidus neurons in a rat model of Parkinson’s disease. J. Physiol. 590, 5861–5875. doi:10.1113/jphysiol.2012.241331

PubMed Abstract | CrossRef Full Text | Google Scholar

Morrison, A., Diesmann, M., and Gerstner, W. (2008). Phenomenological models of synaptic plasticity based on spike timing. Biol. Cybern. 98, 459–478. doi:10.1007/s00422-008-0233-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Nagumo, J., Arimoto, S., and Yoshizawa, S. (1962). An active pulse transmission line simulating nerve axon. Proc. IRE 50, 2061–2070. doi:10.1109/JRPROC.1962.288235

CrossRef Full Text | Google Scholar

Nambu, A. (2011). Somatotopic organization of the primate basal ganglia. Front. Neuroanat. 5, 26. doi:10.3389/fnana.2011.00026

PubMed Abstract | CrossRef Full Text | Google Scholar

Neumann, W.-J., Degen, K., Schneider, G.-H., Brücke, C., Huebl, J., Brown, P., et al. (2016). Subthalamic synchronized oscillatory activity correlates with motor impairment in patients with Parkinson’s disease. Mov. Disord. 31, 1748–1751. doi:10.1002/mds.26759

PubMed Abstract | CrossRef Full Text | Google Scholar

Pamukcu, A., Cui, Q., Xenias, H. S., Berceau, B. L., Augustine, E. C., Fan, I., et al. (2020). Parvalbumin+ and npas1+ pallidal neurons have distinct circuit topology and function. J. Neurosci. 40, 7855–7876. doi:10.1523/JNEUROSCI.0361-20.2020

PubMed Abstract | CrossRef Full Text | Google Scholar

Pariz, A., Trotter, D., Hutt, A., and Lefebvre, J. (2023). Selective control of synaptic plasticity in heterogeneous networks through transcranial alternating current stimulation (tACS). PLoS Comput. Biol. 19, e1010736. doi:10.1371/journal.pcbi.1010736

PubMed Abstract | CrossRef Full Text | Google Scholar

Pfeifer, K. J., Kromer, J. A., Cook, A. J., Hornbeck, T., Lim, E. A., Mortimer, B. J. P., et al. (2021). Coordinated reset vibrotactile stimulation induces sustained cumulative benefits in Parkinson’s disease. Front. Physiol. 12, 624317. doi:10.3389/fphys.2021.624317

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Richardson, K. A., Gluckman, B. J., Weinstein, S. L., Glosch, C. E., Moon, J. B., Gwinn, R. P., et al. (2003). In vivo modulation of hippocampal epileptiform activity with radial electric fields. Epilepsia 44, 768–777. doi:10.1046/j.1528-1157.2003.35402.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Richter, E. O., Hoque, T., Halliday, W., Lozano, A. M., and Saint-Cyr, J. A. (2004). Determining the position and size of the subthalamic nucleus based on magnetic resonance imaging results in patients with advanced Parkinson disease. J. Neurosurg. 100, 541–546. doi:10.3171/jns.2004.100.3.0541

PubMed Abstract | CrossRef Full Text | Google Scholar

Rosenblum, M., Pikovsky, A., Kurths, J., Schäfer, C., and Tass, P. A. (2001). “Chapter 9 Phase synchronization: from theory to data analysis,” in Neuro-informatics and neural modelling, handbook of biological physics. Editors S. Gielen, and F. Moss (Amsterdam: Elsevier), 4, 279–321. doi:10.1016/s1383-8121(01)80012-9

CrossRef Full Text | Google Scholar

Rubin, J., Lee, D. D., and Sompolinsky, H. (2001). Equilibrium properties of temporally asymmetric hebbian plasticity. Phys. Rev. Lett. 86, 364–367. doi:10.1103/PhysRevLett.86.364

PubMed Abstract | CrossRef Full Text | Google Scholar

Schmalz, J., and Kumar, G. (2019). Controlling synchronization of spiking neuronal networks by harnessing synaptic plasticity. Front. Comput. Neurosci. 13, 61. doi:10.3389/fncom.2019.00061

PubMed Abstract | CrossRef Full Text | Google Scholar

Schwab, B. C., Misselhorn, J., and Engel, A. K. (2019). Modulation of large-scale cortical coupling by transcranial alternating current stimulation. Brain Stimul. 12, 1187–1196. doi:10.1016/j.brs.2019.04.013

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Shen, K.-Z., Zhu, Z.-T., Munhall, A., and Johnson, S. W. (2003). Synaptic plasticity in rat subthalamic nucleus induced by high-frequency stimulation. Synapse 50, 314–319. doi:10.1002/syn.10274

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Spix, T. A., Nanivadekar, S., Toong, N., Kaplow, I. M., Isett, B. R., Goksen, Y., et al. (2021). Population-specific neuromodulation prolongs therapeutic benefits of deep brain stimulation. Science 374, 201–206. doi:10.1126/science.abi7852

PubMed Abstract | CrossRef Full Text | Google Scholar

Tass, P. A. (2003a). Desynchronization by means of a coordinated reset of neural sub-populations: a novel technique for demand-controlled deep brain stimulation. Prog. Theor. Phys. Supp. 150, 281–296. doi:10.1143/PTPS.150.281

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Tass, P. A., and Hauptmann, C. (2006). Long-term anti-kindling effects induced by short-term, weak desynchronizing stimulation. Nonlinear Phenom. Complex Syst. 9, 298–312.

Google Scholar

Tass, P. A., and Hauptmann, C. (2007). Therapeutic modulation of synaptic connectivity with desynchronizing brain stimulation. Int. J. Psychophysiol. 64, 53–61. doi:10.1016/j.ijpsycho.2006.07.013

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Vinding, M. C., Waldthaler, J., Eriksson, A., Manting, C. L., Ferreira, D., Ingvar, M., et al. (2024). Oscillatory and non-oscillatory features of the magnetoencephalic sensorimotor rhythm in Parkinson’s disease. NPJ Parkinson’s Dis. 10, 51. doi:10.1038/s41531-024-00669-3

CrossRef Full Text | Google Scholar

Wang, J., Fergus, S. P., Johnson, L. A., Nebeck, S. D., Zhang, J., Kulkarni, S., et al. (2022). Shuffling improves the acute and carryover effect of subthalamic coordinated reset deep brain stimulation. Front. Neurol. 13, 716046. doi:10.3389/fneur.2022.716046

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar