Long-Term Desynchronization by Coordinated Reset Stimulation in a Neural Network Model With Synaptic and Structural Plasticity

Several brain disorders are characterized by abnormal neuronal synchronization. To specifically counteract abnormal neuronal synchrony and, hence, related symptoms, coordinated reset (CR) stimulation was computationally developed. In principle, successive epochs of synchronizing and desynchronizing stimulation may reversibly move neural networks with plastic synapses back and forth between stable regimes with synchronized and desynchronized firing. Computationally derived predictions have been verified in pre-clinical and clinical studies, paving the way for novel therapies. However, as yet, computational models were not able to reproduce the clinically observed increase of desynchronizing effects of regularly administered CR stimulation intermingled by long stimulation-free epochs. We show that this clinically important phenomenon can be computationally reproduced by taking into account structural plasticity (SP), a mechanism that deletes or generates synapses in order to homeostatically adapt the firing rates of neurons to a set point-like target firing rate in the course of days to months. If we assume that CR stimulation favorably reduces the target firing rate of SP, the desynchronizing effects of CR stimulation increase after long stimulation-free epochs, in accordance with clinically observed phenomena. Our study highlights the pivotal role of stimulation- and dosing-induced modulation of homeostatic set points in therapeutic processes.

Introduction

High-frequency deep brain stimulation (HF DBS) is the standard treatment of medically refractory movement disorders such as Parkinson’s disease (PD) (Benabid et al., 1991; Krack et al., 2003; Deuschl et al., 2006) and is also being tested in other psychiatric diseases (Cleary et al., 2015). HF DBS aims to permanently deliver electrical charge-balanced pulses at high frequencies (>100 Hz) to target areas such as the thalamic ventralis intermedius (VIM) nucleus or the subthalamic nucleus (STN) through chronically implanted depth electrodes (Benabid et al., 1991, 2009). However, the precise mechanism of action of standard HF DBS is not yet sufficiently understood (Johnson et al., 2008; Gradinaru et al., 2009; Deniau et al., 2010). Standard HF DBS has acute clinical (Temperli et al., 2003) and electrophysiological (Kühn et al., 2008) effects, which are present only during stimulation. A number of reversible as well as non-reversible adverse events (AE) can arise from DBS treatment, which is itself an invasive treatment; these can range from neurological AEs such as gait disturbances and speech problems, to psychiatric AEs such as depression (Baizabal-Carvallo and Jankovic, 2016; Buhmann et al., 2017). Therefore, any DBS treatment regime should be minimized in duration as much as possible, in order to also minimize the AEs that can result (Buhmann et al., 2017).

Abnormal neuronal synchronization is a hallmark of PD and may occur in different frequency ranges (Lenz et al., 1994; Nini et al., 1995; Llinás et al., 1999; Brown, 2002). Different stimulation protocols showed that trains of HF DBS may in fact modulate synaptic plasticity (Shen et al., 2003; Prescott et al., 2009; Milosevic et al., 2018). There are several different mechanisms of synaptic plasticity, of which spike-timing-dependent plasticity (STDP) is the most well-known and widely accepted (see, e.g., Lisman, 2017). STDP is a fundamental synaptic plasticity mechanism which modifies synaptic strengths based on the relative timing of pre- and postsynaptic spike pairs (Gerstner et al., 1996; Markram et al., 1997; Bi and Poo, 1998, 2001; Song et al., 2000). According to the STDP mechanism, the synaptic strength is potentiated when the postsynaptic spike follows the presynaptic spike; conversely, the synaptic strength is depressed when the postsynaptic spike advances the presynaptic spike (Markram et al., 1997).

Coordinated reset (CR) stimulation is a computationally developed patterned multichannel stimulation (Tass, 2003a,b) that employs STDP to induce long-lasting therapeutic effects (Tass and Majtanik, 2006). CR stimulation specifically counteracts abnormal neuronal synchrony by desynchronization (Tass, 2003a,b), which causes a reduction of the rate of coincidences and, mediated by STDP, a reduction of the synaptic weights (Tass and Majtanik, 2006). Consequently, as shown computationally, neuronal networks are shifted from unfavorable attractors (with strong neuronal synchrony and synaptic connectivity) to more favorable attractors (with weak synchrony and reduced synaptic weights) (Tass and Majtanik, 2006). Accordingly, in neuronal network models a long-lasting desynchronization can be achieved. The term long-lasting desynchronization refers to desynchronization which exceeds cessation of stimulation (see Tass and Majtanik, 2006; Manos et al., 2018a; Kromer and Tass, 2020). Such effects have been demonstrated pre-clinically and clinically in a variety of studies (e.g., Tass et al., 2009, 2012a,b; Adamchic et al., 2012, 2013; Pfeifer et al., 2021). For CR-DBS delivered to the STN, computationally predicted desynchronizing effects (Tass, 2003a,b), cumulative effects (Hauptmann and Tass, 2009) and long-lasting effects (Tass and Majtanik, 2006) were verified in pre-clinical studies in a 1-methyl-4-phenyl-1,2,3,6-tetrahydropyridine (MPTP) monkey PD model (Tass et al., 2012b; Wang et al., 2016) as well as in a proof-of-concept study in PD patients (Adamchic et al., 2014a).

In accordance with computational predictions (Popovych and Tass, 2012; Tass and Popovych, 2012), CR stimulation can also be realized by means of sensory stimulation modalities and induce long-lasting desynchronization. In a proof-of-concept study in patients with chronic subjective tinnitus, acoustic CR stimulation caused a significant and sustained reduction of tinnitus symptoms (Tass et al., 2012a) combined with a significant reduction of abnormal neuronal synchrony (Adamchic et al., 2012, 2013; Tass et al., 2012a), abnormal cross-frequency coupling (Adamchic et al., 2014b) and abnormal effective connectivity (Silchenko et al., 2013). In a study of PD patients, vibrotactile CR (vCR) stimulation (Tass, 2017) delivered to the fingertips during three consecutive days for 4 h per day caused improvements in gait and bradykinesia for up to one month post-stimulation (Syrkin-Nikolau et al., 2018). A case series in three PD patients treated with daily regular or noisy vCR fingertip stimulation during 6+ months revealed continuous improvement of PD motor symptoms as assessed after medication withdrawal (i.e., off medication) (Pfeifer et al., 2021). Regular vCR uses periodic stimulus delivery, whereas for noisy vCR stimulus onsets had a moderate jitter in time (Pfeifer et al., 2021). In addition, in a pilot study six PD patients received daily noisy vCR stimulation during 3 months (Pfeifer et al., 2021). Motor evaluations and at-rest electroencephalographic (EEG) recordings, performed off medication before and after the 3-month vCR therapy, revealed a statistically and clinically significant motor improvement along with a significant decrease of cortical sensorimotor high beta power (21–30 Hz). Both studies with sensory CR performed so far (Tass et al., 2012a; Pfeifer et al., 2021) provide preliminary evidence indicating that the effects of repetitively administered acoustic as well as vibrotactile CR stimulation may increase in the course of the treatment, in this way enabling longer treatment pauses and reduced dosage regimens. However, in the computational studies performed as yet, memory-type effects, mediating an increase of CR efficacy with repetitive stimulation delivery, have not been described. We hypothesize that structural plasticity, acting on a slower time scale than synaptic plasticity (Fauth and Tetzlaff, 2016), might play a role in modifying CR effects in the course of repeated delivery of CR stimuli.

Apart from synaptic plasticity, structural plasticity also appears to play an important role in PD. Structural plasticity allows neurons to establish new or delete pre-existing synaptic connections; this occurs via the extension or retraction of axons and dendrites, or by modifying the number of axonal boutons or dendritic spines (Butz et al., 2009). Changes in structural connectivity and rewiring of connections is crucial for memory formation (Xu et al., 2009) as well as for adaptive and maladaptive processes in response to central and peripheral lesions of the nervous system (Keck et al., 2008; Yamahachi et al., 2009; Brown et al., 2010) and neurodegeneration (Deller et al., 2006). A study of the gray matter value in the basal ganglia (BG) of patients with symptomatic PD revealed that BG undergo a form of progressive atrophy that increases with disease duration and severity (Reetz et al., 2009). Furthermore, it was shown that balance training induces morphometric changes in PD patients (Sehm et al., 2014): a voxel-based morphometry revealed performance improvement-correlated gray matter changes in the right anterior precuneus, left inferior parietal cortex, left ventral premotor cortex, bilateral anterior cingulate cortex, and left middle temporal gyrus as well as time-dependent gray matter changes in the right cerebellum. In contrast, only learning-dependent gray matter changes in the left hippocampus were observed in healthy controls (Sehm et al., 2014). In a single PD patient case study, it was shown that long-term HF DBS gave rise to significant localized structural changes in sensory-motor, prefrontal, limbic, and olfactory brain regions together with a restoration of functional connectivity (van Hartevelt et al., 2014). In addition, a computational analysis revealed that the observed changes of structural weights are compliant with a hypothesized Hebbian-like mechanism (van Hartevelt et al., 2015). For a more extensive discussion on the role of structural plasticity in the brain please see the corresponding section of the Supplementary Material and van Ooyen and Butz-Ostendorf (2017).

Accordingly, we set out to computationally integrate a structural plasticity mechanism into our modeling approach; thus, in addition to STDP we also used the built-in structural plasticity module (Diaz-Pier et al., 2016) of the NEST platform (Bos et al., 2015). NEST is a simulator for neuronal modeling, and one which can incorporate synaptic plasticity modules. Structural plasticity has two main areas of application in NEST. The primary purpose is to model the neurobiological phenomenon of morphological transformations that a neuron undergoes over time, leading to the creation or deletion of synapses. For a more complete description of the model see Butz and van Ooyen (2013) and for the details on the exact implementation see Diaz-Pier et al. (2016). Specifically, the algorithm serves a homeostatic purpose by self-generating or self-eliminating connections in order to reach and maintain the target firing rate of the network. New synapses are plastic, and their weight-values are drawn from the same probability distribution defined in the initialization of the simulation. In order to achieve this, the algorithm follows a set of homeostatic rules which dictate how the connectivity should be modified in order to achieve the desired levels of electrical activity. The algorithm works at the single neuron level, so that each neuron follows this rule independently. Contrary to STDP, the structural plasticity mechanism allows connections to be deleted and created between neurons, and not only alter their weights. By combining both plasticity mechanisms, we increased the sampling capacities and the plasticity potential of the system to explore different connectivity configurations, allowing long-term structural changes and short-term learning in the same simulation.

We hypothesize that structural plasticity may modify long-term effects of CR stimulation. Motivated by experimental findings indicating an increased efficacy of sensory CR stimulation in the course of the therapy (Tass et al., 2012a; Pfeifer et al., 2021), we incorporate structural plasticity in our computational model and consider different treatment scenarios, depending on whether CR stimulation may induce changes of the set point of structural plasticity. The efficacy of CR stimulation after long pauses without stimulation may strongly increase provided CR stimulation causes a reduction of the set point. The computational approach presented here may provide predictions that contribute to optimized and experimentally testable dosage regimens, which specifically aim at modulating homeostatic set points of structural plasticity mechanisms.

Materials and Methods

The Terman-Rubin Neuron Model

The Terman and Rubin (Terman et al., 2002; Rubin and Terman, 2004) single-compartment conductance-based model is used for the description of the neuronal activity of the individual STN and globus pallidus external (GPe) neurons. The membrane potential is given by the following equation:

$$c_{\text{m}}\frac{\text{d}v}{\text{d}t}=-I_{\text{L}}-I_{\text{K}}-I_{\text{Na}}-I_{\text{T}}-I_{\text{Ca}}-I_{\text{ahp}}-I_{\text{syn}}+I_{\text{stim}}+I_{\text{noise}}(1)$$

Spiking activity is caused by the sodium (Na⁺) and potassium (K⁺) ionic currents I_Na, I_K, I_T, and I_Cadescribe the low-threshold T-type and high-threshold Ca²⁺ current, respectively. I_ahp represents a Ca²⁺-activated, voltage-independent after-hyperpolarization K⁺ current and I_L the leak current. In addition, the STN and GPe neurons are influenced by synaptic inputs I_syn. Surrounding brain areas contribute with I_noise while there is an external stimulation current I_stim (only in the STN population) which models the deep brain stimulation (DBS). Other ionic currents (in pA/μm²) are described by Terman et al. (2002); Rubin and Terman (2004), and Ebert et al. (2014).

$$I_{\text{L}}=g_{\text{L}}\left[v-v_{\text{L}}\right](2)$$

$$I_{\text{K}}=g_{\text{K}}{n}^4\left[v-v_{\text{K}}\right](3)$$

$$I_T=g_{\text{T}}{a}_{\mathrm{\infty }}^3\left(v\right)b_{\mathrm{\infty }}^2\left(r\right)\left[v-v_{\text{Ca}}\right](4)$$

$$I_{\text{Ca}}=g_{\text{Ca}}{s}_{\mathrm{\infty }}^2\left(v\right)h\left[v-v_{\text{Ca}}\right](5)$$

$$I_{\text{ahp}}=g_{\text{ahp}}\left[v-v_{\text{K}}\right]\frac{\left[\text{Ca}\right]}{\left[\text{Ca}\right]+k_1}(6)$$

$$\frac{\text{d}[\text{Ca}]}{\text{d}t}=\text{Γ}\left(-I_{\text{Ca}}-I_{\text{T}}-k_{\text{Ca}}[\text{Ca}]\right)(7)$$

Subthalamic nucleus and GPe neurons are described by similar equations and they differ only in a few parameter values (see Supplementary Tables 1, 2) as well as in the form of the low threshold Ca²⁺ current:

$$I_{\text{T}}=g_{\text{T}}{a}_{\mathrm{\infty }}^3\left(v\right)r\left[v-v_{\text{Ca}}\right](8)$$

for the GPe neurons (i.e., the $b_{\mathrm{\infty }}^2$ term is omitted). Moreover, neurons are not identical, i.e., their reverse potential parameters and ion channel maximum conductances are drawn from Gaussian probability distributions with 10% standard deviation around the mean value (except the parameter v_L where we draw values with 0.15% standard deviation). All the mean parameter values and their units are summarized in Supplementary Tables 1, 2. The slowly operating gating variables n, h, and r depend on time and voltage and their first-order kinetics are governed by differential equations of the form:

$$\frac{\text{d}X}{\text{d}t}={\mathrm{\phi }}_X\left[\frac{X_{\mathrm{\infty }}\left(v\right)-X}{{\mathrm{\tau }}_X\left(v\right)}\right],(9)$$

Where X:n,h,r with:

$${\mathrm{\tau }}_X\left(v\right)={\mathrm{\tau }}_X^0+\frac{{\mathrm{\tau }}_X^1}{1+\text{exp}\left[-\frac{v-{\mathrm{\theta }}_X^{\mathrm{\tau }}}{o_X^{\mathrm{\tau }}}\right]},(10)$$

where θ^τ is the voltage at which the time constant is midway between its maximum and minimum values, and o^τ is the slope factor for the voltage dependence of the time constant. Activation gating for the rapidly activating channels (m, a, and s) was treated as instantaneous. For all gating variables X = n,m,h,a,r, or s, the steady-state voltage dependence was determined using:

$$X_{\mathrm{\infty }}\left(v\right)=\frac{1}{1+\exp \left[-\frac{v-{\mathrm{\theta }}_X}{o_X}\right]},(11)$$

where θ_X is the half activation (or inactivation) voltage for gating variable X, and o_X is the slope factor for that variable. For the T current inactivation variable b, we used:

$$b_{\mathrm{\infty }}\left(r\right)=\frac{1}{1+\text{exp}\left[\frac{r-{\mathrm{\theta }}_b}{o_b}\right]}-\frac{1}{1+\text{exp}\left[-\frac{{\mathrm{\theta }}_b}{o_b}\right]},(12)$$

See Terman et al. (2002) and Rubin and Terman (2004) for more information on the model and the variables’ units. One main difference of our study compared to Ebert et al. (2014) is the initialization of the STN synaptic weights. In Ebert et al. (2014), the initial mean synaptic weight was set at ${\overline{w}}_{\text{ss}}^{\text{STN}}$ = 0.018, in order to tune the system at a strongly synchronized state. Here, we chose a different approach. We set ${\overline{w}}_{\text{ss}}^{\text{STN}}$ = 0.0025, at a much lower value, and drove the system to a similar strongly synchronized state by applying periodic stimulation (PS) to the STN neurons.

In our NEST implementation, which is based on NEST 2.10, we have incorporated an additional feature which allows us to enable and disable synaptic and structural plasticity. We have also modified the topology module in order to allow plastic synapses to be created. The Terman-Rubin neuron model (STN and GPe populations) was implemented in the NESTML platform (Plotnikov et al., 2016; Blundell et al., 2018; Perun et al., 2018).

The Network Model

Our network model consists of two interacting nuclei, the STN and the GPe. The STN neuron population is linked with the cortex in an excitatory manner while the GPe population is linked with the striatum in an inhibitory manner, see Figure 1A. The connectivity matrix of the STN-GPe network can be expressed as the combination of several sub-networks:

FIGURE 1

Figure 1. Network, synaptic, and structural plasticity setup. (A) Structure of the model network. Excitatory synaptic connections are shown in red (with arrows) and inhibitory connections in blue (with squares). The connections whose synaptic weights do not vary in time are in solid lines. The red dashed line represents time dependent connections between the subthalamic nucleus (STN) neurons, modified by spike time-dependent plasticity (STDP) and structural plasticity. The external electrical stimulation delivered to the STN is indicated by the gray solid line. The Poissonian spike train input from the cortex to the STN is constantly active. (B) The three-dimensional STN-GPe model (left brain hemisphere). Segmented MRI slices taken from a Parkinson’s disease (PD) patient prior to deep brain stimulation (DBS) surgery. Each ellipsoid volume contains 1,000 randomly distributed point-like neurons. Note, neurons are not overlapping. (C) STDP rule for STN neurons. Synaptic connection potentiation is highlighted in red and depression in blue. (D) Growth rate curves determining the rate of creation and deletion of synaptic elements in the structural plasticity model for STN-STN connections (red) and STN-GPe connections (blue). The parameters which define the shape of the curve are the growth rate ν_SP, which defines the peak of the curve, and the target firing rate ε. The light red band indicates the region where the model operates during the simulations reported in this study. These Gaussian curves correspond to Eq. (19) which represents the change in synaptic elements dz/dt as a function of the firing rate FR(t). In this case, the red line has a growth rate ν_SP of 0.00008 synaptic elements/ms and the blue line of 0.00002 synaptic elements/ms. The target firing rate ε in both curves is 3.0 Hz. Finally, the minimal initial firing rate η is 0.0 for both curves.

$$W=\left(\begin{array}{cc}\hfill w_{\text{ss}}\hfill & \hfill w_{gs}\hfill \\ \hfill w_{\text{sg}}\hfill & \hfill w_{\text{gg}}\hfill \end{array}\right)(13)$$

Denoting the connectivity matrices for the STN-STN (ss), STN-GPe (sg), GPe-STN (gs), and GPe-GPe (gg) sub-networks, respectively. Each W_ij value corresponds to the connection strength between neuron i and j while self-connections are excluded, i.e., W_ii = 0,∀i. Our neuronal network is a weighted and directed graph, and hence, this matrix is not necessarily symmetric.

The STN population projects long range excitatory connections on the GPe population, while the GPe population projects long range inhibitory connections to the STN neurons (Shink et al., 1996). The connections between STN and GPe were established following Baufreton et al. (2009) and using connectivity values of 2% (i.e., 200 synapses per neuron, connecting to 20 randomly chosen targets in the distant population, respectively). Following Fujimoto and Kita (1993) and Ebert et al. (2014), the transmission delay for connections between the STN and the GPe was set to δ_sg = δ_gs = 4.0ms. Recent findings on the synaptic organization of the STN population identifying its neurons as parallel processing units (Steiner et al., 2019) might help to further refine computational models of the BG.

Network and Neuron Coupling Description

The network setup and 3D anatomical arrangement of the neurons was introduced in Ebert et al. (2014) where a sample of 10,000 STN and 10,000 GPe neuron coordinates was drawn. The coordinates used for the 3D implementation originated from magnetic resonance imaging (MRI) data taken from a PD patient before DBS surgery (both STN and GPe coordinates are measured in the left-brain hemisphere while there is no actual overlap of the single neurons, Figure 1B). The rationale for choosing equal size (# of neurons) for both populations (even though the GPe area is larger than STN) is that the GPe has a lower density of neurons than the STN (see Levesque and Parent, 2005). The two populations form two ellipsoids of different size, inside of which coordinates for the neuron positions are randomly drawn. For more details on the 3D implementation with two set of coordinates and ellipsoid dimensions and its motivation, we refer to Ebert et al. (2014).

The general connectivity properties are largely based on the values used and presented in Ebert et al. (2014), where 10,000 STN and 10,000 GPe neurons were considered [more details regarding the parameter values can also be found in Gillies and Willshaw (1998) and Holgado et al. (2010)]. The main difference in this study is that we have considered a smaller in size neural network (1,000 STN and 1,000 GPe neurons instead of 2×10,000) in order to trade-off between a large-scale model and the corresponding computation time. Hence, the each STN neuron extends connections to p_ss = 70% of the entire STN population (instead of p_ss = 7% that was used for the large-scale system). The connection probability between STN neurons at a certain distance is given by the exponential function p(x) = e^−x/c_d with c_d = 0.5 and x the Euclidean distance between neurons. This probability is a non-negative function with values in the interval [0, 1]. Note, it is not a probability density which would have an integral equal to 1. The corresponding parameter values for the inhibitory GPe neurons are p_gg = 10% (instead of p_gg = 1% that was used for the large-scale system),δ_gg =  4.0 ms and c_d = 0.63. The initial synaptic weights for both STN and GPe connections are drawn from a Gaussian probability distribution around a mean value ${\overline{w}}_{\text{ss}}^{\text{STN}}=\mathrm{\hspace{0.17em}0.0025}$ with standard deviation ${\mathrm{\sigma }}_{\text{ss}}^{\text{STN}}=\mathrm{\hspace{0.17em}0.000125}$ (no distance dependence). The two structures are connected to each other as the STN affects GPe via excitatory input while the GPe exerts an inhibitory impact on the STN. There is no distance dependence for their connectivity. We adopted a probability connectivity (see Baufreton et al., 2009) of 20% (instead of 2% that was used for the large-scale system), i.e., 200 fixed randomly picked synapses with transmission delay of δ_sg = δ_gs = 6.0ms for connections between the STN and the GPe, as published in Fujimoto and Kita (1993) and Kita et al. (2005). Their corresponding initial synaptic weights are drawn again from a Gaussian probability distribution around a mean value ${\overline{w}}_{\text{sg}}^{\text{STN}}=\mathrm{\hspace{0.17em}0.006}$ and ${\overline{w}}_{\text{gs}}^{\text{STN}}=\mathrm{\hspace{0.17em}0.003}$ with standard deviation ${\mathrm{\sigma }}_{\text{sg}}^{\text{STN}}=\mathrm{\hspace{0.17em}0.0003}$ and ${\mathrm{\sigma }}_{\text{sg}}^{\text{STN}}=\mathrm{\hspace{0.17em}0.00015}$, respectively.

Using the aforementioned described parameter values, we achieved a good agreement between the two networks of different size concerning their dynamics, as assessed by macroscopic quantities, degree of synchrony (as measured by the order parameter throughout this study) and mean firing rate, we adjusted all relative connections in a similar manner. The comparison between the large network (2 × 10,000 neurons) with the reduced size one (2 × 1,000) can be found in Supplementary Figure 1 showing similar overall system evolution without any significant qualitative differences.

Synaptic Currents

The postsynaptic currents are described with α -functions (Dayan and Abbott, 2000; Gerstner and Kistler, 2002):

$$\mathrm{\alpha }(t)=\frac{t-t_k}{{\mathrm{\tau }}_{\text{syn}}^2}{\text{e}}^{-\frac{t-t_k}{{\mathrm{\tau }}_{\text{syn}}}},t_k\le t<t_{k+1}(14)$$

Where t_k denotes the spike time. The total synaptic input current to a postsynaptic neuron i received from presynaptic neurons j is then given by:

$$I_{\text{syn},i}(t)=\sum _j{W}_{ij}\left(v_i\left(t\right)-v_{\text{syn}}\right)\mathrm{\alpha }\left(t\right),(15)$$

Where W_ij is the synaptic weight (coupling strength) between presynaptic neurons j and the postsynaptic neuron i, and v_syn is the reversal potential for excitatory or inhibitory ss, sg, gs, and gg connections. Therefore, v_syn depends on the types of connected neurons and on whether the connection is excitatory or inhibitory. In our model, there exist four types of synapses for ss, sg, gs, and gg connections. All values are given in Supplementary Tables 1, 2.

Noise Inputs

Each GPe neuron receives external inhibitory input from the striatum which is described by an additional constant negative input current I_app = −7.0pAin Eq. (1) (note that for the excitatory STN neurons I_app =  0pA) and external noise I_noise described by Poissonian spike trains (form other surrounding brain structures) with a frequency of ${\text{f}}_{\text{P}}^{\text{GPe}}=\mathrm{\hspace{0.33em}40}\text{Hz}$:

$$I_{\text{noise},i}=\sum _{j=0}^{N_t}{w}_{\text{noise}}\left(v_{j,i}\left(t\right)-v_{\text{noise}}\right)\mathrm{\alpha }\left(t\right).(16)$$

N_t is an integer random variable following the Poisson distribution $P_{\mathrm{\lambda }}\left(k\right)=\frac{{\mathrm{\lambda }}^k}{k!}{\text{e}}^{-\mathrm{\lambda }}$, with parameter given by the simulated time interval λ=f_Pδ, where δ is equal to the time resolution of the simulation (0.1 ms in our simulations). In a similar manner, we model the excitatory input from the cortex to all STN neurons with a frequency of ${\text{f}}_{\text{P}}^{\text{GPe}}=\mathrm{\hspace{0.33em}20}\text{Hz}$. Each neuron receives a different random Poissonian spike train. These noise frequency values are chosen with respect to the respective STN and GPe relative firing rates, i.e., the mean firing rate of GPe is about twice the one of the STN [see Ebert et al. (2014) and references therein for motivation and justification]. The noise input synaptic weight is set at w_noise =  0.2mS/cm² while the time constants for the α -function are set at τ_noise = 1.0 ms and the reversal potential is v_noise = 0mV for both STN and GPe neurons.

Spike-Timing-Dependent Plasticity

The synaptic coupling weights (w_ij) of the STN neurons evolve dynamically and depend on the time difference ($\mathrm{\Delta }{t}_{ij}=t_j^f-t_i^f$) between the firing (onset) of the spikes of the post-synaptic neuron and pre-synaptic neuron (denoted by $t_i^f$ and $t_j^f$, respectively). When Δt_ij > 0 the rule implements synaptic potentiation due to causal relationships and when Δt_ij < 0 it implements synaptic depression. We use the following STDP rule as implemented in NEST:

$$\mathrm{\Delta }{w}_{ij}\left(\mathrm{\Delta }{t}_{ij}\right)=\{\begin{array}{c}\mathrm{\lambda }{e}^{-\frac{|\mathrm{\Delta }{t}_{ij}|}{{\mathrm{\tau }}_{+}}},\mathrm{\Delta }{t}_{ij}>0\hfill \\ -\mathrm{\lambda }\mathrm{\gamma }{e}^{-\frac{|\mathrm{\Delta }{t}_{ij}|}{{\mathrm{\tau }}_{-}}},\mathrm{\Delta }{t}_{ij}\le 0\hfill \end{array},(17)$$

where τ₋ =  27.5ms,τ₊ =  12.0ms (time constants for the synaptic weight change of depression and potentiation, respectively),λ =  0.002(learning rate of the synaptic connection) andγ = 1.4 is the ratio between depression and potentiation in the synaptic learning rule (Figure 1C). The (de)synchronized dynamics are stable with the above rule and parameter values resulting in multistability. We restrict the synaptic weights $w_{ij}^{\text{ss}}$ (within the STN population) on the interval[0,0.02], in this way avoiding a non-physiological unbounded increase or decrease. More detailed analysis and parameter exploration of this particular STPD rule and its motivation can also be found in Ebert et al. (2014).

Structural Plasticity

Using the structural plasticity framework in NEST, a network will self-generate synapses in order to stably achieve a desired profile of electrical activity, a measure that is experimentally more easily accessible than detailed connectivity data. By progressively and slowly changing the connections between neurons in the network and the weight of these connections for all populations simultaneously, the structural plasticity algorithm is able to find a stable configuration with the desired firing rate profile. The structural plasticity implementation in NEST is based on the model proposed by Butz and van Ooyen (2013) and described in technical detail by Diaz-Pier et al. (2016). In this plasticity framework, neurons have contact points called synaptic elements which increase or decrease in number according to simple homeostatic rules. This rule is originally based on the intracellular calcium concentration (Butz and van Ooyen, 2013; Diaz-Pier et al., 2016) which is modulated by the network’s firing rate. Synaptic elements typically represent axonal boutons and dendritic spines. When new synaptic elements become available, they can be used to create new synapses with other compatible elements. If the contact points are eliminated, the synapses formed earlier are destroyed. Connectivity in the network is updated on a much slower timescale than the electrical activity of neurons. In order to have a computationally efficient simulation, connectivity in the network is updated between 100 and 1,000 times more slowly than the electrical activity in the neurons. The homeostatic rules applied to the synaptic elements are intended to keep the mean firing rate stable. The firing rate FR(t) is calculated in NEST by each neuron following this equation:

$$\frac{\text{dFR}}{\text{d}t}=\{\begin{array}{c}-\frac{\text{FR}\left(t\right)}{{\mathrm{\tau }}_{\text{SP}}}+\mathrm{\beta },\mathrm{if}\mathrm{the}\text{neuron}\mathrm{fires}\hfill \\ -\frac{\text{FR}\left(t\right)}{{\mathrm{\tau }}_{\text{SP}}},\mathrm{otherwise}\hfill \end{array},(18)$$

FR is the firing rate of the neuron which is directly proportional to the internal calcium concentration of the neuron at any point in time (Diaz-Pier et al., 2016), β is the calcium intake constant which specifies the increment of the internal calcium concentration every time the neuron spikes and τ_SP is the calcium concentration decay time constant. This means that the firing rate is internally calculated by NEST for each neuron model using a low pass filtering technique on the spiking activity corresponding to the instantaneous mean value of the firing rate. For more details on how this calculation is performed as well as the impact of the update interval and other structural plasticity parameters in the evolution of the simulations see Nowke et al. (2018). The firing rate is always calculated by NEST, even when structural plasticity is not enabled. In NEST it is possible to use different types of growth curves for the synaptic elements like Linear, Gaussian, and Sigmoidal. Additionally, users can implement their own growth curves to match their specific questions. In the current study, a Gaussian curve as defined in Figure 1D describes the growth rate of connection points for neurons, i.e.:

$$\frac{\text{d}z}{\text{d}t}={\mathrm{\nu }}_{\text{SP}}\left[2{\text{e}}^{-{\left(\frac{\text{FR}\left(t\right)-\mathit{\text{ξ}}}{\mathrm{\zeta }}\right)}^2}-1\right],(19)$$

where z is the number of synaptic elements, $\mathit{\text{ξ}}=\frac{\mathrm{\eta }+\mathrm{\epsilon }}{2}$, $\mathrm{\zeta }=\frac{\mathrm{\epsilon }-\mathrm{\eta }}{2\sqrt{ln2}}$, the maximum amplitude is the growth rate ν_SP, the target firing rate ε is the right intersection with the x-axis and the minimum firing rate η is the left intersection with the x-axis. In our application there are no negative firing rates and η is omitted from Figure 1D as we are only interested in the curve on the positive side of the x-axis. Internally, each neuron has a count of the number of connected synaptic elements (z_connected). When z is higher than z_connected, new connections can be made and when z is lower than z_connected, some connections must be deleted. This will be done until z = z_connected for each update in connectivity. For an analysis of the dynamical stability caused by structural changes induced by structural plasticity in other network models we refer to Nowke et al. (2018).

Spike time-dependent plasticity and structural plasticity act on very different time scales. Let us consider the time scale if only structural plasticity is turned ON, whereas STDP is turned OFF. The growth rate of the structural plasticity mechanism used in our simulations reflects a timescale of network formation which corresponds to a timescale of days or weeks (Butz and van Ooyen, 2013; Butz et al., 2014). Note, arbitrary units were introduced in the structural plasticity model mechanism because there is no exact experimental information about the speed of creation or deletion of synapses in the brain regions of interest available (Butz and van Ooyen, 2013; Butz et al., 2014). There are experimental studies showing that the time it takes for boutons and spines to form ranges from days to weeks, as shown, e.g., in the context of early development (Kirov et al., 2004). The structural plasticity model is not dimensionless in time (Butz and van Ooyen, 2013; Butz et al., 2014). However, due to the lack of detailed experimental data the growth rate is not adapted to experimentally pre-defined values. Rather, the growth rate in this model of structural plasticity can be used to change the speed at which effects in the structure of the network are visible and it is chosen so that the total duration of the simulations, with time-steps in milliseconds, do not require to sum up to weeks, otherwise they would be computationally intractable. In our simulations we have shown 20 min of simulation which do not represent 20 min of biological time but a timescale of days or weeks. However, the massive reduction of computation time for the structural plasticity-induced dynamics is only possible if STDP is turned OFF while structural plasticity is turned ON (see below). Due to the pronounced structural changes induced in the network by SP, it is necessary to split the simulation and keep the STDP part separated from the SP part, in order to ensure numerical accuracy in the results.

Coordinated Reset Stimulation

In order to model the invasive DBS electrical CR stimulation delivered to the STN population, we consider short biphasic current-controlled pulses as previously used in Ebert et al. (2014) [see Cogan (2008) for original motivation]. In more detail, a current pulse P(t) for stimulation of neuronal tissue consists of a cathodal and an anodal phase with current amplitudes and durations that result in an overall zero net charge for the entire pulse (to guarantee charge-balance):

$$P\left(t\right)=\{\begin{array}{c}\mathrm{\kappa },t_l\le t<t_l+\mathrm{\omega }\hfill \\ -\mathrm{\kappa }/p_s,t_l+\mathrm{\omega }\le t<t_l+\mathrm{\omega }(1+p_s).\hfill \\ 0,\text{else}\hfill \end{array}(20)$$

t_l denotes the onset times of the current pulses, κ is the amplitude and ω is the width of the cathodal pulse. p_s determines the duration and amplitude of the charge-balancing anodal pulse part which prevents any permanent charge transfer into the neuronal tissue that could possibly damage the tissue. In this study, we used a fixed electrode position at the center of the STN population which was found to be most adequate for optimal CR stimulation performance in Ebert et al. (2014). In order to model the electric field produced by the external stimulus, we used the setup of a Medtronic DBS lead model 3389 which has four separate cylindrical contacts made of a Pt-Ir alloy with a typical length of 1.5 mm (Coffey, 2009). We used a specific equation which adequately approximates the overall distance dependent decay of stimulation strength as used in Ebert et al. (2014; see also Richardson et al., 2003):

$$S\left(d_{il}\right)=\frac{1}{d_{il}{l}_c\sqrt{1+4{(d_{il}/l_c)}^2}},(21)$$

Where d_il is the distance between neuron i and the location of the stimulation contact l and l_c is the length of the electrode contacts. By setting d_min = 0.7 mm we avoid any singularity that could occur at d_il = 0. Hence, we prohibit possible neuron coordinates within a cylindrical volume with 1.4 mm diameter around the electrode axis.

We used M =  4 stimulation sites while all neurons receive input from all four sites. Figure 2 shows an example of our CR stimulation signal. During each stimulation period T = 125ms, each site is activated via an electric burst only once, and not at the exact same onset time with any other site. This order of activation varies randomly in every cycle (during the CR ON period). Three cycles of CR stimulation (ON-cycles) were followed by two cycles without stimulation (OFF-cycles) following Lysyansky et al. (2011) and Adamchic et al. (2014a). The stimulation period within a burst is set at T_p = 7.69ms, the CR stimulation amplitude κ = −3.3 is, the width of the cathodal pulse ω= 200μs and the duration and amplitude of the charge-balancing anodal pulse p_s = 8. The CR stimulation signal delivered at neuron i from stimulation contact l is given by the following equation:

FIGURE 2

Figure 2. Temporal evolution of coordinated reset (CR) simulation signals. Successive 3 ON and 2 OFF cycles are indicated by vertical dashed lines and the temporal distance between two successive vertical lines correspond to the period Tof each cycle. Each stimulation site is activated exactly once during an ON cycle, where the sequence of activation across stimulation sites is random.

$$I_{\text{stim}}\left(t\right)=\sum _{l=1}^{M}S\left(d_{il}\right){\mathrm{\rho }}_l\left(t\right)P\left(t\right),(22)$$

ρ_l(t) is an indicator function representing the spatiotemporal activation of stimulation sites (set at 1 when the site is activated and 0 otherwise). For more details on the CR stimulation setup see Ebert et al. (2014).

Macroscopic Measurement

In order to measure the degree of neuronal synchronization within the STN and GPe ensembles we use the order parameter (Kuramoto, 1984; Daido, 1992; Tass, 1999):

$$R\left(t\right)=\left|N^{-1}\sum _j{\text{e}}^{i{\mathrm{\phi }}_j\left(t\right)}\right|,(23)$$

Where i here denotes the unit imaginary unit $\sqrt{-1}$, φ_j(t) =  2π(t−t_{j, m})/(t_{j, m + 1}−t_{j, m})for t_{j, m}≤t < t_{j, m + 1} is a linear approximation of the phase of neuron j between its m^th and (m + 1)^th spikes at spiking times t_j,m and t_j,m+1. R(t) is influenced by the synaptic weights, as the latter are time dependent due to the STDP. The order parameter R measures the extent of phase synchronization in the neuronal ensemble and takes values between 0 (absence of in-phase synchronization) and 1 (perfect in-phase synchronization).

Results

Simulation Description and Protocols

We present the general simulation setups and protocol designs we employed throughout this study.

• We start in a desynchronized state with sufficiently weak synapses and large numbers of possible synaptic connections within the STN population (i.e., large number of dendritic spines) and corresponding firing rate of ∼4 Hz.

• To generate a synchronized reference state, we deliver sufficiently strong kindling-like PS to increase synaptic weights, so that the network is shifted to a strongly synaptically connected and strongly synchronized stable state. In several model networks kindling is achieved more rapidly than anti-kindling (Tass and Majtanik, 2006). The kindling process is achieved by PS of comparably short duration (T_PS =  2.5min).

• The PS-induced synchronized state is stable, so that both STN and GPe remain in a stable synchronized state after cessation of PS.

• CR stimulation is then delivered to the network in the synchronized state, which causes a desynchronization and a reduction of its synchrony (as measured by the order parameter).

• If CR stimulation is delivered for a sufficiently long time, during the CR-off period the network relaxes to a desynchronized state with down-regulated synaptic weights but still a large number of possible synaptic connections (i.e., large number of dendritic spines).

• Up to this point STDP is active for STN neurons.

• From this point on, two options are being explored. The system may either evolve only with STDP activated (Figure 3 and Supplementary Figure 2) or, alternatively, with structural plasticity activated between CR intervals (Figures 4, 5 and Supplementary Figures 1, 3, 4). In the latter case, we consider different scenarios that evolve depending on whether the target firing rate of the structural plasticity mechanism is modified by CR stimulation.

FIGURE 3

Figure 3. Periodic Stimulation (PS)-CR sequences without structural plasticity. (A) Time evolution of the order parameter R(t)averaged over a sliding window (10 ms) for STN (red solid line) and GPe (blue solid line) neurons. (B) Spiking rates for STN and GPe neurons for three identical CR periods. Light blue bands denote the PS (2.5 min) intervals and light gray the CR intervals (10 min). STDP is active throughout the entire trial. The synchronization-desynchronization process caused by the first PS-CR sequence is repeated by the two subsequent PS-CR sequences, indicating its reproducibility in the absence of PS. (C) PS-CR sequence with one long CR epoch followed by two PS-CR sequences with insufficient CR duration (6 and 5 min, respectively) without structural plasticity.

FIGURE 4

Figure 4. Three identical PS-CR sequences, with structural plasticity with 3.0 Hz target firing rate. (A) A sliding window of 100 simulation time steps at each different epoch when STDP or structural plasticity was active for STN (red solid line) and GPe (blue solid line) neurons and the average value for the whole network (cyan solid line). The horizontal back dot-dashed lines indicate the plateaus’ gradual decrease of the level of synchronization degree after each SP period interval. The vertical black dot-dashed depict the gradual improvement of the CR’ performance after each consecutive stimulation and SP cycle. The lines’ respective location is set at the beginning of the plateaus immediately after the initial sharp decrease. The equally sized horizontal red lines and the circle are used to provide a visual aid to this effect. (B) Spiking rates for STN and GPe neurons, respectively. PS, CR, STDP, and structural plasticity activation interval period. The color coding and time units are defined in the legends of panels (A,C). (C) Time course of the total number of synaptic excitatory connections in the STN population.

FIGURE 5

Figure 5. Periodic Stimulation-Coordinated Reset sequence with long CR epoch followed by two PS-CR sequences with reduced CR duration, with structural plasticity (SP) with 3.0 Hz target firing rate. (A) A sliding window of 100 simulation time steps at each different epoch when STDP or structural plasticity was active for STN (red solid line) and GPe (blue solid line) neurons and the average value for the whole network (cyan solid line). (B) Spiking rates for STN and GPe neurons, respectively. PS, CR, STDP, and structural plasticity activation interval period color coding and units are defined in Figure 3.

This elementary stimulation protocol comprises a PS-induced transition from a desynchronized to a synchronized state and a CR-induced reverse transition. In this PS-CR sequence, PS serves as a test stimulation assessing the network’s susceptibility to synchronizing stimulation. We use different combinations and variations of this elementary stimulation protocol. First, to demonstrate reversibility, we simply repeat the same PS-CR sequence again using the same CR-ON time period and the exact same number of connections. Second, after the first PS-CR sequence, we turn on structural plasticity in the desynchronized state for a computational time interval corresponding to arbitrary time units. In our implementation, STDP and structural plasticity operate at different time scales simulating biological time which ranges from milliseconds up to minutes when STDP is ON and up to several days when structural plasticity is ON. It is important to note that due to software constraints only one of the two plasticity mechanisms is allowed to be ON in a given time interval. Otherwise, the growth rate in the structural plasticity model could not be used to change the speed at which structural plasticity effects occur, which would render the model computationally intractable. In this way, the simulation does not mix the different physiological components and effects, taking place on different time scales. The structural plasticity model efficiently simulates relatively large biological time intervals in feasible computing times. Furthermore, it allows the system to evolve smoothly enough, maintaining its dynamical stability during the synapse rewiring process. The simulated structural plasticity period of approximately 20 min computing time corresponds to a biological time that ranges from days to weeks (Butz et al., 2014). This enables a slow dynamical evolution of the network’s topology, governed by the structural plasticity homeostatic rule with pre-defined or stimulation-induced target parameters related to non-pathological neuronal spiking rates. The subsequent second PS-CR sequence acts like a probe, detecting PS-induced alterations of the network’s synchronizability and desynchronizability.

In this study, we configured the structural plasticity algorithm as follows: we use Gaussian growth curves and set ν_SP = 0.00008 synaptic elements/ms for the STN-STN connections, ν_SP = 0.00002 synaptic elements/ms for the STN-GPe connections and ε = 3.0Hz for both pre-synaptic elements. We note that the units of the growth rate ν_SP are defined here in synaptic elements/ms of simulation time. Biologically realistic values of the growth rate ν_SP would be at least two orders of magnitude slower than the ones used in this manuscript. The use of these higher values allows us to have faster changes in the connectivity while keeping the normal simulation step of NEST and reduce the simulation time representing the structural changes which would be observed in extended biological times. This setup is meant to drive the target STN firing rate toward a non-pathological rate of about 3.0 Hz. In our implementation, structural plasticity is applied only to excitatory connections originating from the STN population (see “Materials and Methods” section for a detailed description). The growth rate has positive values because the synapses are excitatory. The curve specifies that when neurons have a firing rate higher than 3.0 Hz, synapses will be deleted; if neurons have a firing rate lower than 3.0 Hz, synapses will be created. Rewiring occurs when a synaptic element is deleted together with its corresponding synapse, but the matching synaptic element in the other neuron remains available for a new connection with a new matching element. We provide a rewiring factor in the form of Gaussian curves for the post-synaptic elements of both connections of ν_SP = −0.005 synaptic elements/ms. This relatively higher value of ν_SP compared to the ones in the pre-synaptic growth curves allows the formation of several potential new contact points to which deleted connections can be rewired. In our study, we observed a firing rate of approximately 3.6 Hz at the moment structural plasticity is enabled after each PS-CR sequence. Note, the firing rate is the STN network’s actual firing rate which is measured by the structural plasticity algorithm and taken into account for the generation/deletion of synaptic elements, based on the structural plasticity homeostatic rule. The chosen configuration assumes that as a result of a CR epoch the target firing rate of the STN population is reset to 3.0 Hz, which reflects a “healthy” firing rate. This means that during the simulation, synapses will slowly be deleted until the desired 3.0 Hz target is reached. We have experimented with a broad range of parameter combinations before choosing these values. The results discussed here use a configuration within a range of parameters which show memory effects, see Supplementary Material for more details. The major rationale was to design a structural plasticity mechanism whose time evolution causes sufficient synaptic deletion in order to impact on the network’s macroscopic dynamics (as measured by the order parameter) within computationally feasible times.

As already mentioned, the structural plasticity algorithm measures and takes into account the current STN firing rate and modifies the synaptic connectivity until the target (3.0 Hz) is reached. With the selected parameters, on each period of simulation with structural plasticity approximately 5.33% of the total STN-STN connections and 4.8% of the total STN-GPe connections are effectively deleted. We then apply a second, identical PS-CR sequence or a modified PS-CR sequence. To this end, we either keep the same CR-ON period or gradually shorten it to 6 and 5 min. After a sufficient amount of time STDP does no longer change the order parameter and mean firing rate. Based on a series of simulations it turned out that a 10 min post-CR window with STDP ON is sufficient for the network to attain a stable state as characterized by stable order parameter and firing rate. We then turn on structural plasticity and turn off STDP, since an increased growth rate is required to speed structural plasticity simulations up to achieve feasible computation times.

Repeated Delivery of PS-CR Sequences With STPD Only

Anti-kindling requires CR stimulation to be delivered at sufficient amplitude and duration (Tass and Majtanik, 2006; Hauptmann and Tass, 2009; Popovych and Tass, 2012; Ebert et al., 2014). Figure 3A shows the overall desynchronization effect for sufficiently strong CR stimulation amplitude κ=−3.3 (see Ebert et al., 2014 for details). We synchronize and desynchronize the network repetitively by delivering three PS-CR sequences (see “Results” Section: Simulation Description and Protocols). The network starts in a desynchronized state, and then gets synchronized by PS. PS epochs are highlighted by light blue color bands, ranging from $t_{\text{start}}^{PS}=\mathrm{\hspace{0.33em}2}$ min to $t_{\text{stop}}^{PS}=\mathrm{\hspace{0.33em}4.5}$min. After PS delivery the network first evolves spontaneously; that is, without further stimulation. Next, the PS-CR sequence is completed by administration of CR stimulation. CR epochs are indicated by light gray color bands, ranging from $t_{\text{start}}^{CR}=\mathrm{\hspace{0.33em}7}$ min to $t_{\text{stop}}^{CR}=17$min. A 10 min CR epoch is considered a “long” CR period. Subsequently, the network again evolves spontaneously, thereby remaining in a desynchronized state. Afterward, two more PS-CR sequences are delivered. Note, at the beginning of each PS-CR sequence, we reset the time back to zero, while the simulation proceeds in a continuous manner. Resetting the time scale is performed in order to facilitate successive comparisons with altered plasticity mechanisms and/or CR stimulation durations. In Figure 3B, we show the corresponding time evolution of the spiking rates of the STN (red curve) and GPe (blue curve) neurons, starting at ∼4 Hz/∼8 Hz (healthy model state), increasing to ∼15 Hz/∼21 Hz (pathological model state) and up to ∼23 Hz/∼34 Hz (during the CR epoch) and decreasing back to ∼4 Hz/∼8 Hz (after the PS-CR sequence, in the healthy model state). This process is repeated for each of the three PS-CR sequences in a reproducible manner, based on the macroscopic quantities used for assessment of overall synchrony and mean spiking rate. In the Supporting Information Section, we present a similar simulation with three PS-CR sequences with insufficiently strong CR stimulation amplitude (κ = −2.0) (Supplementary Figure 2); despite appropriate stimulation duration, the administered CR stimulation is not sufficient to desynchronize the STN-GPe circuit. The transient behavior of all macrovariables, STN R(t), GPe R(t) as well as STN and GPe mean firing rates, during the CR epoch depends on model parameters, e.g., on the plasticity ratio between depression and potentiationγ -value and CR stimulation intensity κ -value. For instance, for smaller γ -values (e.g., γ≲1.3 and fixed κ = −3.3) and identical CR stimulation duration there is no re-increase of the order parameter of GPe.

In general, in stimulated neural networks with STDP the values of macrovariables, such as order parameter R(t), mean synaptic weight, mean firing rate, may substantially depend on whether stimulation is ON and specifics of stimulus parameters and protocols (Manos et al., 2018a,b). In the PS epoch vs. the pause between the PS and the CR epoch the system approaches stable states that differ with respect to their order parameters and mean firing rates. In the model under consideration, depending on model and stimulus parameters different transients may occur. We do not have a theory or mechanistic explanation for the transitory decrease and re-increase of GPe R(t) observed for the parameters under consideration.

Next, we consider the impact of insufficiently short CR epochs (Figure 3C). To this end, the first PS-CR sequence with sufficient CR duration (10 min) is followed by two PS-CR sequences with insufficient CR duration (5 and 6 min, respectively). The short CR epoch does not cause a full-blown desynchronization. Rather, it shifts the network to a stable state with somewhat reduced synchronization. The effects of the second PS-CR sequence are reproduced by the third PS-CR sequence. This illustrates that the STN-GPe network is multistable. Apart from a fully synchronized state and a fully desynchronized state, we observe stable states with intermediate levels of synchronization. The short CR epochs only shift the network to such an intermediate state, not to a fully desynchronized state.

Impact of Structural Plasticity on Three PS-CR Sequences

To study whether structural plasticity effects may accumulate despite interspersed PS-CR sequences, we administered three identical PS-CR sequences with long CR epochs and keep the target firing rate at 3.0 Hz (Figures 4A,B). Here, we also present the average (combined STN and GPe populations) order parameter value for the whole network is shown by the cyan solid line. The horizontal black dot-dashed lines indicate the plateaus’ gradual decrease of the level of synchronization degree after each structural plasticity period interval. In addition, lower levels of the order parameters and mean firing rates are reached faster during the third CR epoch compared to the second CR epoch, indicating cumulative structural plasticity effects. The vertical black dot-dashed depict the gradual improvement of the CR’ performance after each consecutive stimulation and SP cycle. The lines’ respective location is set at the beginning of the plateaus immediately after the initial sharp decrease. The equally sized horizontal red lines and the circle are used to provide a visual aid to this effect. The time course of the total number of synaptic excitatory connections in the STN population is shown in Figure 4C.

The cumulative effect of structural plasticity is further illustrated by delivering one PS-CR sequence with long CR epoch (10 min) followed by two PS-CR sequences with shorter CR duration (5 and 6 min) (Figure 5). We observe the same stepwise decrease of the post-PS STN synchronization (Figure 4A) and firing rates (Figure 5B). During the third CR epoch STN and GPe synchronization reaches lower levels and approaches steady state levels more rapidly compared to during the second CR epoch (Figure 5A).

Several model parameters are drawn from Gaussian distributions (see Supplementary Material). To check whether the model’s dynamics is robust with respect to variations of these parameters, for all relevant findings we have run several simulations, even if we presented only one. In fact, a relevant limitation of the structural plasticity model mechanism used in our manuscript is its computational cost which is even demanding for the JURECA supercomputer center at Jülich Research Center. To illustrate the robustness of our findings, we ran 20 random initializations for the simulation with structural plasticity and one PS-CR period (see Supplementary Figure 4).

Structural plasticity with sufficiently small target firing rate reduces the maximum values of STN synchronization. This is because the maximally achievable synchronization, as assessed by the order parameter, depends on the number of excitatory synapses. The more excitatory synapses can attain maximum connectivity values the greater the order parameter, see, e.g., Taylor (2012) and Townsend et al. (2020) and references therein. For illustration, let us consider the effect of a plain removal of a fraction of randomly selected STN-STN connections after initializing the network in a synchronized state (Supplementary Figure 5). To this end, we initialized our model in a synchronized state by delivering PS, with STDP turned ON and structural plasticity permanently turned OFF. We then turn off STDP, randomly deleted between 5 and 20% percentage of the excitatory STN-STN connections and evaluated the amount of synchronization in a 2 min time window. With increasing percentage of deleted STN-STN connections the amount of STN synchronization continuously decreases (Supplementary Figure 5).

Discussion

This is the first computational study investigating the impact of structural plasticity on the outcome of CR stimulation. We implemented structural plasticity in a STN-GPe network with plastic STN-STN synapses governed by STDP. In the absence of structural plasticity, repeated administration of PS-CR sequences leads to reversible transitions between desynchronized and synchronized states (Figures 3A,B). By the same token, delivering PS-CR sequences with CR epochs of insufficient duration causes transitions between synchronized states of different amount of synchrony (Figure 3C and Supplementary Figure 2). To date, changes of susceptibility to CR stimulation after stimulation-free pauses have not been tested in computational studies. However, such memory-like phenomena were previously observed in a clinical study with vibrotactile CR stimulation (Pfeifer et al., 2021) and acoustic CR stimulation (Tass et al., 2012a). In general, effects of this kind may allow the further reduction of the total stimulation duration when designing therapeutic treatment regimes. In an example such as DBS, this could further reduce the occurrence or severity of AE in treated patients, which could significantly improve their quality of life. It is not straightforward to relate our measures, such as the order parameter, to clinical outcome scores. However, it was shown that CR-DBS-induced changes of UPDRS motor score (normalized to baseline) were significantly correlated with changes of peak beta power of the local field potential (LFP) (normalized to baseline) measured by the very STN contacts used for stimulation (Adamchic et al., 2014a). This correlation was obtained off medication, i.e., after proper withdrawal of Parkinson’s medication. Furthermore, Pfeifer et al. (2021) revealed that after 3 months of vCR therapy the off medication high beta band power in the sensorimotor cortex decreased significantly. Hence, after 3 months of vCR the sensorimotor cortex’ ability to resynchronize after profound medication withdrawal was limited. The 3-month interval corresponds to clinical findings obtained in a study with pre-defined timing of clinical visits (Pfeifer et al., 2021) and, hence, can only provide a first time scale approximation. The approximative 3-month time scale is also in agreement with findings obtained in a clinical proof of concept study with acoustic CR stimulation for the treatment of tinnitus (Tass et al., 2012a).

Main Findings and Conclusion

The goal of our study was to demonstrate that clinically observed phenomena can be reproduced by incorporating structural plasticity as opposed to taking into account STDP alone. Our model qualitatively reproduces previously observed memory-like phenomena (Tass et al., 2012a; Pfeifer et al., 2021) if we assume that CR stimulation epochs cause a reduction of the target firing rate. Comparable memory effects of previous treatment sessions come into play when structural plasticity is activated with appropriate target firing rates (Figures 4, 5). On one hand, the PS epochs in the PS-CR sequences may serve as a standardized model process accounting for detrimental influences on patients during stimulation-free intervals, leading to a re-increase of symptoms (see Tass et al., 2012a). On the other hand, the PS epochs may be considered as standardized probes, enabling the testing of the network’s resistance to synchronizing stimulation protocols. Furthermore, we showed that simulation setups which account only for STDP (Figure 3) are not able to exhibit this type of cumulative effect observed in clinical trials, acting on time scales of months (see, e.g., Figure 8 in Pfeifer et al., 2021).

In this study we have taken into account STDP and structural plasticity by two separate mechanisms. STDP changes synaptic weights according to the timing relationship of the corresponding pre- and postsynaptic neurons on a fast time scale (Gerstner et al., 1996; Markram et al., 1997; Bi and Poo, 1998, 2001; Song et al., 2000; Ebert et al., 2014), whereas structural plasticity adapts the single neuron’s firing rate to a target firing rate by creating and deleting synapses on a slow time scale (Butz and van Ooyen, 2013; Diaz-Pier et al., 2016), as discussed above. So far, in other modeling studies certain types of structural plasticity mechanisms and STDP were simultaneously implemented at a more theoretical level, making connections to Hamiltonian sampling (Yu et al., 2016) or Bayesian inference mechanisms (Kappel et al., 2015). In contrast, in our study, we separated the two mechanisms for computational convenience, efficiency and performance, since in a first approximation we assume that no new synapses will be created or pre-existing synapses will be deleted during the short long CR epochs. In contrast, we simulate the slow structural plasticity effects initiated by the brief CR epochs. However, in future studies firing pattern-dependent mechanisms such as STDP may also be included in the structural plasticity model in order to account for selective, firing pattern-guided creation or deletion of synapses (Clopath et al., 2010; Legenstein and Maass, 2011; Perin et al., 2011; Zheng et al., 2013). The actual structural plasticity model will certainly impact on the design of optimal stimulation patterns which cause long-lasting desynchronization by therapeutic synaptic rewiring, achieved by adequately controlling collective firing patterns as well as firing rates. Note, in our model we incorporated plasticity mechanisms, STDP and structural plasticity, in the directly stimulated target area, the STN, only.

The numerically implemented SP mechanism in our simulations provides a plausible explanation of the physiological memory-like effects. However, this approach requires an appropriate parameter tuning in order to demonstrate such effects within computationally reasonable time scales. Assuming that the target firing rate is sufficiently low (3.0 Hz), characterizing a healthy STN firing state, structural plasticity continuously down-regulates the firing rate by deleting excitatory synapses. Consequently, during the PS epoch the STN re-synchronization is less pronounced and does not reach the initial levels observed by the first PS epoch (Figure 4A). The level of synchronization of the GPe network, which provides an inhibitory output, during and directly after the PS epoch, does not change substantially compared to the first PS epoch (Figure 4A). However, during the subsequent CR epoch the desynchronization of STN and GPe are both accelerated in a similar manner compared to the first CR epoch (Figure 4A). Accordingly, a long-lasting desynchronization of STN and GPe can be achieved by means of a considerably shorter CR epoch (Figure 4A). An important prerequisite of the structural plasticity effect is that after the CR epoch the STN firing rate relaxes to a value (3.6 Hz) sufficiently close to the structural plasticity’s target firing rate (3.0 Hz). This enables structural plasticity to take over efficiently and down-regulate the number of excitatory connections within the STN, together with the firing rate using the homeostatic rule. The choice of a Gaussian curve in the structural plasticity model employed in this study (Bos et al., 2015; Diaz-Pier et al., 2016) was motivated by experimental data (Lipton and Kater, 1989; Al-Mohanna et al., 1992) which suggests that the generation of new spines starts slowly when the activity in the neuron is low and progressively increases to a maximum. Then it starts to decay until a homeostatic equilibrium is reached (Butz and van Ooyen, 2013). During development the degree of plasticity during critical periods also increases slowly until a maximum level, and then decreases until a degree of maturity is achieved (see box 1 in Hensch, 2005). Note, a reduction of approximately 5.76% of STN-STN connections together with the rewiring of connections during the structural plasticity epoch are sufficient to increase the network’s resistance to PS and also its increased susceptibility to CR stimulation; thus, the structural plasticity effects accumulate, even though intersected by a PS-CR sequence (Figure 5).

In this study with our novel NEST based code, we combine the power of structural and synaptic plasticity providing a modeling test environment that enables us to capture memory-like phenomena of susceptibility to CR stimulation. From the modeling and simulation perspectives, there are several challenges related to this setup; the most significant is that simulations of detailed networks are computationally demanding, especially when the two types of plasticity are involved. This is mainly due to the variations in the transmission and generation of spiking events in the network, as well as the additional calculations that take place at the synapse level. Simulations of the size that we report in this study can have a simulation to real time factor of about 10–20. This means that simulations of 1 s take 10–20 s to be simulated computationally. Accordingly, exploring slow structural plasticity induced changes of a simulated network within 1 month of biological time, could take a normal simulation of up to several (∼20) months duration. Obviously, such a scenario is significantly resource-intensive. In contrast, using our NEST simulations, we can accelerate the effects of structural plasticity so that simulations can be reduced from months to a couple of hours. With this solution, we are able to simulate at feasible computing times both the synaptic plasticity, which takes place during PS and CR epochs, followed by a stabilization period and the subsequent activation of structural plasticity, running at reasonably high speed by means of an accelerated model. During this accelerated simulation period, we disable synaptic plasticity in order to avoid inconsistencies that would be induced due to the fact that the model operates at different time scales. By combining both types of plasticity in such a temporally interleaved way, we are able to explore the immediate stimulation-induced effects of synaptic plasticity as well as the long-term effects mediated by structural plasticity within a single consistent, efficient, and tractable simulation model.

Future studies could be designed to investigate the development of desynchronizing stimulation protocols that are adapted to specific ranges of target firing rates. It is therefore crucial to understand stimulation-induced changes of the firing rate. In fact, with very few exceptions (see, e.g., Lysyansky et al., 2011), the vast majority of computational studies devoted to CR stimulation to date have not focused on stimulation-induced changes of the firing rate. Our study highlights that not only the CR epochs may be crucial for the long-term stimulation outcome; to the best of our knowledge, our study is the first to also demonstrate that pauses in between effective CR epochs also play an active and critical role. CR stimulation and, in general, all kinds of efficient desynchronizing stimulation protocols may initiate STN-GPe circuits in favorable states, so that structural plasticity makes networks more susceptible to desynchronizing stimulation over time. In a previous computational study in a neural network with STDP, but no structural plasticity, it was shown that spaced CR stimulation (i.e., CR stimulation intermingled with sufficiently long stimulation-free pauses) may significantly improve efficacy; this is particularly important if CR stimulation is delivered at otherwise ineffectively weak intensities (Popovych et al., 2015). Together with our current study, this illustrates the importance and active role of pauses in therapeutic processes employing plasticity principles. Stimulus-free intervals should not just be considered periods without intervention. Rather they may be an integral part of a therapeutic process enabling and potentiating stimulus effects which should be adequately addressed by dose-finding studies. In future studies we plan to incorporate plasticity mechanisms in additional brain areas. However, given the huge amount of computation time required to networks with STDP and structural plasticity, one might modify and further adapt structural plasticity model mechanisms to given high-performance computing architectures.

In summary, in our proof of concept-like computational study we added structural plasticity to an existing STN-GPe model which, so far, only comprised STDP in the STN (Ebert et al., 2014). The STN is an excitatory population (Shink et al., 1996). To model STDP for STN’s efferent excitatory synapses, we have used a standard STDP model which has been used in numerous other studies (see, e.g., Gerstner et al., 1996; Gerstner and Kistler, 2002; Dan and Poo, 2004; Sjöström et al., 2008). In future studies both structural plasticity and STDP should also be taken into account for the GPe. Unlike STN, GPe is an inhibitory population (Shink et al., 1996). The structural plasticity mechanism used in this manuscript (Butz and van Ooyen, 2013; Butz et al., 2014) can also be applied to inhibitory synapses. However, to model STDP in inhibitory synapses, appropriate rules of rules of inhibitory STDP have to be taken into account (Vogels et al., 2013).

Code and Scripts Availability

The NEST version used for the simulations presented in this manuscript can be found here: https://github.com/sdiazpier/nest-simulator.git (branch dbs_sp) and the model scripts can be found in this repository: https://github.com/sdiazpier/dbs_sp.git.

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

Author Contributions

TM, SD-P, and PT conceived and designed the experiments, analyzed the data, and wrote the manuscript. TM and SD-P performed the experiments and contributed reagents, materials, and analysis tools. All authors contributed to the article and approved the submitted version.

Funding

This research was partially funded by the Helmholtz Association through the Helmholtz Portfolio Theme “Supercomputing and Modeling for the Human Brain.” This project was also received funding from the European Union’s Horizon 2020 Research and Innovation Program under grant agreement nos. 720270 (HBP SGA1) and 785907 (HBP SGA2). TM was also supported by the Labex MME DII (ANR-11-LBX-0023-01) French national funding program.

Conflict of Interest

PT works as consultant for Boston Scientific Neuromodulation and Gretap AG and is inventor or co-inventor on a number of patents for invasive and non-invasive neuromodulation.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

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.

Acknowledgments

We are grateful to the computing time granted by the JARA-HPC Vergabegremium and provided on the JARA-HPC Partition part of the supercomputer JURECA at Forschungszentrum Jülich (VSR computation time grant JINB33). We gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this project by providing computing time through the John von Neumann Institute for Computing (NIC) on the GCS Supercomputer JUWELS at Jülich Supercomputing Centre (JSC). We would like to thank Hans E. Plesser, Abigail Morisson, Alex Peyser, Jochen Martin Eppler, and Dimitri Plotnikov for many fruitful discussions, for helping us adjust our code and model modules in the recent distribution of NEST as well as in optimizing it for Forschungszentrum Jülich JURECA High Performance Computing Architectures. TM would also like to thank the entire NEST community for providing solutions for every problem encountered. PT gratefully acknowledges support by the John A. Blume Foundation and by the Binns Family Foundation. We gratefully acknowledge editorial support by Christine Plant (Stanford Neurosurgery).

Supplementary Material

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

References

Adamchic, I., Hauptmann, C., Barnikol, U. B., Pawelczyk, N., Popovych, O., Barnikol, T. T., et al. (2014a). 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

Adamchic, I., Hauptmann, C., and Tass, P. A. (2012). Changes of oscillatory activity in pitch processing network and related tinnitus relief induced by acoustic CR neuromodulation. Front. Syst. Neurosci. 6:18. doi: 10.3389/fnsys.2012.00018

PubMed Abstract | CrossRef Full Text | Google Scholar

Adamchic, I., Langguth, B., Hauptmann, C., and Tass, P. A. (2014b). Abnormal cross-frequency coupling in the tinnitus network. Front. Neurosci. 8:284. doi: 10.3389/fnins.2014.00284

PubMed Abstract | CrossRef Full Text | Google Scholar

Adamchic, I., Toth, T., Hauptmann, C., and Tass, P. A. (2013). Reversing pathologically increased EEG power by acoustic coordinated reset neuromodulation. Hum. Brain Mapp. 35, 2099–2118. doi: 10.1002/hbm.22314

PubMed Abstract | CrossRef Full Text | Google Scholar

Al-Mohanna, F. A., Cave, J., and Bolsover, S. R. (1992). A narrow window of intracellular calcium concentration is optimal for neurite outgrowth in rat sensory neurons. Dev. Brain Res. 70, 287–290. doi: 10.1016/0165-3806(92)90209-F

CrossRef Full Text | Google Scholar

Baizabal-Carvallo, J. F., and Jankovic, J. (2016). Movement disorders induced by deep brain stimulation. Parkinsonism Relat. Disord. 25, 1–9. doi: 10.1016/j.parkreldis.2016.01.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Baufreton, J., Kirkham, E., Atherton, J. F., Menard, A., Magill, P. J., Bolam, J. P., et al. (2009). Sparse but selective and potent synaptic transmission from the globus pallidus to the subthalamic nucleus. J. Neurophysiol. 102, 532–545. doi: 10.1152/jn.00305.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Benabid, A. L., Chabardes, S., Mitrofanis, J., and Pollak, P. (2009). Deep brain stimulation of the subthalamic nucleus for the treatment of Parkinson’s disease. Lancet Neurol. 8, 67–81. doi: 10.1016/S1474-4422(08)70291-6

CrossRef Full Text | Google Scholar

Benabid, A. L., Pollak, P., Gervason, C., Hoffmann, D., Gao, D. M., Hommel, M., et al. (1991). Long-term suppression of tremor by chronic stimulation of the ventral intermediate thalamic nucleus. Lancet 337, 403–406. doi: 10.1016/0140-6736(91)91175-t

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

Bi, G. Q., and Poo, M. M. (2001). Synaptic modification by correlated activity: Hebb’s postulate revisited. Annu. Rev. Neurosci. 24, 139–166. doi: 10.1146/annurev.neuro.24.1.139

PubMed Abstract | CrossRef Full Text | Google Scholar

Blundell, I., Plotnikov, D., Eppler, J. M., and Morrison, A. (2018). Automatically selecting a suitable integration scheme for systems of differential equations in neuron models. Front. Neuroinform. 12:50. doi: 10.3389/fninf.2018.00050

PubMed Abstract | CrossRef Full Text | Google Scholar

Bos, H., Morrison, A., Peyser, A., Hahne, J., Helias, M., and Kunkel, S. (2015). NEST 2.10.0. Available online at: http://dx.doi.org/10.5281/zenodo.44222 (accessed December 31, 2015).

Google Scholar

Brown, C. E., Boyd, J. D., and Murphy, T. H. (2010). Longitudinal in vivo imaging reveals balanced and branch-specific remodeling of mature cortical pyramidal dendritic arbors after stroke. J. Cereb. Blood Flow Metab. 30, 783–791. doi: 10.1038/jcbfm.2009.241

PubMed Abstract | CrossRef Full Text | Google Scholar

Brown, P. (2002). Oscillatory nature of human basal ganglia activity: relationship to the pathophysiology of Parkinson’s disease. Mov. Disord. 18, 357–363. doi: 10.1002/mds.10358

PubMed Abstract | CrossRef Full Text | Google Scholar

Buhmann, C., Huckhagel, T., Engel, K., Gulberti, A., Hidding, U., Poetter-Nerger, M., et al. (2017). Adverse events in deep brain stimulation: a retrospective long-term analysis of neurological, psychiatric and other occurrences. PLoS One 12:e0178984. doi: 10.1371/journal.pone.0178984

PubMed Abstract | CrossRef Full Text | Google Scholar

Butz, M., Steenbuck, I., and van Ooyen, A. (2014). Homeostatic structural plasticity increases the efficiency of small-world networks. Front. Synaptic Neurosci. 6:7. doi: 10.3389/fnsyn.2014.00007

PubMed Abstract | CrossRef Full Text | Google Scholar

Butz, M., and van Ooyen, A. (2013). A simple rule for dendritic spine and axonal bouton formation can account for cortical reorganization after focal retinal lesions. PLoS Comput. Biol. 9:e1003259. doi: 10.1371/journal.pcbi.1003259

PubMed Abstract | CrossRef Full Text | Google Scholar

Butz, M., Wörgötter, F., and van Ooyen, A. (2009). Activity-dependent structural plasticity. Brain Res. Rev. 60, 287–305. doi: 10.1016/j.brainresrev.2008.12.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Cleary, D. R., Ozpinar, A., Raslan, A. M., and Ko, A. L. (2015). Deep brain stimulation for psychiatric disorders: where we are now. Neurosurg. Focus 38:E2. doi: 10.3171/2015.3.FOCUS1546

PubMed Abstract | CrossRef Full Text | Google Scholar

Clopath, C., Büsing, L., Vasilaki, E., and Gerstner, W. (2010). Connectivity reflects coding: a model of voltage-based STDP with homeostasis. Nat. Neurosci. 13, 344–352. doi: 10.1038/nn.2479

PubMed Abstract | CrossRef Full Text | Google Scholar

Coffey, R. J. (2009). Deep brain stimulation devices: a brief technical history and review. Artif. Organs 33, 208–220. doi: 10.1111/j.1525-1594.2008.00620.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Cogan, S. F. (2008). Neural stimulation and recording electrodes. Annu. Rev. Biomed. Eng. 10, 275–309. doi: 10.1146/annurev.bioeng.10.061807.160518

PubMed Abstract | CrossRef Full Text | Google Scholar

Daido, H. (1992). Order function and macroscopic mutual entrainment in uniformly coupled limit-cycle oscillators. Prog. Theor. Phys. Suppl. 88, 1213–1218. doi: 10.1143/ptp/88.6.1213

CrossRef Full Text | Google Scholar

Dan, Y., and Poo, M. M. (2004). Spike timing-dependent plasticity of neural circuits. Neuron 44, 23–30. doi: 10.1016/j.neuron.2004.09.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Dayan, P., and Abbott, L. F. (2000). Theoretical Neuroscience. Cambridge MA: The MIT Press.

Google Scholar

Deller, T., Haas, C. A., Freiman, T. M., Phinney, A., Jucker, M., and Frotscher, M. (2006). “Lesion-induced axonal sprouting in the central nervous system,” in Brain Repair, ed. M. Bähr (Boston, MA: Springer US), 101–121. doi: 10.1007/0-387-30128-3_6

CrossRef Full Text | Google Scholar

Deniau, J.-M., Degos, B., Bosch, C., and Maurice, N. (2010). Deep brain stimulation mechanisms: beyond the concept of local functional inhibition. Eur. J. Neurosci. 32, 1080–1091. doi: 10.1111/j.1460-9568.2010.07413.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Deuschl, G., Schade-Brittinger, C., Krack, P., Volkmann, J., Schäfer, H., Bötzel, K., et al. (2006). A Randomized trial of deep-brain stimulation for Parkinson’s disease. N. Engl. J. Med. 355, 896–908. doi: 10.1056/NEJMoa060281

PubMed Abstract | CrossRef Full Text | Google Scholar

Diaz-Pier, S., Naveau, M., Butz-Ostendorf, M., and Morrison, A. (2016). Automatic generation of connectivity for large-scale neuronal network models through structural plasticity. Front. Neuroanat. 10:57. doi: 10.3389/fnana.2016.00057

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

Fauth, M., and Tetzlaff, C. (2016). Opposing effects of neuronal activity on structural plasticity. Front. Neuroanat. 10:75. doi: 10.3389/fnana.2016.00075

PubMed Abstract | CrossRef Full Text | Google Scholar

Fujimoto, K., and Kita, H. (1993). Response characteristics of subthalamic neurons to the stimulation of the sensorimotor cortex in the rat. Brain Res. 23, 185–192. doi: 10.1016/0006-8993(93)90872-k

CrossRef Full Text | Google Scholar

Gerstner, W., Kempter, R., van Hemmen, J. L., and Wagner, H. (1996). A neuronal learning rule for sub-millisecond temporal coding. Nature 383, 76–78. doi: 10.1038/383076a0

PubMed Abstract | CrossRef Full Text | Google Scholar

Gerstner, W., and Kistler, W. M. (2002). Spiking Neuron Models. Cambridge: Cambridge University Press.

Google Scholar

Gillies, A. J., and Willshaw, D. J. (1998). A massively connected subthalamic nucleus leads to the generation of widespread pulses. Proc. R. Soc. Lond. B 265, 2101–2109. doi: 10.1098/rspb.1998.0546

PubMed Abstract | CrossRef Full Text | Google Scholar

Gradinaru, V., Mogri, M., Thompson, K. R., Henderson, J. M., and Deisseroth, K. (2009). Optical deconstruction of Parkinsonian neural circuitry. Science 324, 354–359. doi: 10.1126/science.1167093

PubMed Abstract | CrossRef Full Text | Google Scholar

Hauptmann, C., and Tass, P. A. (2009). Cumulative and after-effects of short and weak coordinated reset stimulation: a modeling study. J. Neural Eng. 6:016004. doi: 10.1088/1741-2560/6/1/016004

CrossRef Full Text | Google Scholar

Hensch, T. K. (2005). Critical period plasticity in local cortical circuits. Nat. Rev. Neurosci. 6, 877–888. doi: 10.1038/nrn1787

PubMed Abstract | CrossRef Full Text | Google Scholar

Holgado, A. J. N., Terry, J. R., and Bogacz, R. (2010). Conditions for the generation of beta oscillations in the subthalamic nucleus-globus pallidus network. J. Neurosci. 30, 12340–12352. doi: 10.1523/jneurosci.0817-10.2010

PubMed Abstract | CrossRef Full Text | Google Scholar

Johnson, M. D., Miocinovic, S., McIntyre, C. C., and Vitek, J. L. (2008). Mechanisms and targets of deep brain stimulation in movement disorders. Neurotherapeutics 5, 294–308. doi: 10.1016/j.nurt.2008.01.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Kappel, D., Habenschuss, S., Legenstein, R., and Maass, W. (2015). Network plasticity as Bayesian inference. PLoS Comput. Biol. 11:e1004485. doi: 10.1371/journal.pcbi.1004485

PubMed Abstract | CrossRef Full Text | Google Scholar

Keck, T., Mrsic-Flogel, T. D., Vaz Afonso, M., Eysel, U. T., Bonhoeffer, T., and Hübener, M. (2008). Massive restructuring of neuronal circuits during functional reorganization of adult visual cortex. Nat. Neurosci. 11, 1162–1167. doi: 10.1038/nn.2181

PubMed Abstract | CrossRef Full Text | Google Scholar

Kirov, S. A., Goddard, C. A., and Harris, K. M. (2004). Age-dependence in the homeostatic upregulation of hippocampal dendritic spine number during blocked synaptic transmission. Neuropharmacology 47, 640–648. doi: 10.1016/j.neuropharm.2004.07.039

PubMed Abstract | CrossRef Full Text | Google Scholar

Kita, H., Tachibana, Y., Nambu, A., and Chiken, S. (2005). Balance of monosynaptic excitatory and disynaptic inhibitory responses of the globus pallidus induced after stimulation of the subthalamic nucleus in the monkey. J. Neurosci. 25, 8611–8619. doi: 10.1523/jneurosci.1719-05.2005

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., and Tass, P. A. (2020). Long-lasting desynchronization by decoupling stimulation. Phys. Rev. Res. 2:033101.

Google Scholar

Kühn, A. A., Kempf, F., Brücke, C., Gaynor Doyle, L., Martinez-Torres, I., Pogosyan, A., et al. (2008). High-frequency stimulation of the subthalamic nucleus suppresses oscillatory β activity in patients with Parkinson’s disease in parallel with improvement in motor performance. J. Neurosci. 28, 6165–6173. doi: 10.1523/jneurosci.0282-08.2008

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuramoto, Y. (1984). Chemical Oscillations, Waves, and Turbulence. Berlin: Springer.

Google Scholar

Legenstein, R., and Maass, W. (2011). Branch-specific plasticity enables self-organization of nonlinear computation in single neurons. J. Neurosci. 31, 10787–10802. doi: 10.1523/jneurosci.5684-10.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

Lenz, F. A., Kwan, H. C., Martin, R. L., Tasker, R. R., Dostrovsky, J. O., and Lenz, Y. E. (1994). Single unit analysis of the human ventral thalamic nuclear group – tremor-related activity in functionally identified cells. Brain 117, 531–543. doi: 10.1093/brain/117.3.531

PubMed Abstract | CrossRef Full Text | Google Scholar

Levesque, 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

Lipton, S. A., and Kater, S. B. (1989). Neurotransmitter regulation of neuronal outgrowth, plasticity and survival. Trends Neurosci. 12, 265–270. doi: 10.1016/0166-2236(89)90026-X

CrossRef Full Text | Google Scholar

Lisman, J. (2017). Glutamatergic synapses are structurally and biochemically complex because of multiple plasticity processes: long-term potentiation, long-term depression, short-term potentiation and scaling. Philos. Trans. R. Soc. B Biol. Sci. 372:20160260. doi: 10.1098/rstb.2016.0260

PubMed Abstract | CrossRef Full Text | Google Scholar

Llinás, R. R., Ribary, U., Jeanmonod, D., Kronberg, E., and Mitra, P. P. (1999). Thalamocortical dysrhythmia: a neurological and neuropsychiatric syndrome characterized by magnetoencephalography. Proc. Natl. Acad. Sci. U.S.A. 96, 15222–15227. doi: 10.1073/pnas.96.26.15222

PubMed Abstract | CrossRef Full Text | Google Scholar

Lysyansky, B., Popovych, O. V., and Tass, P. A. (2011). Desynchronizing anti-resonance effect of m:n ON–OFF coordinated reset stimulation. J. Neural Eng. 8:036019. doi: 10.1088/1741-2560/8/3/036019

CrossRef Full Text | Google Scholar

Manos, T., Zeitler, M., and Tass, P. A. (2018a). 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

Manos, T., Zeitler, M., and Tass, P. A. (2018b). Short-term dosage regimen for stimulation-induced long-lasting desynchronization. Front. Physiol. 9:376. doi: 10.3389/fphys.2018.00376

PubMed Abstract | CrossRef Full Text | Google Scholar

Markram, H., Lübke, J., Frotscher, M., and Sakmann, B. (1997). Regulation of synaptic efficacy by coincidence of postsynaptic APs and EPSPs. Science 275, 213–215. doi: 10.1126/science.275.5297.213

PubMed Abstract | CrossRef Full Text | Google Scholar

Milosevic, L., Kalia, S. K., Hodaie, M., Lozano, A. M., Fasano, A., Popovic, M. R., et al. (2018). Neuronal inhibition and synaptic plasticity of basal ganglia neurons in Parkinson’s disease. Brain 141, 177–190. doi: 10.1093/brain/awx296

PubMed Abstract | CrossRef Full Text | Google Scholar

Nini, A., Feingold, A., Slovin, H., and Bergman, H. (1995). Neurons in the globus pallidus do not show correlated activity in the normal monkey, but phase-locked oscillations appear in the MPTP model of Parkinsonism. J. Neurophysiol. 74, 1800–1805. doi: 10.1152/jn.1995.74.4.1800

PubMed Abstract | CrossRef Full Text | Google Scholar

Nowke, C., Diaz-Pier, S., Weyers, B., Hentschel, B., Morrison, A., Kuhlen, T. W., et al. (2018). Toward rigorous parameterization of underconstrained neural network models through interactive visualization and steering of connectivity generation. Front. Neuroinform. 12:32. doi: 10.3389/fninf.2018.00032

PubMed Abstract | CrossRef Full Text | Google Scholar

Perin, R., Berger, T. K., and Markram, H. (2011). A synaptic organizing principle for cortical neuronal groups. Proc. Natl. Acad. Sci. U.S.A. 108, 5419–5424. doi: 10.1073/pnas.1016051108

PubMed Abstract | CrossRef Full Text | Google Scholar

Perun, K., Rumpe, B., Plotnikov, D., Trensch, G., Eppler, J. M., Blundell, I., et al. (2018). Reengineering NestML with Python and MontiCore. Master Thesis. Aachen: RWTH Aachen University. doi: 10.5281/zenodo.1319653

CrossRef Full Text | Google Scholar

Pfeifer, K. J., Kromer, J. A., Cook, A. J., Hornbeck, T., Lim, E. A., Mortimer, B., 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

Plotnikov, D., Blundell, I., Ippen, T., Eppler, J. M., Morrison, A., and Rumpe, B. (2016). “NESTML: a modeling language for spiking neurons,” in Proceedings of the Modellierung 2016, eds A. Oberweis, and R. Reussner (Bonn: Gesellschaft für Informatik e.V.).

Google Scholar

Popovych, O., 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

Popovych, O. V., Xenakis, M. N., and Tass, P. A. (2015). The spacing principle for unlearning abnormal neuronal synchrony. PLoS One 10:e0117205. doi: 10.1371/journal.pone.0117205

PubMed Abstract | CrossRef Full Text | Google Scholar

Prescott, I. A., Dostrovsky, J. O., Moro, E., Hodaie, M., Lozano, A. M., and Hutchison, W. D. (2009). Levodopa enhances synaptic plasticity in the substantia nigra pars reticulata of Parkinson’s disease patients. Brain 132, 309–318. doi: 10.1093/brain/awn322

PubMed Abstract | CrossRef Full Text | Google Scholar

Reetz, K., Gaser, C., Klein, C., Hagenah, J., Büchel, C., Gottschalk, S., et al. (2009). Structural findings in the basal ganglia in genetically determined and idiopathic Parkinson’s disease. Mov. Disord. 24, 99–103. doi: 10.1002/mds.22333

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

Rubin, J. E., and Terman, D. (2004). High frequency stimulation of the subthalamic nucleus eliminates pathological thalamic rhythmicity in a computational model. J. Comput. Neurosci. 16, 211–235. doi: 10.1023/b:jcns.0000025686.47117.67

CrossRef Full Text | Google Scholar

Sehm, B., Taubert, M., Conde, V., Weise, D., Classen, J., Dukart, J., et al. (2014). Structural brain plasticity in Parkinson’s disease induced by balance training. Neurobiol. Aging 35, 232–239. doi: 10.1016/j.neurobiolaging.2013.06.021

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

Shink, E., Bevan, M. D., Bolam, J. P., and Smith, Y. (1996). The subthalamic nucleus and the external pallidum: two tightly interconnected structures that control the output of the basal ganglia in the monkey. Neuroscience 73, 335–357. doi: 10.1016/0306-4522(96)00022-x

CrossRef Full Text | Google Scholar

Silchenko, A. N., Adamchic, I., Hauptmann, C., and Tass, P. A. (2013). Impact of acoustic coordinated reset neuromodulation on effective connectivity in a neural network of phantom sound. Neuroimage 77, 133–147. doi: 10.1016/j.neuroimage.2013.03.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Sjöström, P. J., Rancz, E. A., Roth, A., and Häusser, M. (2008). Dendritic excitability and synaptic plasticity. Physiol. Rev. 88, 769–840. doi: 10.1152/physrev.00016.2007

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

Steiner, L. A., Barreda Tomás, F. J., Planert, H., Alle, H., Vida, I., and Geiger, J. R. P. (2019). Connectivity and dynamics underlying synaptic control of the subthalamic nucleus. J. Neurosci. 39, 2470–2481. doi: 10.1523/JNEUROSCI.1642-18.2019

PubMed Abstract | CrossRef Full Text | Google Scholar

Syrkin-Nikolau, J., Neuville, R., O’Day, J., Anidi, C., Miller Koop, M., Martin, T., et al. (2018). Coordinated reset vibrotactile stimulation shows prolonged improvement in Parkinson’s disease. Mov. Disord. 33, 179–180. doi: 10.1002/mds.27223

PubMed Abstract | CrossRef Full Text | Google Scholar

Tass, P. A. (1999). Phase Resetting in Medicine and Biology – Stochastic Modelling and Data Analysis. Berlin: Springer.

Google Scholar

Tass, P. A. (2003a). 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. (2003b). Desynchronization by means of a coordinated reset of neural sub-populations. Prog. Theor. Phys. Suppl. 150, 281–296. doi: 10.1143/ptps.150.281

PubMed Abstract | CrossRef Full Text | Google Scholar

Tass, P. A. (2017). Vibrotactile coordinated reset stimulation for the treatment of neurological diseases: concepts and device specifications. Cureus 9:e1535. doi: 10.7759/cureus.1535

PubMed Abstract | CrossRef Full Text | Google Scholar

Tass, P. A., Adamchic, I., Freund, H.-J., von Stackelberg, T., and Hauptmann, C. (2012a). Counteracting tinnitus by acoustic coordinated reset neuromodulation. Restor. Neurol. Neurosci. 30, 137–159. doi: 10.3233/rnn-2012-110218

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., and Popovych, O. (2012). Unlearning tinnitus-related cerebral synchrony with acoustic coordinated reset stimulation: theoretical concept and modelling. Biol. Cybern. 106, 27–36. doi: 10.1007/s00422-012-0479-5

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Tass, P. A., Silchenko, A., Hauptmann, C., Barnikol, U. B., and Speckmann, E. J. (2009). Long-lasting desynchronization in rat hippocampal slice induced by coordinated reset stimulation. Phys. Rev. E 80:011902.

Google Scholar

Taylor, R. (2012). There is no non-zero stable fixed point for dense networks in the homogeneous Kuramoto model. J. Phys. A Math. Theor. 45:055102. doi: 10.1088/1751-8113/45/5/055102

CrossRef Full Text | Google Scholar

Temperli, P., Ghika, J., Villemure, J. G., Burkhard, P. R., Bogousslavsky, J., and Vingerhoets, F. J. G. (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

Terman, D., Rubin, J. E., Yew, A. C., and Wilson, C. J. (2002). Activity patterns in a model for the subthalamopallidal network of the basal ganglia. J. Neurosci. 22, 2963–2976. doi: 10.1523/jneurosci.22-07-02963.2002

PubMed Abstract | CrossRef Full Text | Google Scholar

Townsend, A., Stillman, M., and Strogatz, S. H. (2020). Dense networks that do not synchronize and sparse ones that do. Chaos 30:083142. doi: 10.1063/5.0018322

CrossRef Full Text | Google Scholar

van Hartevelt, T. J., Cabral, J., Deco, G., Møller, A., Green, A. L., Aziz, T. Z., et al. (2014). Neural plasticity in human brain connectivity: the effects of long term deep brain stimulation of the subthalamic nucleus in Parkinson’s disease. PLoS One 9:e86496. doi: 10.1371/journal.pone.0086496

PubMed Abstract | CrossRef Full Text | Google Scholar

van Hartevelt, T. J., Cabral, J., Møller, A., FitzGerald, J. J., Green, A. L., Aziz, T. Z., et al. (2015). Evidence from a rare case study for Hebbian-like changes in structural connectivity induced by long-term deep brain stimulation. Front. Behav. Neurosci. 9:167. doi: 10.3389/fnbeh.2015.00167

PubMed Abstract | CrossRef Full Text | Google Scholar

van Ooyen, A., and Butz-Ostendorf, M. (2017). The Rewiring Brain A Computational Approach to Structural Plasticity in the Adult Brain. London: Academic Press.

Google Scholar

Vogels, T., Froemke, R., Doyon, N., Gilson, M., Haas, J., Liu, R., et al. (2013). Inhibitory synaptic plasticity: spike timing-dependence and putative network function. Front. Neural Circuits 7:119. doi: 10.3389/fncir.2013.00119

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

Xu, T., Yu, X., Perlik, A. J., Tobin, W. F., Zweig, J. A., Tennant, K., et al. (2009). Rapid formation and selective stabilization of synapses for enduring motor memories. Nature 462, 915–919. doi: 10.1038/nature08389

PubMed Abstract | CrossRef Full Text | Google Scholar

Yamahachi, H., Marik, S. A., McManus, J. N. J., Denk, W., and Gilbert, C. D. (2009). Rapid axonal sprouting and pruning accompany functional reorganization in primary visual cortex. Neuron 64, 719–729. doi: 10.1016/j.neuron.2009.11.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, Z., Kappel, D., Legenstein, R., Song, S., Chen, F., and Maass, W. (2016). CaMKII Activation Supports Reward-Based Neural Network Optimization Through Hamiltonian Sampling (arXiv:1606.00157v3).

Google Scholar

Zheng, P., Dimitrakakis, C., and Triesch, J. (2013). Network self-organization explains the statistics and dynamics of synaptic connection strengths in cortex. PLoS Comput. Biol. 9:e1002848. doi: 10.1371/journal.pcbi.1002848

PubMed Abstract | CrossRef Full Text | Google Scholar