Dr Thanos Manos
email: thanos.manos@cyu.fr

Sandra Diaz-Pier

Peter A Tass

Long-term desynchronization by coordinated reset stimulation in a neural network model with synaptic and structural plasticity

Keywords: Coordinated reset neuromodulation, Desynchronization, Spike time-dependent plasticity (STDP), Structural plasticity, Anti-kindling

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 stimulationfree 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 stimulationand 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) [START_REF] Benabid | Long-term suppression of tremor by chronic stimulation of the ventral intermediate thalamic nucleus[END_REF][START_REF] Deuschl | A Randomized Trial of Deep-Brain Stimulation for Parkinson's Disease[END_REF][START_REF] Krack | Five-Year Follow-up of Bilateral Stimulation of the Subthalamic Nucleus in Advanced Parkinson's Disease[END_REF] and is also being tested in other psychiatric diseases [START_REF] Cleary | Deep brain stimulation for psychiatric disorders: where we are now[END_REF]. 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 [START_REF] Benabid | Deep brain stimulation of the subthalamic nucleus for the treatment of Parkinson's disease[END_REF][START_REF] Benabid | Long-term suppression of tremor by chronic stimulation of the ventral intermediate thalamic nucleus[END_REF]. However, the precise mechanism of action of standard HF DBS is not yet sufficiently understood [START_REF] Deniau | Deep brain stimulation mechanisms: beyond the concept of local functional inhibition[END_REF][START_REF] Gradinaru | Optical Deconstruction of Parkinsonian Neural Circuitry[END_REF][START_REF] Johnson | Mechanisms and targets of deep brain stimulation in movement disorders[END_REF]. Standard HF DBS has acute clinical [START_REF] Temperli | How do parkinsonian signs return after discontinuation of subthalamic DBS?[END_REF] and electrophysiological [START_REF] Kühn | High-Frequency Stimulation of the Subthalamic Nucleus Suppresses Oscillatory β Activity in Patients with Parkinson's Disease in Parallel with Improvement in Motor Performance[END_REF] 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 & Jankovic, 2017;[START_REF] Buhmann | Adverse events in deep brain stimulation: A retrospective long-term analysis of neurological, psychiatric and other occurrences[END_REF]. Therefore, any DBS treatment regime should be minimized in duration as much as possible, in order to also minimize the AEs that can result [START_REF] Buhmann | Adverse events in deep brain stimulation: A retrospective long-term analysis of neurological, psychiatric and other occurrences[END_REF].

Abnormal neuronal synchronization is a hallmark of Parkinson's disease (PD) and may occur in different frequency ranges [START_REF] Brown | Oscillatory nature of human basal ganglia activity: Relationship to the pathophysiology of Parkinson's disease[END_REF][START_REF] Lenz | Single unit analysis of the human ventral thalamic nuclear group -Tremor-related activity in functionally identified cells[END_REF][START_REF] Llinás | Thalamocortical dysrhythmia: A neurological and neuropsychiatric syndrome characterized by magnetoencephalography[END_REF][START_REF] Nini | 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[END_REF]. Different stimulation protocols showed that trains of HF DBS may in fact modulate synaptic plasticity [START_REF] Milosevic | Neuronal inhibition and synaptic plasticity of basal ganglia neurons in Parkinson's disease[END_REF][START_REF] Prescott | Levodopa enhances synaptic plasticity in the substantia nigra pars reticulata of Parkinson's disease patients[END_REF][START_REF] Shen | Synaptic plasticity in rat subthalamic nucleus induced by high-frequency stimulation[END_REF]. 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. [START_REF] Lisman | Glutamatergic synapses are structurally and biochemically complex because of multiple plasticity processes: long-term potentiation, long-term depression, short-term potentiation and scaling[END_REF]. STDP is a fundamental synaptic plasticity mechanism which modifies synaptic strengths based on the relative timing of pre-and postsynaptic spike pairs [START_REF] Bi | Synaptic modifications in cultured hippocampal neurons: Dependence on spike timing, synaptic strength, and postsynaptic cell type[END_REF][START_REF] Bi | Synaptic modification by correlated activity: Hebb's postulate revisited[END_REF][START_REF] Gerstner | A neuronal learning rule for sub-millisecond temporal coding[END_REF][START_REF] Markram | Regulation of synaptic efficacy by coincidence of postsynaptic APs and EPSPs[END_REF][START_REF] Song | Competitive Hebbian learning through spike-timingdependent synaptic plasticity[END_REF]. 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 [START_REF] Markram | Regulation of synaptic efficacy by coincidence of postsynaptic APs and EPSPs[END_REF].

Coordinated reset (CR) stimulation is a computationally developed patterned multichannel stimulation (Tass, 2003;Tass, 2003) that employs STDP to induce long-lasting therapeutic effects [START_REF] Tass | Long-term anti-kindling effects of desynchronizing brain stimulation: a theoretical study[END_REF]. CR stimulation specifically counteracts abnormal neuronal synchrony by desynchronization (Tass, 2003;Tass, 2003), which causes a reduction of the rate of coincidences and, mediated by STDP, a reduction of the synaptic weights [START_REF] Tass | Long-term anti-kindling effects of desynchronizing brain stimulation: a theoretical study[END_REF]. 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) [START_REF] Tass | Long-term anti-kindling effects of desynchronizing brain stimulation: a theoretical study[END_REF]. Accordingly, in neuronal network models a longlasting desynchronization can be achieved. The term long-lasting desynchronization refers to desynchronization which exceeds cessation of stimulation (see Manos et al., 2018a;[START_REF] Tass | Long-term anti-kindling effects of desynchronizing brain stimulation: a theoretical study[END_REF][START_REF] Kromer | Long-lasting desynchronization by decoupling stimulation[END_REF]. Such effects have been demonstrated pre-clinically and clinically in a variety of studies e.g. [START_REF] Tass | Long-lasting desynchronization in rat hippocampal slice induced by coordinated reset stimulation[END_REF], Adamchic et al., 2012;[START_REF] Adamchic | Reversing pathologically increased EEG power by acoustic coordinated reset neuromodulation[END_REF]Tass et al., 2012a;Tass et al., 2012b[START_REF] Pfeifer | Coordinated Reset Vibrotactile Stimulation Induces Sustained Cumulative Benefits in Parkinson's Disease[END_REF]. For CR-DBS delivered to the subthalamic nucleus (STN), computationally predicted desynchronizing effects (Tass, 2003;Tass, 2003), cumulative effects (Hauptmann & Tass, 2009) and long-lasting effects [START_REF] Tass | Long-term anti-kindling effects of desynchronizing brain stimulation: a theoretical study[END_REF] were verified in preclinical studies in a 1-methyl-4-phenyl-1,2,3,6-tetrahydropyridine (MPTP) monkey PD model (Tass et al., 2012;[START_REF] Wang | 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[END_REF] as well as in a proof-of-concept study in PD patients (Adamchic et al., 2014a).

In accordance with computational predictions (Popovych & Tass, 2012;Tass & 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., 2012) combined with a significant reduction of abnormal neuronal synchrony [START_REF] Adamchic | Changes of oscillatory activity in pitch processing network and related tinnitus relief induced by acoustic CR neuromodulation[END_REF][START_REF] Adamchic | Reversing pathologically increased EEG power by acoustic coordinated reset neuromodulation[END_REF]Tass et al., 2012), abnormal cross-frequency coupling [START_REF] Adamchic | Abnormal cross-frequency coupling in the tinnitus network[END_REF] and abnormal effective connectivity [START_REF] Silchenko | Impact of acoustic coordinated reset neuromodulation on effective connectivity in a neural network of phantom sound[END_REF]. In a study of Parkinson's disease patients, vibrotactile CR (vCR) stimulation [START_REF] Tass | Vibrotactile Coordinated Reset Stimulation for the Treatment of Neurological Diseases: Concepts and Device Specifications[END_REF] delivered to the fingertips during three consecutive days for four hours per day caused improvements in gait and bradykinesia for up to one-month poststimulation [START_REF] Syrkin-Nikolau | Coordinated reset vibrotactile stimulation shows prolonged improvement in Parkinson's disease[END_REF]. 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) [START_REF] Pfeifer | Coordinated Reset Vibrotactile Stimulation Induces Sustained Cumulative Benefits in Parkinson's Disease[END_REF]. Regular vCR uses periodic stimulus delivery, whereas for noisy vCR stimulus onsets had a moderate jitter in time [START_REF] Pfeifer | Coordinated Reset Vibrotactile Stimulation Induces Sustained Cumulative Benefits in Parkinson's Disease[END_REF]. In addition, in a pilot study six PD patients received daily noisy vCR stimulation during three months [START_REF] Pfeifer | Coordinated Reset Vibrotactile Stimulation Induces Sustained Cumulative Benefits in Parkinson's Disease[END_REF]. 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;[START_REF] Pfeifer | Coordinated Reset Vibrotactile Stimulation Induces Sustained Cumulative Benefits in Parkinson's Disease[END_REF] 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 [START_REF] Fauth | Opposing Effects of Neuronal Activity on Structural Plasticity[END_REF], 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 preexisting synaptic connections; this occurs via the extension or retraction of axons and dendrites, or by modifying the number of axonal boutons or dendritic spines [START_REF] Butz | Activity-dependent structural plasticity[END_REF]. Changes in structural connectivity and rewiring of connections is crucial for memory formation [START_REF] Xu | Rapid formation and selective stabilization of synapses for enduring motor memories[END_REF] as well as for adaptive and maladaptive processes in response to central and peripheral lesions of the nervous system [START_REF] Brown | Longitudinal in vivo imaging reveals balanced and branch-specific remodeling of mature cortical pyramidal dendritic arbors after stroke[END_REF][START_REF] Keck | Massive restructuring of neuronal circuits during functional reorganization of adult visual cortex[END_REF][START_REF] Yamahachi | Rapid Axonal Sprouting and Pruning Accompany Functional Reorganization in Primary Visual Cortex[END_REF] and neurodegeneration [START_REF] Deller | Lesion-Induced Axonal Sprouting in the Central Nervous System[END_REF]. A study of the grey 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 [START_REF] Reetz | Structural findings in the basal ganglia in genetically determined and idiopathic Parkinson's disease[END_REF]. Furthermore, it was shown that balance training induces morphometric changes in PD patients [START_REF] Sehm | Structural brain plasticity in Parkinson's disease induced by balance training[END_REF]: 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 [START_REF] Sehm | Structural brain plasticity in Parkinson's disease induced by balance training[END_REF]. In a single PD patient case study, it was shown that long-term HF DBS gave rise to significant localized structural changes in sensorymotor, prefrontal, limbic, and olfactory brain regions together with a restoration of functional connectivity [START_REF] Van Hartevelt | Neural plasticity in human brain connectivity: the effects of long term deep brain stimulation of the subthalamic nucleus in Parkinson's disease[END_REF]. In addition, a computational analysis revealed that the observed changes of structural weights are compliant with a hypothesized Hebbian-like mechanism [START_REF] Van Hartevelt | Evidence from a rare case study for Hebbian-like changes in structural connectivity by long-term deep brain stimulation[END_REF]. 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 & [START_REF] Van Ooyen | The Rewiring Brain A Computational Approach to Structural Plasticity in the Adult Brain[END_REF].

Accordingly, we set out to computationally integrate a structural plasticity mechanism into our modelling approach; thus, in addition to STDP we also used the built-in structural plasticity module [START_REF] Diaz-Pier | Automatic Generation of Connectivity for Large-Scale Neuronal Network Models through Structural Plasticity[END_REF] 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 [START_REF] Butz | A Simple Rule for Dendritic Spine and Axonal Bouton Formation Can Account for Cortical Reorganization after Focal Retinal Lesions[END_REF] and for the details on the exact implementation see [START_REF] Diaz-Pier | Automatic Generation of Connectivity for Large-Scale Neuronal Network Models through Structural Plasticity[END_REF]. 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;[START_REF] Pfeifer | Coordinated Reset Vibrotactile Stimulation Induces Sustained Cumulative Benefits in Parkinson's Disease[END_REF], 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 [START_REF] Rubin | High frequency stimulation of the subthalamic nucleus eliminates pathological thalamic rhythmicity in a computational model[END_REF][START_REF] Terman | Activity Patterns in a Model for the Subthalamopallidal Network of the Basal Ganglia[END_REF] 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:

𝑐 m 𝑑𝑣 𝑑𝑡 = -𝐼 L -𝐼 K -𝐼 Na -𝐼 T -𝐼 Ca -𝐼 ahp -𝐼 syn + 𝐼 stim + 𝐼 noise (1)
Spiking activity is caused by the sodium (Na + ) and potassium (K + ) ionic currents 𝐼 Na , 𝐼 K . 𝐼 T and 𝐼 Ca describe the low-threshold T-type and high-threshold Ca 2+ current respectively. 𝐼 ahp represents a Ca 2+ -activated, voltage-independent after-hyperpolarization K + current and 𝐼 L the leak current. In addition, the STN and GPe neurons are influenced by synaptic inputs 𝐼 syn . Surrounding brain areas contribute with 𝐼 noise while there is an external stimulation current 𝐼 stim (only in the STN population) which models the deep brain stimulation (DBS). Other ionic currents (in pA/μm 2 ) are described by [see [START_REF] Ebert | Coordinated reset stimulation in a large-scale model of the STN-GPe circuit[END_REF][START_REF] Rubin | High frequency stimulation of the subthalamic nucleus eliminates pathological thalamic rhythmicity in a computational model[END_REF][START_REF] Terman | Activity Patterns in a Model for the Subthalamopallidal Network of the Basal Ganglia[END_REF]]:

𝐼 L = 𝑔 L [𝑣 -𝑣 L ] (2) 𝐼 K = 𝑔 K 𝑛 4 [𝑣 -𝑣 K ] (3)
𝐼 𝑻 = 𝑔 T 𝑎 ∞ 3 (𝑣)𝑏 ∞ 2 (𝑟)[𝑣 -𝑣 Ca ] (4) 𝐼 Ca = 𝑔 Ca 𝑠 ∞ 2 (𝑣)ℎ[𝑣 -𝑣 Ca ] (5) 𝐼 ahp = 𝑔 ahp [𝑣 -𝑣 K ] [Ca] [Ca] + 𝑘 1 (6) 𝑑[Ca] 𝑑𝑡 = Γ(-𝐼 Ca -𝐼 T -𝑘 Ca [Ca]) (7)
STN and GPe neurons are described by similar equations and they differ only in a few parameter values (see Supplementary Tables 1 and2) as well as in the form of the low threshold Ca 2+ current:

𝐼 T = 𝑔 T 𝑎 ∞ 3 (𝑣) 𝑟 [𝑣 -𝑣 Ca ]
(8) for the GPe neurons (i.e., the 𝑏 ∞ 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 𝑣 L where we draw values with 0.15% standard deviation). All the mean parameter values and their units are summarized in Supplementary Tables 1 and2. The slowly operating gating variables 𝑛, ℎ, and 𝑟 depend on time and voltage and their first-order kinetics are governed by differential equations of the form:

𝑑𝑋 𝑑𝑡 = 𝜑 𝑋 [ 𝑋 ∞ (𝑣) -𝑋 𝜏 𝑋 (𝑣) ] , (9)
where 𝑋: 𝑛, ℎ, 𝑟 with:

𝜏 𝑋 (𝑣) = 𝜏 𝑋 0 + 𝜏 𝑋 1 1 + exp [- 𝑣 -𝜃 𝑋 𝜏 𝑜 𝑋 𝜏 ] , (10)
where 𝜃 𝜏 is the voltage at which the time constant is midway between its maximum and minimum values, and 𝜊 𝜏 is the slope factor for the voltage dependence of the time constant. Activation gating for the rapidly activating channels (𝑚, 𝑎, and 𝑠) was treated as instantaneous. For all gating variables 𝑋 = 𝑛, 𝑚, ℎ, 𝑎, 𝑟, or 𝑠, the steady-state voltage dependence was determined using:

𝑋 ∞ (𝑣) = 1 1 + exp [- 𝑣 -𝜃 𝑋 𝜊 𝑋 ] , ( 11
)
where 𝜃 𝑋 is the half activation (or inactivation) voltage for gating variable 𝑋, and 𝑜 𝑋 is the slope factor for that variable. For the 𝑇 current inactivation variable 𝑏, we used:

𝑏 ∞ (𝑟) = 1 1 + exp [ 𝑟 -𝜃 𝑏 𝜊 𝑏 ] - 1 1 + exp [- 𝜃 𝑏 𝜊 𝑏 ] , (12)
[see [START_REF] Terman | Activity Patterns in a Model for the Subthalamopallidal Network of the Basal Ganglia[END_REF] and [START_REF] Rubin | High frequency stimulation of the subthalamic nucleus eliminates pathological thalamic rhythmicity in a computational model[END_REF] for more information on the model and the variables' units]. One main difference of our study compared to [START_REF] Ebert | Coordinated reset stimulation in a large-scale model of the STN-GPe circuit[END_REF] is the initialization of the STN synaptic weights. In [START_REF] Ebert | Coordinated reset stimulation in a large-scale model of the STN-GPe circuit[END_REF], the initial mean synaptic weight was set at 𝑤 ̅ ss STN =0.018, in order to tune the system at a strongly synchronized state. Here, we chose a different approach. We set 𝑤 ̅ ss 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 [START_REF] Blundell | Automatically selecting a suitable integration scheme for systems of differential equations in neuron models[END_REF][START_REF] Perun | Reengineering NestML with Python and MontiCore[END_REF][START_REF] Plotnikov | NESTML: a modeling language for spiking neurons[END_REF].

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:

𝑊 = ( 𝑤 ss 𝑤 𝑔𝑠 𝑤 sg 𝑤 gg ) ( 13
)
denoting the connectivity matrices for the STN-STN (ss), STN-GPe (sg), GPe-STN (gs), and GPe-GPe (gg) sub-networks, respectively. Each 𝑊 𝑖𝑗 value corresponds to the connection strength between neuron 𝑖 and j while self-connections are excluded, i.e. 𝑊 𝑖𝑖 = 0, ∀𝑖. 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 [START_REF] Shink | The subthalamic nucleus and the external pallidum: Two tightly interconnected structures that control the output of the basal ganglia in the monkey[END_REF]. The connections between STN and GPe were established following [START_REF] Baufreton | Sparse but selective and potent synaptic transmission from the globus pallidus to the subthalamic nucleus[END_REF] and using connectivity values of 2% (i.e., 200 synapses per neuron, connecting to 20 randomly chosen targets in the distant population respectively). Following [START_REF] Ebert | Coordinated reset stimulation in a large-scale model of the STN-GPe circuit[END_REF][START_REF] Fujimoto | Response characteristics of subthalamic neurons to the stimulation of the sensorimotor cortex in the rat[END_REF], the transmission delay for connections between the STN and the GPe was set to 𝛿 sg = 𝛿 gs = 4.0 ms. Recent findings on the synaptic organization of the STN population identifying its neurons as parallel processing units [START_REF] Steiner | Connectivity and Dynamics Underlying Synaptic Control of the Subthalamic Nucleus[END_REF] might help to further refine computational models of the basal ganglia. 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 STN neurons, modified by STDP and structural plasticity. The external electrical stimulation delivered to the STN is indicated by the grey 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 PD patient prior to 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 equation ( 19) which represents the change in synaptic elements 𝑑𝑧/𝑑𝑡 as a function of the firing rate FR(𝑡). 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.

Network and Neuron Coupling Description

The network setup and 3D anatomical arrangement of the neurons was introduced in [START_REF] Ebert | Coordinated reset stimulation in a large-scale model of the STN-GPe circuit[END_REF] 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 [START_REF] Levesque | GABAergic interneurons in human subthalamic nucleus[END_REF]]. 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 [START_REF] Ebert | Coordinated reset stimulation in a large-scale model of the STN-GPe circuit[END_REF].

The general connectivity properties are largely based on the values used and presented in [START_REF] Ebert | Coordinated reset stimulation in a large-scale model of the STN-GPe circuit[END_REF], where 10,000 STN and 10,000 GPe neurons were considered [more details regarding the parameter values can also be found in [START_REF] Gillies | A massively connected subthalamic nucleus leads to the generation of widespread pulses[END_REF] and [START_REF] Holgado | Conditions for the generation of beta oscillations in the subthalamic nucleus-globus pallidus network[END_REF]]. 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 𝑝 ss = 70% of the entire STN population (instead of 𝑝 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 𝑝(𝑥) = e -𝑥 𝑐 𝑑 ⁄ with 𝑐 𝑑 = 0.5 and 𝑥 the Euclidean distance between neurons. This probability is a nonnegative 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 𝑝 gg = 10% (instead of 𝑝 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 𝑤 ̅ ss STN =0.0025 with standard deviation 𝜎 ss STN = 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 [START_REF] Baufreton | Sparse but selective and potent synaptic transmission from the globus pallidus to the subthalamic nucleus[END_REF]] 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.0 ms for connections between the STN and the GPe, as published in [START_REF] Fujimoto | Response characteristics of subthalamic neurons to the stimulation of the sensorimotor cortex in the rat[END_REF] and [START_REF] Kita | Balance of monosynaptic excitatory and disynaptic inhibitory responses of the globus pallidus induced after stimulation of the subthalamic nucleus in the monkey[END_REF]. Their corresponding initial synaptic weights are drawn again from a Gaussian probability distribution around a mean value 𝑤 ̅ sg STN = 0.006 and 𝑤 ̅ gs STN = 0.003 with standard deviation 𝜎 sg STN = 0.0003 and 𝜎 sg STN = 0.00015 respectively.

Using the aforementioned described parameter values, we achieved a good agreement between the two networks of different size concerning theirdynamics, 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 [START_REF] Dayan | Theoretical Neuroscience[END_REF][START_REF] Gerstner | Spiking neuron models[END_REF]:

𝛼(𝑡) = 𝑡 -𝑡 𝑘 𝜏 syn 2 𝑒 - 𝑡-𝑡 𝑘 𝜏 syn , 𝑡 𝑘 ≤ 𝑡 < 𝑡 𝑘+1 ( 14
)
where 𝑡 𝑘 denotes the spike time. The total synaptic input current to a postsynaptic neuron 𝑖 received from presynaptic neurons 𝑗 is then given by:

𝐼 syn,𝑖 (𝑡) = ∑ 𝑊 𝑖𝑗 (𝑣 𝑖 (𝑡) -𝑣 syn )𝛼(𝑡), 𝑗 (15)
where 𝑊 𝑖𝑗 is the synaptic weight (coupling strength) between presynaptic neurons j and the postsynaptic neuron 𝑖, and 𝑣 syn is the reversal potential for excitatory or inhibitory ss, sg, gs, and gg connections. Therefore, 𝑣 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 𝐼 app = -7.0 pA in Equation (1) (note that for the excitatory STN neurons 𝐼 app = 0 pA) and external noise 𝐼 noise described by Poissonian spike trains (form other surrounding brain structures) with a frequency of f P GPe = 40 Hz:

𝐼 noise,𝑖 = ∑ 𝑤 noise (𝑣 𝑗,𝑖 (𝑡) -𝑣 noise )𝛼(𝑡) 𝑁 𝑡 𝑗=0 . ( 16
)
𝑁 𝑡 is an integer random variable following the Poisson distribution 𝑃 𝜆 (𝑘) = 𝜆 𝑘 𝑘! 𝑒 -𝜆 , with parameter given by the simulated time interval 𝜆 = 𝑓 𝑃 𝛿, 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 f P GPe = 20 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 [START_REF] Ebert | Coordinated reset stimulation in a large-scale model of the STN-GPe circuit[END_REF] and references therein for motivation and justification). The noise input synaptic weight is set at 𝑤 noise = 0.2 mS/cm 2 while the time constants for the 𝛼-function are set at 𝜏 noise = 1.0 ms and the reversal potential is 𝑣 noise = 0 mV for both STN and GPe neurons.

Spike-Timing-Dependent Plasticity

The synaptic coupling weights (𝑤 𝑖𝑗 ) of the STN neurons evolve dynamically and depend on the time difference (𝛥𝑡 𝑖𝑗 = 𝑡 𝑗 𝑓 -𝑡 𝑖 𝑓 ) between the firing (onset) of the spikes of the post-synaptic neuron and pre-synaptic neuron (denoted by 𝑡 𝑖 𝑓 and 𝑡 𝑗 𝑓 respectively). When 𝛥𝑡 𝑖𝑗 > 0 the rule implements synaptic potentiation due to causal relationships and when 𝛥𝑡 𝑖𝑗 < 0 it implements synaptic depression. We use the following STDP rule as implemented in NEST:

𝛥𝑤 𝑖𝑗 (𝛥𝑡 𝑖𝑗 ) = { 𝜆𝑒 - |𝛥𝑡 𝑖𝑗 | 𝜏 + , 𝛥𝑡 𝑖𝑗 > 0 -𝜆𝛾𝑒 - |𝛥𝑡 𝑖𝑗 | 𝜏 -, 𝛥𝑡 𝑖𝑗 ≤ 0 , ( 17
)
where 𝜏 -= 27.5 ms, 𝜏 + = 12.0 ms, (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 𝑤 𝑖𝑗 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 [START_REF] Ebert | Coordinated reset stimulation in a large-scale model of the STN-GPe circuit[END_REF].

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 [START_REF] Butz | A Simple Rule for Dendritic Spine and Axonal Bouton Formation Can Account for Cortical Reorganization after Focal Retinal Lesions[END_REF] and described in technical detail by [START_REF] Diaz-Pier | Automatic Generation of Connectivity for Large-Scale Neuronal Network Models through Structural Plasticity[END_REF]. 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 [START_REF] Butz | A Simple Rule for Dendritic Spine and Axonal Bouton Formation Can Account for Cortical Reorganization after Focal Retinal Lesions[END_REF][START_REF] Diaz-Pier | Automatic Generation of Connectivity for Large-Scale Neuronal Network Models through Structural Plasticity[END_REF]

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 [START_REF] Diaz-Pier | Automatic Generation of Connectivity for Large-Scale Neuronal Network Models through Structural Plasticity[END_REF], 𝛽 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 [START_REF] Nowke | Toward Rigorous Parameterization of Underconstrained Neural Network Models Through Interactive Visualization and Steering of Connectivity Generation[END_REF]. 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.:

𝑑𝑧 𝑑𝑡 = 𝜈 SP [2𝑒 -( FR(𝑡)-𝜉 𝜁 ) 2 -1], ( 19
)
where 𝑧 is the number of synaptic elements,

𝜉 = 𝜂+𝜀 2 , 𝜁 = 𝜀-𝜂 2√𝑙𝑛2
, the maximum amplitude is the growth rate 𝜈 SP , the target firing rate 𝜀 is the right intersection with the 𝑥 -axis and the minimum firing rate 𝜂 is the left intersection with the 𝑥 -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 𝑥 -axis.

Internally, each neuron has a count of the number of connected synaptic elements (𝑧_connected). When z is higher than 𝑧_connected, new connections can be made and when 𝑧 is lower than 𝑧_connected, some connections must be deleted. This will be done until 𝑧 = 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 [START_REF] Nowke | Toward Rigorous Parameterization of Underconstrained Neural Network Models Through Interactive Visualization and Steering of Connectivity Generation[END_REF].

STDP 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 [START_REF] Butz | Homeostatic structural plasticity increases the efficiency of small-world networks[END_REF][START_REF] Butz | A Simple Rule for Dendritic Spine and Axonal Bouton Formation Can Account for Cortical Reorganization after Focal Retinal Lesions[END_REF]. 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 [START_REF] Butz | Homeostatic structural plasticity increases the efficiency of small-world networks[END_REF][START_REF] Butz | A Simple Rule for Dendritic Spine and Axonal Bouton Formation Can Account for Cortical Reorganization after Focal Retinal Lesions[END_REF]. 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 [START_REF] Kirov | Age-dependence in the homeostatic upregulation of hippocampal dendritic spine number during blocked synaptic transmission[END_REF]. The structural plasticity model is not dimensionless in time [START_REF] Butz | Homeostatic structural plasticity increases the efficiency of small-world networks[END_REF][START_REF] Butz | A Simple Rule for Dendritic Spine and Axonal Bouton Formation Can Account for Cortical Reorganization after Focal Retinal Lesions[END_REF]. 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 minutes of simulation which do not represent 20 minutes 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 [START_REF] Ebert | Coordinated reset stimulation in a large-scale model of the STN-GPe circuit[END_REF] [see [START_REF] Cogan | Neural stimulation and recording electrodes[END_REF] 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):

𝑃(𝑡) = { 𝜅, 𝑡 𝑙 ≤ 𝑡 < 𝑡 𝑙 + 𝜔 -𝜅/𝑝 𝑠 , 𝑡 𝑙 + 𝜔 ≤ 𝑡 < 𝑡 𝑙 + 𝜔(1 + 𝑝 𝑠 ) 0, else . ( 20
)
𝑡 𝑙 denotes the onset times of the current pulses, 𝜅 is the amplitude and 𝜔 is the width of the cathodal pulse. 𝑝 𝑠 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 [START_REF] Ebert | Coordinated reset stimulation in a large-scale model of the STN-GPe circuit[END_REF]. 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 [START_REF] Coffey | Deep brain stimulation devices: A brief technical history and review[END_REF]. We used a specific equation which adequately approximates the overall distance dependent decay of stimulation strength as used in [START_REF] Ebert | Coordinated reset stimulation in a large-scale model of the STN-GPe circuit[END_REF] [see also [START_REF] Richardson | In vivo modulation of hippocampal epileptiform activity with radial electric fields[END_REF]]:

𝑆(𝑑 𝑖𝑙 ) = 1 𝑑 𝑖𝑙 𝑙 𝑐 √1 + 4(𝑑 𝑖𝑙 /𝑙 𝑐 ) 2 , ( 21
)
where 𝑑 𝑖𝑙 is the distance between neuron i and the location of the stimulation contact 𝑙 and 𝑙 𝑐 is the length of the electrode contacts. By setting 𝑑 min = 0.7 mm we avoid any singularity that could occur at 𝑑 𝑖𝑙 = 0. Hence, we prohibit possible neuron coordinates within a cylindrical volume with 1.4 mm diameter around the electrode axis.

We used 𝑀 = 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 𝑇 = 125 ms, 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 [START_REF] Adamchic | Coordinated reset neuromodulation for Parkinson's disease: proof-of-concept study[END_REF][START_REF] Lysyansky | Desynchronizing anti-resonance effect of m:n ON--OFF coordinated reset stimulation[END_REF]. The stimulation period within a burst is set at 𝑇 p = 7.69 ms, 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 𝑝 s = = 8. The CR stimulation signal delivered at neuron i from stimulation contact 𝑙 is given by the following equation:

𝐼 stim (𝑡) = ∑ 𝑆(𝑑 𝑖𝑙 )𝜌 𝑙 (𝑡)𝑃(𝑡) 𝑀 𝑙=1 , ( 22
)
𝜌 𝑙 (𝑡) 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 [START_REF] Ebert | Coordinated reset stimulation in a large-scale model of the STN-GPe circuit[END_REF].

Macroscopic measurement

In order to measure the degree of neuronal synchronization within the STN and GPe ensembles we use the order parameter [START_REF] Daido | Order function and macroscopic mutual entrainment in uniformly coupled limitcycle oscillators[END_REF][START_REF] Kuramoto | Chemical oscillations, waves, and turbulence[END_REF][START_REF] Tass | Phase resetting in medicine and biology -stochastic modelling and data analysis[END_REF]:

𝑅(𝑡) = |𝑁 -1 ∑ 𝑒 𝑖𝜑 𝑗 (𝑡) 𝑗 | , ( 23
)
where 𝑖 here denotes the unit imaginary unit √-1, 𝜑 𝑗 (𝑡) = 2𝜋(𝑡 -𝑡 𝑗,𝑚 ) (𝑡 𝑗,𝑚+1 -𝑡 𝑗,𝑚 ⁄ ) for 𝑡 𝑗,𝑚 ≤ 𝑡 < 𝑡 𝑗,𝑚+1 is a linear approximation of the phase of neuron j between its 𝑚 𝑡ℎ and (𝑚 + 1) 𝑡ℎ spikes at spiking times 𝑡 𝑗,𝑚 and 𝑡 𝑗,𝑚+1 . 𝑅(𝑡) 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 periodic stimulation (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 [START_REF] Tass | Long-term anti-kindling effects of desynchronizing brain stimulation: a theoretical study[END_REF]. The kindling process is achieved by PS of comparably short duration (𝑇 PS = 2.5 min). • 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 and4). 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.

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 minutes computing time corresponds to a biological time that ranges from days to weeks [START_REF] Butz | Homeostatic structural plasticity increases the efficiency of small-world networks[END_REF]. 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.0 Hz 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 towards 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 [START_REF] Ebert | Coordinated reset stimulation in a large-scale model of the STN-GPe circuit[END_REF]Hauptmann & Tass, 2009;Popovych & Tass, 2012;[START_REF] Tass | Long-term anti-kindling effects of desynchronizing brain stimulation: a theoretical study[END_REF]. Figure 3A shows the overall desynchronization effect for sufficiently strong CR stimulation amplitude 𝜅 = -3.3 (see [START_REF] Ebert | Coordinated reset stimulation in a large-scale model of the STN-GPe circuit[END_REF] 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 𝑡 start 𝑃𝑆 = 2 min to 𝑡 stop 𝑃𝑆 = 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 grey color bands, ranging from 𝑡 start 𝐶𝑅 = 7 min to 𝑡 stop 𝐶𝑅 = 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. Afterwards, 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 𝑅(𝑡), GPe 𝑅(𝑡) 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 𝑅(𝑡), 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;Manos et al., 2018b). 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 𝑅(𝑡) 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 (6 and 5 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 (Figure 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 dotdashed 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 (6 and 5 min) (Figure 5). We observe the same stepwise decrease 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. [START_REF] Taylor | There is no non-zero stable fixed point for dense networks in the homogeneous Kuramoto model[END_REF][START_REF] Townsend | Dense networks that do not synchronize and sparse ones that do[END_REF] 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 (Figure 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 [START_REF] Pfeifer | Coordinated Reset Vibrotactile Stimulation Induces Sustained Cumulative Benefits in Parkinson's Disease[END_REF] 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 [START_REF] Adamchic | Coordinated reset neuromodulation for Parkinson's disease: proof-of-concept study[END_REF]. This correlation was obtained off medication, i.e., after proper withdrawal of Parkinson's medication. Furthermore, [START_REF] Pfeifer | Coordinated Reset Vibrotactile Stimulation Induces Sustained Cumulative Benefits in Parkinson's Disease[END_REF] 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 [START_REF] Pfeifer | Coordinated Reset Vibrotactile Stimulation Induces Sustained Cumulative Benefits in Parkinson's Disease[END_REF] 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 conclusions:

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 [START_REF] Pfeifer | Coordinated Reset Vibrotactile Stimulation Induces Sustained Cumulative Benefits in Parkinson's Disease[END_REF]Tass et al., 2012a) 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 and5). 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 [START_REF] Pfeifer | Coordinated Reset Vibrotactile Stimulation Induces Sustained Cumulative Benefits in Parkinson's Disease[END_REF]].

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 [START_REF] Bi | Synaptic modifications in cultured hippocampal neurons: Dependence on spike timing, synaptic strength, and postsynaptic cell type[END_REF][START_REF] Bi | Synaptic modification by correlated activity: Hebb's postulate revisited[END_REF][START_REF] Ebert | Coordinated reset stimulation in a large-scale model of the STN-GPe circuit[END_REF][START_REF] Gerstner | A neuronal learning rule for sub-millisecond temporal coding[END_REF][START_REF] Markram | Regulation of synaptic efficacy by coincidence of postsynaptic APs and EPSPs[END_REF][START_REF] Song | Competitive Hebbian learning through spike-timingdependent synaptic plasticity[END_REF], 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 [START_REF] Butz | A Simple Rule for Dendritic Spine and Axonal Bouton Formation Can Account for Cortical Reorganization after Focal Retinal Lesions[END_REF][START_REF] Diaz-Pier | Automatic Generation of Connectivity for Large-Scale Neuronal Network Models through Structural Plasticity[END_REF], as discussed above. So far, in other modelling studies certain types of structural plasticity mechanisms and STDP were simultaneously implemented at a more theoretical level, making connections to Hamiltonian sampling [START_REF] Yu | Hamiltonian synaptic sampling in a model for reward-gated network plasticity[END_REF] or Bayesian inference mechanisms [START_REF] Kappel | Network Plasticity as Bayesian Inference[END_REF]. 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 [START_REF] Focus ; Clopath | Connectivity reflects coding: a model of voltage-based STDP with homeostasis[END_REF][START_REF] Legenstein | Branch-Specific Plasticity Enables Self-Organization of Nonlinear Computation in Single Neurons[END_REF][START_REF] Perin | A synaptic organizing principle for cortical neuronal groups[END_REF][START_REF] Zheng | Network Self-Organization Explains the Statistics and Dynamics of Synaptic Connection Strengths in Cortex[END_REF]. 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;[START_REF] Diaz-Pier | Automatic Generation of Connectivity for Large-Scale Neuronal Network Models through Structural Plasticity[END_REF] was motivated by experimental data [START_REF] Al-Mohanna | A narrow window of intracellular calcium concentration is optimal for neurite outgrowth in rat sensory neurons[END_REF][START_REF] Lipton | Neurotransmitter regulation of neuronal outgrowth, plasticity and survival[END_REF] 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 [START_REF] Butz | A Simple Rule for Dendritic Spine and Axonal Bouton Formation Can Account for Cortical Reorganization after Focal Retinal Lesions[END_REF]. 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 [START_REF] Hensch | Critical period plasticity in local cortical circuits [Review Article[END_REF]]. 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 using our novel NEST based code, we combine the power of structural and synaptic plasticity providing a modelling test environment that enables us to capture memory-like phenomena of susceptibility to CR stimulation. From the modelling 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 second take 10 to 20 seconds 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. [START_REF] Lysyansky | Desynchronizing anti-resonance effect of m:n ON--OFF coordinated reset stimulation[END_REF]], 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 [START_REF] Popovych | The Spacing Principle for Unlearning Abnormal Neuronal Synchrony[END_REF]. 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 [START_REF] Ebert | Coordinated reset stimulation in a large-scale model of the STN-GPe circuit[END_REF]. The STN is an excitatory population [START_REF] Shink | The subthalamic nucleus and the external pallidum: Two tightly interconnected structures that control the output of the basal ganglia in the monkey[END_REF]. 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. [START_REF] Dan | Spike Timing-Dependent Plasticity of Neural Circuits[END_REF][START_REF] Gerstner | A neuronal learning rule for sub-millisecond temporal coding[END_REF][START_REF] Gerstner | Spiking neuron models[END_REF][START_REF] Sjöström | Dendritic Excitability and Synaptic Plasticity[END_REF]]. In future studies both structural plasticity and STDP should also be taken into account for the GPe. Unlike STN, GPe is an inhibitory population [START_REF] Shink | The subthalamic nucleus and the external pallidum: Two tightly interconnected structures that control the output of the basal ganglia in the monkey[END_REF]. The structural plasticity mechanism used in this paper [START_REF] Butz | Homeostatic structural plasticity increases the efficiency of small-world networks[END_REF][START_REF] Butz | A Simple Rule for Dendritic Spine and Axonal Bouton Formation Can Account for Cortical Reorganization after Focal Retinal Lesions[END_REF] 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 [START_REF] Vogels | Inhibitory synaptic plasticity: spike timing-dependence and putative network function[END_REF].

Figure 1 .

 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 STN neurons, modified by STDP and structural plasticity. The external electrical stimulation delivered to the STN is indicated by the grey 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 PD patient prior to 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 equation (19) which represents the change in synaptic elements 𝑑𝑧/𝑑𝑡 as a function of the firing rate FR(𝑡). 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.

Figure 2 .

 2 Figure 2. Temporal evolution of 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 𝑇 of each cycle. Each stimulation site is activated exactly once during an ON cycle, where the sequence of activation across stimulation sites is random.

Figure 3 .

 3 Figure 3. Periodic Stimulation (PS)-Coordinated Reset (CR) sequences without structural plasticity. (A) Time evolution of the order parameter 𝑅(𝑡) 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 grey the CR intervals (10 min). STDP is active throughout the entire trial. The synchronizationdesynchronization 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 .

 4 Figure 4. Three identical Periodic Stimulation (PS)-Coordinated Reset (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 panel (A) and (C). (C) Time course of the total number of synaptic excitatory connections in the STN population.

Figure 5 .

 5 Figure 5. PS-CR 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.

  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 to 1000 times more slowly than the electrical activity in the neurons. The homeostatic rules applied to the

	𝑑FR 𝑑𝑡	=	{	--	FR(𝑡) 𝜏 SP 𝜏 SP FR(𝑡)	+ 𝛽, ,	if the neuron fires otherwise	,

synaptic elements are intended to keep the mean firing rate stable. The firing rate FR(𝑡) is calculated in NEST by each neuron following this equation:

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). This research was partially funded by the Helmholtz Association through the Helmholtz Portfolio Theme "Supercomputing and Modeling for the Human Brain". This project also received funding from the European Union's Horizon 2020 research and innovation program under grant agreement No 720270 (HBP SGA1) and No 785907 (HBP SGA2). The authors 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 also 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. TM is also supported by the Labex MME DII (ANR-11-LBX-0023-01) French national funding program. PAT 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).

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.

Conflict of Interest

PAT works as consultant for Boston Scientific Neuromodulation and Gretap AG and is inventor or coinventor on a number of patents for invasive and non-invasive neuromodulation. The other 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.

Author