Coordinated reset stimulation in a large-scale model of the STN-GPe circuit

Abstract

Synchronization of populations of neurons is a hallmark of several brain diseases. Coordinated reset (CR) stimulation is a model-based stimulation technique which specifically counteracts abnormal synchrony by desynchronization. Electrical CR stimulation, e.g., for the treatment of Parkinson's disease (PD), is administered via depth electrodes. In order to get a deeper understanding of this technique, we extended the top-down approach of previous studies and constructed a large-scale computational model of the respective brain areas. Furthermore, we took into account the spatial anatomical properties of the simulated brain structures and incorporated a detailed numerical representation of 2 · 10⁴ simulated neurons. We simulated the subthalamic nucleus (STN) and the globus pallidus externus (GPe). Connections within the STN were governed by spike-timing dependent plasticity (STDP). In this way, we modeled the physiological and pathological activity of the considered brain structures. In particular, we investigated how plasticity could be exploited and how the model could be shifted from strongly synchronized (pathological) activity to strongly desynchronized (healthy) activity of the neuronal populations via CR stimulation of the STN neurons. Furthermore, we investigated the impact of specific stimulation parameters especially the electrode position on the stimulation outcome. Our model provides a step forward toward a biophysically realistic model of the brain areas relevant to the emergence of pathological neuronal activity in PD. Furthermore, our model constitutes a test bench for the optimization of both stimulation parameters and novel electrode geometries for efficient CR stimulation.

1. Introduction

Several brain disorders are associated with abnormal neuronal synchronization. For example, epileptic seizures are characterized by an abnormally synchronized firing of neurons during the epileptic process (Lehnertz et al., 2009). Parkinson's disease (PD) is another pathological condition of the central nervous system, which is associated with abnormal neuronal synchronization in specific brain regions (Lenz et al., 1994; Nini et al., 1995; Plenz and Kital, 1999). The occurrence of resting tremor with frequencies between 4 and 9 Hz in PD is correlated to pathological synchronization of particular neuronal populations located in the thalamus and basal ganglia (Deuschl et al., 2000; Bevan et al., 2002; Mallet et al., 2008; Smirnov et al., 2008; Tass et al., 2010), which actually seem to drive the tremor (Smirnov et al., 2008; Tass et al., 2010). In contrast, in the healthy state these populations fire in an uncorrelated manner (Nini et al., 1995; Brown et al., 2001; Fahn, 2008).

An established therapy for patients with medically refractory PD is electrical deep brain stimulation (DBS) (Benabid et al., 1991; Rodriguez-Oroz et al., 2005; Deuschl et al., 2006). Today the field of application of DBS is expanding. For instance, DBS has recently been evaluated in clinical trials in patients suffering from treatment-resistant depression, and first-in-man studies have been conducted for the treatment of Tourette's syndrome (Anderson et al., 2012; Cannon et al., 2012). Despite positive clinical results, the understanding of the mechanisms of action underlying the therapeutic effects of DBS is still a matter of debate (Grill and McIntyre, 2001; Volkmann, 2004; Anderson et al., 2012).

For DBS depth electrodes are implanted in particular target areas during high precision stereotactic surgery (Voges et al., 2002; Coffey, 2009). For the treatment of PD, the STN is the most commonly used target for DBS (Rodriguez-Oroz et al., 2005; Deuschl et al., 2006). This small lens-shaped brain structure is part of the basal ganglia, which are important for voluntary movement control (Nini et al., 1995). The electrodes mostly used for implantation comprise four cylindrical contacts, which are made of a Pt-Ir alloy. Each contact can individually be activated for stimulation (Coffey, 2009). Extension cables connect the electrodes to the pacemaker, which is implanted below the patient's collarbone and serves as counter-pole during stimulation. The standard DBS algorithm is high-frequency (HF) permanent stimulation, where a biphasic current pulse train is continuously administered via the implanted electrode contact(s). Stimulation frequencies of 100 − 180 Hz and amplitudes varying between 0 and 10 V have clinically proven to be effective (Moro et al., 2002; Volkmann et al., 2002). The mechanisms of DBS are still matter of debate (McIntyre et al., 2004c). One hypothesis is, that HF DBS leads to a suppression of activity in the target structure and thus mimics lesioning effects in a reversible way (McIntyre et al., 2004c; Volkmann, 2004; Kühn et al., 2008; Montgomery and Gale, 2008; Pyragas et al., 2013). However, there is evidence for other non-local mechanisms (McIntyre et al., 2004c; Gradinaru et al., 2009). Rubin et al. (2012) gives a good review of the complex mechanisms of DBS based on experimental and numerical observations.

HF DBS treatment allows for a significant reduction of medication (Deuschl et al., 2006). However, there is a relevant non-responder rate, for instance, tremor suppression is non-effective in 15% of the patients, and in 73% of the patients HF DBS is non-effective with respect to axial symptoms (Limousin et al., 1999; Volkmann, 2004). In fact, not all symptoms are equally reduced (Rodriguez-Oroz et al., 2005), some patients show side effects such as worsening speech, reduced verbal fluency, and weight gain affect an important number of patients (Voges et al., 2002). Some of these side effects may result from current spread to remote structures or improper electrode placement, but the underlying reasons for the occurrence of side effects are still not sufficiently understood (Voges et al., 2002).

To overcome limitations of HF DBS, with methods from non-linear dynamics and statistical physics, stimulation techniques have been developed, which specifically counteract pathological neuronal synchronization by desynchronization (Tass, 1999). Coordinated reset (CR) stimulation is a robust desynchronization technique which can be applied in an open loop as well as closed loop mode (Tass, 2003b,c).

CR aims at a therapeutic reshaping of neuronal connectivity by harnessing synaptic plasticity (e.g., spike timing-dependent plasticity, STDP) (Gerstner et al., 1996; Markram et al., 1997). Based on Shen et al. (2003) we hypothesize that STDP also occurs within the STN population. The goal of CR stimulation is to decrease the rate of coincidences, which – due to STDP – causes a decrease of the synaptic weights and, in turn, an unlearning of pathological connectivity and synchrony (Tass and Majtanik, 2006). Put otherwise, CR DBS aims at long-lasting effects that persist after cessation of stimulation (Tass and Majtanik, 2006; Tass and Hauptmann, 2007). First experimental in vitro studies in rat hippocampal slice and in vivo experiments in Parkinsonian monkeys confirmed that CR causes a long-lasting desynchronization (Tass et al., 2009) and reduction of symptoms (Tass et al., 2012b). By the same token, long-lasting and cumulative aftereffects of CR stimulation were observed in a proof of concept study in PD patients (Adamchic et al., 2014a). Computational studies showed that CR stimulation can be delivered both invasively and non-invasively (Popovych and Tass, 2012; Tass and Popovych, 2012). Accordingly, CR has successfully been applied to the treatment of tinnitus via acoustic stimulation (Tass and Popovych, 2012; Tass et al., 2012a; Silchenko et al., 2013; Adamchic et al., 2014b).

The biophysical mechanism of action of DBS is not yet sufficiently understood (Grill and McIntyre, 2001; Volkmann, 2004). In order to investigate how electrical stimulation affects neuronal tissue, computational models of the respective systems are useful tools (Rubin and Terman, 2004; Miocinovic et al., 2006; So et al., 2012). In the present study we focus on the effects of different electrical stimulation algorithms on the collective activity of a model neuronal network. So far, the mechanisms of CR stimulation have been studied utilizing neuron models with reduced complexity and a limited number of modeled neurons (Hauptmann et al., 2005; Tass and Hauptmann, 2009; Guo and Rubin, 2011; Lysyansky et al., 2011). Therefore, in this study we consider a sufficient number of neurons required to appropriately sample the stimulated volume in space and establish a computational platform for time consuming simulations taking into account the slow STDP-mediated dynamics.

The high number and density of neurons gave us the opportunity to investigate local and propagation effects of synchronization within the network. Hence, we constructed a large-scale model of the two structures hypothesized to be responsible for the pathological activity associated with PD. Our model contains a network of in total 2 · 10⁴ neurons. In contrast to previous models, where the nodes of the simulated networks were placed on one- or two-dimensional (regular) lattices (Hauptmann et al., 2005; Guo and Rubin, 2011; Lysyansky et al., 2011), and similar to a previous study (Hauptmann and Tass, 2010) we arranged the neurons within a spatial network in order to replicate the three-dimensional structure of the simulated parts of the brain. We used a conductance based, and biophysically realistic model associated with physical dimensions for the individual neurons (Terman et al., 2002; Rubin and Terman, 2004).

With increased complexity, however, it becomes challenging to fully understand the dynamics of the system. Hence, we mainly based our model on experimentally constrained parameters and continued the top-down approach of previous studies (Hauptmann et al., 2005; Hauptmann and Tass, 2007; Maistrenko et al., 2007; Tass and Hauptmann, 2007, 2009; Hauptmann and Tass, 2009, 2010; Lysyansky et al., 2011; Popovych and Tass, 2012) by increasing the complexity of the considered model gradually. The model presented in this study can be considered as a step forward toward a biophysically realistic model of a major target area for DBS.

2. Methods

We used Neural Simulation Tool (NEST) (version 2.0.0-rc4) (Gewaltig and Diesmann, 2007) to implement and perform the simulations presented in this study. The simulation code of the model network was implemented using the Python interface for NEST (Eppler et al., 2008). NEST allows for an implementation of three dimensional models of neuronal networks via the topology module (Plesser and Austvoll, 2009). We performed all simulations on the high-performance computer JUROPA at the Research Center Jülich, Germany. This supercomputer is composed of 3288 compute nodes, while each node comprises two Intel Xeon X55704 quad-core processors, resulting in 26304 processors in total available for computation. The supercomputer is equipped with a main memory capacity of 79 TB and provides 274.8 FLOPS performance, measured by the LINPACK-benchmark. Accordingly, each compute node has access to 24 GB of memory. This high amount of memory available per node allows for simulations with complex individual neurons and an online analysis of simulation data. The various computation nodes are linked via an InfiniBand quad data rate fat tree topology. The used Message Passing Interface (MPI) implementation was Parastation-MPI 5.0.25. The differential equations describing the dynamics of our model network were solved numerically utilizing an embedded Runge-Kutta-Fehlberg method of fourth order with adaptive step size control, which is implemented as part of the GNU Scientific Library (GSL) (version 1.14) (Fehlberg, 1970; Galassi et al., 2009).

2.1. Model equations

In order to model the neuronal activity of the individual STN and GPe neurons, we used the model developed by Terman and Rubin (Terman et al., 2002; Rubin and Terman, 2004), which is a single-compartment conductance based model.

The Na⁺ and K⁺ ionic currents I_Na, I_K are responsible for the spike production. Besides those spike-producing currents a low-threshold T-type Ca²⁺ current I_T, a high-threshold Ca²⁺ current I_Ca along with a K⁺ current I_ahp, which kicks in after hyperpolarization, enter the equation. The intrinsic neuron dynamics is additionally influenced by synaptic inputs from connected STN and GPe neurons I_syn, the input from surrounding brain structures I_noise, and an external stimulation current I_stim. Via this stimulation current the influence of DBS on each single neuron is modeled. The various ionic currents are all expressed in pA/μm² and described by several equations:

The equations for both STN and GPe neurons are very similar. Besides different parameter values, the only difference is the slightly simpler form of the low threshold Ca²⁺ current I_T = g_Ta³_∞ (v) r [v − v_Ca] for the GPe neurons. In our model the parameters for the reversal potentials and maximum conductances of the respective ion channels were Gaussian distributed around the mean values given in Tables 1, 2. We chose a standard deviation equal to 5% of the mean value. Thus, not all neurons were identical, which led to a more realistic scenario. A schematic diagram of our model network is illustrated in Figure 1. We focussed on a detailed description of the connections within and between the two nuclei, the STN and the GPe. Hence, our model consists of two mutually interacting sub-networks representing the STN and the GPe. The connectivity matrix of the entire network could be written as a combination of several sub-networks: w_ss, w_gs, w_sg, and w_gg represent the connection matrices for the STN-STN (ss), STN-GPe (sg), GPe-STN (gs), and GPe-GPe (gg) sub-networks, respectively. Since we considered the simulated neuronal network as a weighted and directed graph, the matrix was not necessarily symmetric and each entry W_ij corresponded to the connection strength between neuron i and j. We excluded self-connections in our model, which resulted in: W_ii = 0∀_i.

Figure 1

2.2. Three-dimensional model of the STN and the GPe

In our network we modeled 10⁴ neurons each within the STN and the GPe. We had to make a trade-off between a physiologically realistic number of neurons in our system and the required computation time. Because the GPe fills a larger volume and has a lower density of neurons than the STN, it was reasonable to choose the same number of neurons for both nuclei (Levesque and Parent, 2005). In order to incorporate anatomical constraints conceding the spatial location of the STN and the GPe into our model, we made use of magnetic resonance imaging (MRI) data taken from a PD patient before DBS surgery. Both nuclei of the left brain hemisphere were segmented on slices of the MRI data by an experienced neurosurgeon. The segmentation contours are illustrated in green in Figure 2. In a second step we manually fitted an ellipsoid [given by (x/a)² + (y/b)² + (z/c)² − 1 ≤ 0] into each of the obtained segmentation contours in order to approximate the anatomical shape of both structures. These two ellipsoid structures were filled with uniformly randomly distributed pointlike neurons. In this way we created a spatial model of the STN and the GPe, which is illustrated in Figure 2. The approximated values of the semi-principal axes for the two ellipsoids were: a_stn = 2.5 mm, b_stn = 6.0 mm, c_stn = 3.0 mm for the STN, and a_gpe = 4.6 mm, b_gpe = 12.3 mm, c_gpe = 3.2 mm for the GPe. Therefore, the resulting STN volume was V_stn = 188.5 mm³, which is in good accordance to a mean volume of 174.5 ± 20.4 mm³ reported in the literature (Levesque and Parent, 2005). The GPe volume is not presented here, because it is not the stimulation target and, thus, in the framework of our simulations the geometric dimensions are not crucial.

Figure 2

Although it is still up for discussion whether synaptic connections between STN neurons exist, there has been substantial experimental and theoretical evidence, which supports our assumption of the existence of excitatory connections among STN neurons (Gillies and Willshaw, 1998, 2004; Gillies et al., 2002; Shen et al., 2003; Shen and Johnson, 2006; Baufreton and Bevan, 2008). We based our connectivity value on a theoretical estimation given in Gillies and Willshaw (1998). Hence, each STN neuron extended connections to p_ss = 7% of the entire STN population, resulting in 700 synapses per neuron and 7 · 10⁶ synapses in total. The synaptic delays for all connections within the STN were extrapolated from values for gg connections reported in Holgado et al. (2010) and set to a homogeneous value of δ_ss = 4.0 ms. According to Holgado et al. (2010) we assumed that the corresponding time delays within both structures are shorter than the delays between GPe and STN. The connection probability between STN neurons declined exponentially with increasing distance (see Figure 3) and was given by p(x) = e^−x/c_d. The basis for our proposed distance dependence was the extent of the dendritic tree surrounding a STN neuron. In accordance with (Hellwig, 2000) we assumed that a synaptic contact is present when an axonal and a dendritic branch are spatially close to each other. Because the density of dendritic branches is lower for larger distances from the soma, the probability of a touch (and thus a synapse) between a dendritic tree and an axon decreases with increasing distance (Hellwig, 2000). The range of the dendrites varies from 0.2 to 0.9 mm between individual STN neurons and has a mean of r_d = 0.543 mm (Kita et al., 1983; Afsharpour, 1985). We adjusted the decay constant c_d = 0.5 such that the connection probability for the mean range of the dendritic tree r_d was still 33% and the probability for connections above the theoretical maximum distance of 1.9 mm was considerably low. While we chose a distance dependent probability for connections between neurons, the initial synaptic weights of the connections were Gaussian distributed around a mean value w_ss(0) with standard deviation σ_ss without any distance dependence.

Figure 3

In contrast to the STN, the connectivity within the GPe was significantly lower with each neuron contacting only 1% of the entire population (Sadek et al., 2007), resulting in 100 synapses per neuron and 10⁶ synapses in total for our model. Connections between GPe neurons were inhibitory. We chose the same distance dependence as reported above, but adjusted the parameter c_d = 0.63 in accordance to the increased reach of the dendritic tree, which ranged to a maximum of 1.5 mm with a mean of r_d = 0.7 mm (Yelnik et al., 1984; Sadek et al., 2007). For the synaptic delay we chose δ_ss = 4.0 ms as suggested in Holgado et al. (2010). Similar to the connections within the STN, the synaptic weights were Gaussian distributed with a mean value of w_gg = 0.25 · 10⁻³ and standard deviation σ_gg = 0.125 · 10⁻³.

Besides the connections within the two simulated parts of the basal ganglia, both structures are also tightly interconnected (Shink et al., 1996). First, the STN has an excitatory influence on the GPe and on the other hand the GPe exerts an inhibitory impact on the STN neurons. These long range connections are established by axons extending from the STN to the GPe and vice versa. Therefore, we did not use distance dependent connection probabilities for links between STN and GPe, but we relied only on the connectivity values of 2% reported in Baufreton et al. (2009). This resulted in a fixed number of 200 synapses per neuron, contacting randomly chosen targets in the distant nucleus, respectively. The transmission delay of δ_sg = δ_gs = 4.0 ms for connections between the STN and the GPe is published in Fujimoto and Kita (1993) and Kita et al. (2005). Likewise to the above described connections, the initial values for the synaptic weight were Gaussian distributed around the mean values w_sg = 6.0 · 10⁻³ and w_gs = 3.0 · 10⁻³ with respective standard deviations σ_sg = 0.3 · 10⁻³ and σ_gs = 0.15 · 10⁻³.

2.3. Synaptic currents

In our model each incoming spike triggers a postsynaptic current with the shape of an α-function (Dayan and Abbott, 2000; Gerstner and Kistler, 2002). The postsynaptic current generated by a spike at time t_k reads:

This yields the following coupling term, which quantifies the total synaptic input current to a postsynaptic neuron i received from presynaptic neurons j:

W_ij denotes the synaptic weight (coupling strength) between presynaptic neurons j and the postsynaptic neuron i. v_syn is the reversal potential, at which a given neurotransmitter causes no net current flow of ions through the ion channels associated with that particular neurotransmitter receptor. 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. The exact values of v_syn and τ_syn for the four different types of connections within our model are given in Tables 1, 2.

2.3.1. Noise inputs

The external inhibitory input from the striatum to the GPe is modeled by a constant negative input current of I_app = −7.0 pA applied to all GPe neurons and additional Poissonian spike trains with a frequency of f^gpe_p = 40 Hz.

with the Poisson distribution P_λ(i) and λ = f_pδ with the observation (simulated) time interval of the Poisson distribution δ between subsequent time steps, which is equal to the time resolution of the simulation. Likewise, we mimick the excitatory connections from the cortex to the STN also by feeding Poissonian distributed spikes with a frequency of f^stn_p = 20 Hz into each STN neuron. We choose f^gpe_p = 2 · f^stn_p, so that the physiological mean firing rate in the GPe is approximately twice the mean firing rate of the STN, as reported from experiments (Steigerwald et al., 2008; Nishibayashi et al., 2011; Wichmann et al., 2011). We do not feed the same spike train to all neurons, but each neuron receives a unique Poissonian spike train. For both the STN and the GPe the connection strength for the noise input is w_noise = 0.2, the value of the time constants for the α-function is τ_noise = 1.0 ms, and the reversal potential is chosen to be equal to v_noise = 0 mV.

2.4. Synaptic plasticity

We consider the synaptic strength of connections between STN neurons to be dynamical variables. Therefore, we introduce a learning rule that included synaptic weights, which are dependent on the timing of the firing of pre- and postsynaptic STN neurons (Gerstner et al., 1996; Markram et al., 1997; Bi and Poo, 2001). Our hypothesis that STN connections are governed by STDP is based on the findings reported in Shen et al. (2003), where long-term depression (LTD) and potentiation (LTP) could be induced by stimulation in rat STN neurons.

We use an additive spike timing-dependent plasticity (STDP) rule, which is illustrated in Figure 4A and defined by the following equation:

Figure 4

In a previous study it has been shown that in a Kuramoto model of coupled phase oscillators multistability only occurs for asymmetric STDP rules (Maistrenko et al., 2007). This asymmetry was substantiated experimentally in different brain regions (e.g., hippocampus Bi and Poo, 1998 and neocortex Feldman, 2000). Therefore, we also used an asymmetric STDP rule with the following parameters: τ₋ = 27.5, τ₊ = 12.0, λ = 2.0 · 10⁻³, and β = 1.1. τ₋ and τ₊ are the time constants for the synaptic weight change of depression and potentiation, respectively. λ quantifies the impact of a single spike on the synaptic weight change and can be considered as the learning rate of the synaptic connection. The parameter β specifies the ratio between depression and potentiation in the synaptic learning rule. The synaptic weights among STN neurons are confined to the interval w^ss_ij(t) ∈ [0, w^ss_max], where w^ss_max = 20 · 10⁻³. This upper boundary prohibits an unphysiological non-saturated increase in synaptic weights. We use hard bounds for the upper and lower boundary in our simulations. We use the same parameters for all synaptic connections among STN neurons. In order to investigate the dynamics of the model system with plasticity involved, we varied the initial mean synaptic weight among STN neurons w_ss(0) and observed the time course of the mean synaptic weight w_ss(t). The initial values for the mean coupling were Gaussian distributed with a standard deviation of 5% of the maximum weight σ_ss = 0.1 · 10⁻³. These simulations show that the system approaches one of the two dominating attractors with strongly synchronized activity (high mean coupling strength) or strongly desynchronized activity (low mean coupling strength). Below a certain threshold value in the initial mean synaptic weight w_ss(0) = 5.0 · 10⁻³ the system monotonically approaches the desynchronized state, and with an initial value above the threshold the system monotonically converges to the strongly synchronized state (cf. Figure 4B). The plot shows only the time course of the mean synaptic weight among STN neurons within the first 350 s of neuronal activity. Nevertheless, we simulated the network activity for 600 s, in order to exclude transient or meta-stable behavior. Although STDP might induce many different multistable states, only the strongly desynchronized and strongly synchronized states predominate, and the system ultimately approaches either one of them. For instance, a strong initial connection strength among STN neurons leads to an inception of synchronized activity. This causes the neurons to fire nearly simultaneously, and this narrow spike timing results in a potentiation of the coupling strength and thus stabilizes the strong coupling and associated synchronization. In contrast, weak initial coupling is associated with desynchronized dynamics within the system and therefore results in a further decrease in synaptic connection strength. Hence, the desynchronized and synchronized dynamics are stable in the model system governed by plasticity. STDP facilitates the stabilization of synchronized and desynchronized activity and thus leads to multistability in the STN population. This observation is in good accordance with results from previous studies (Tass and Majtanik, 2006; Hauptmann and Tass, 2007; Popovych and Tass, 2012). The effects described above are robust against variations in the standard deviation of the initial mean coupling strength of up to 20% of the maximum weight.

2.4.1. Stimulation input

For DBS the electrical stimulation is administered to the neuronal tissue in form of short biphasic current or voltage pulses (Cogan, 2008). A typical current pulse P(t) for stimulation of neuronal tissue has 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): where t_k are 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. This charge-balance is important in order to avoid permanent charge transfer into the neuronal tissue, which might cause permanent tissue damage. For details of the possible electrochemical processes occurring at the electrode-tissue interface we refer to Cogan (2008).

The spatial arrangement of the stimulated neurons and the stimulation sites is crucial if the impact of electrical stimulation on a neuronal population is investigated. The influence of electrical stimulation on an individual neuron strongly depends on the relative positions within the three-dimensional coordinate system. One of the clinically most commonly used electrodes for DBS is the Medtronic DBS lead model 3389. The distal end with four electrode contacts is inserted into the target area. The electrode consists of four separate cylindrical contacts made of a Pt-Ir alloy with a typical length of 1.5 mm (Coffey, 2009). The electric field of such a cylindrical electrode is approximately radial and decreases inversely proportional to the distance. We used an analytic approximation based on the electric field of a line charge of finite length in order to estimate the distance dependence of the stimulation strength (Richardson et al., 2003). Furthermore, we assumed a homogeneous resistivity of 1.25 Ωm for the surrounding neuronal tissue. The approximation reads (Richardson et al., 2003):

The function S(d_ij) models the distance dependent decay of stimulation strength with d_ij representing the distance between neuron i and the location of the stimulation contact j. Additionally, S(d_ij) depends on the length of the electrode contacts l_c.

In order to eliminate the singularity at d_ij = 0, we assume a minimal distance of d_min = 0.7 mm to the electrode for all neurons. We achieve this by prohibiting possible neuron coordinates within a cylindrical volume with 1.4 mm diameter surrounding the electrode axis. This electrode canal is slightly larger than the electrode diameter and is created during the implantation surgery.

We apply CR stimulation (Tass, 2003a,c) to the modeled network of STN neurons. Each neuron receives stimulation input from all M stimulation contacts. The distance dependence of the stimulation strength causes the formation of M different sub-populations, because only the neurons close to a particular stimulation site are influenced sufficiently in order to alter their firing behavior. According to the CR stimulation protocol M different sub-populations are stimulated via M different stimulation sites at equidistant timing points by application of short bursts of current pulses within a stimulation period T, which approximately matches the period of the targeted pathological oscillation. T_p denotes the stimulation period within a burst. For our simulations we chose M = 4, which matches the number of available electrode contacts for the considered electrode type Medtronic 3389. In contrast to the sequential CR algorithm, we randomized the order of activated contacts between stimulation periods (i.e., stimulation cycles) (Figure 5). The only constraint was that the same contact could not be activated twice in a row between subsequent cycles. This modification of the CR stimulation protocol has previously been successfully used in an animal experiment with Parkinsonian monkeys (Tass et al., 2012b). Three cycles of CR stimulation (ON-cycles) were followed by two cycles without stimulation (OFF-cycles). According to a previous computational study this alternating pattern of 3 stimulation ON- and 2 stimulation OFF-cycles is superior to permanent CR DBS (without OFF-cycles), since during the OFF-cycles the CR-induced desynchronization is maximal (Lysyansky et al., 2011). Throughout our simulations we only varied the stimulation amplitude κ. The remaining stimulation parameters were fixed to T_p = 7.69 ms, ω = 200 μs, p_s = 8, and T = 125.0 ms.

Figure 5

Therefore, the effective stimulation signal received by neuron i from stimulation contact j is given by multiplication of the stimulation amplitude with the spatial profile defined by Equation 15 and reads:

The function ρ_j(t) is an indicator function, which controls the activation of the respective stimulation site j. M is the number of activated stimulation contacts, where M = 4 for a typical electrode.

In our model we inject the stimulation current exclusively to the soma of the stimulated cells. Thus, our model omits various proposed effects of DBS (e.g., axonal activation). However, a recent computational study showed, that DBS of local STN cells located within the nucleus is more beneficial than the stimulation of fibers of passage (So et al., 2012). Nevertheless, our investigations revealed that the observed mixture of firing rates within the STN population is consistent with the interplay of effects reported or predicted for HF stimulation (Grill and McIntyre, 2001; McIntyre et al., 2004a; Miocinovic et al., 2006; Hauptmann and Tass, 2007). While some STN neurons showed stimulus-locked firing, mimicking the stimulus-locked somatic or axonal generation of action potentials in STN neurons, other STN neurons were completely or partially inhibited. Similar to Humphries and Gurney (2012) this mixture entirely arose from the interaction between the stimulation pulses fed into STN neurons and the effects induced by the reciprocal connections within the simulated STN-GPe network. This effect has also been reported previously for other simplified models of electrical stimulation (Humphries and Gurney, 2012). However, it is unquestioned that the mechanism of action of DBS includes suppression of synaptic input and direct (including retrograde) activation of axons (Ranck, 1975; Grill and McIntyre, 2001; McIntyre et al., 2004a; Miocinovic et al., 2006; Hauptmann and Tass, 2007; Gradinaru et al., 2009).

2.5. Simulation observables

Previous studies (Kuramoto, 1984; Daido, 1992; Tass, 1999) have shown that the phenomenon of synchronization within a neuronal model system can be quantified by the order parameter R(t). The order parameter R(t) bases on the phases ϕ_j(t) of the individual neurons of the neuronal model system and is defined by: n denotes the number of neurons contained in the system. ψ(t) is the mean phase of the neuronal system and the quantity ϕ_j(t) represents the phase of each individual neuron j, which is defined by linear interpolation between two subsequent spikes at times t_j,k and t_j,k+1 (k = 0, 1, 2, …). Thus, the phase of a neuron ranges in the interval [0, 2π) and is given by the following equation (Pinsky and Rinzel, 1995):

The above definition of a neuronal phase requires an adequate spike detection. We detected spikes by a combined threshold and local maximum search: If there was a local maximum above a certain threshold of the membrane potential, it was considered a spike. We performed the spike detection with an accuracy of 0.1 ms, which was sufficient according to a typical spike width at half-maximal spike amplitude of approximately 2.0 ms. We did not distinguish between single spikes and bursts (episodes of rapid spiking). R(t) quantifies the coherence of the phases and is ranging from 0 to 1, where R(t) = 1 corresponds to a fully synchronized system, while values close to zero correspond to complete absence of phase synchronization. Due to finite size fluctuations R(t) never reaches zero, is the expected value. We measured the degree of synchronization separately for the STN and GPe sub-networks. For example, in the desynchronized regime the time course of the first order parameter R(t) is stable and close to zero in both neuronal populations with time averages of 〈R〉_stn = 0.009 and 〈R〉_gpe = 0.009, which are in good accordance with the expected value of R(t) = = 0.01. In contrast, in the synchronized regime R(t) is stable and close to 1 for the STN sub-network with an average value of 〈R₁〉_stn = 0.97 indicating almost perfect in-phase synchronization. For the GPe population the values of the order parameter are considerably lower with a time average of 〈R₁〉_gpe = 0.67. The weakly coupled GPe neurons are periodically synchronized by the external force of the STN. In order to smooth out short-term fluctuations and highlight long-term trends, we calculated a moving average over the recorded data. Thereby we used a time window of 1000 ms duration.

In order to predict an optimal location for the DBS electrode we additionally studied the spatial distribution of synchronization within the STN volume. Therefore, we introduced a method to measure local synchronization within the STN, based on the order parameter described above. In our approach we subdivided the volume of the STN in equally sized cubic voxels and computed the order parameter R_k(t) for the set of N_k neurons contained within each voxel k separately.

Synchrony was calculated for each voxel individually. The voxels do not overlap. The high density of neurons in our large-scale model guarantees that each voxel contains approximately N_k ≈ 42 neurons. Nevertheless, we only considered voxels containing at least 10 neurons. For voxels with a low number of contained neurons the distinction between partly and fully synchronized states can become difficult due to the finite size fluctuations mentioned above. This also limits the spatial resolution of this synchrony measure. In order to estimate the homogeneity of local synchronization within the STN volume, we define the following local order parameter: r₁ is close to one if the entire STN volume is strongly synchronized and close to zero if it is strongly desynchronized. N_V = 232 denotes the number of voxels considered for the calculation. Note, r₁ is close to one if either the whole STN population is synchronized and constitutes a giant oscillator, or if the STN population is synchronized on the level of single voxels only, whereas between different voxels there might be non-vanishing phase shifts. The standard deviation of this local order parameter σ_r1 incorporated additional information about the homogeneity of the local synchronization. Large values of σ_r1 indicate high discrepancies in the local order parameter between different regions of the STN volume.

3. Results

3.1. Scalability test

The large number of simulated neurons made it necessary to parallelize the numerical computations. In our performance tests, we focused on the strong scaling of our model. Strong scaling characterizes the ability of an application to perform the same sized calculation faster on computers with more processors (Plesser et al., 2007). The network was built according to the scheme described above. Thus, the resulting network contained 2 · 10⁴ nodes connected via 1.2 · 10⁷ synapses. Only the synaptic weights between the STN neurons were governed by a STDP rule, while all remaining connections were constant in time. We focused on the simulation of the network activity and omitted the online data analysis in the scalability tests. The time evolution of this network was simulated for one biological second and the runtime was measured. To test the strong scaling behavior of the simulation program, the problem size was fixed and the number of used processors on the supercomputer was increased from 8 to 2048 processors. During simulations the neurons of the entire simulated network were active with a mean firing rate of 15.8 Hz. The scalability of our simulation is illustrated in Figure 6. The speed-up factor quantifies how much faster the simulation runs on n processors, with runtime T_n, compared to a simulation on 8 processors, (T₈). Thus, the speed-up factor is given by S = T₈/T_n. For up to 128 used processors our scalability test showed a slightly steeper slope than the ideal linear speed-up (see Figure 6B, i.e., for up to 128 processors the speed-up marks “+” are positioned above the ideal line of linear speed-up). A maximal slope for the speed-up of 1.1 was observed (averaged over 10 independent trials). Although the linear speed-up is commonly accepted to be the theoretical maximum, this supralinear behavior of our simulation can be explained by the fast cache memory increasing proportionally to the number of used processors (Morrison et al., 2005). Thus, the ratio between the locally required working memory and the locally available cache memory decreases, which leads to a performance gain (Wilkinson and Allen, 2004). Note, the supralinear speed-up might be more prominent for simulated networks with higher average firing rates (Morrison et al., 2005). This is because the size of the available cache memory has a stronger impact on the efficiency of writing to the spike buffers for higher firing rates of the simulated neurons.

Figure 6

The computational workload of the individual neurons in a simulated neuronal network has a significant impact on the scaling behavior of the simulation (Plesser et al., 2007). This is due to the fact that a computationally expensive model for individual neurons or a larger number of neurons per processor leads to an increase in the ratio of local computation to communication costs. The overhead due to the communication in a distributed simulation grows with an increasing number of used processors (Plesser et al., 2007). In our performance test this effect balanced the performance gain from efficient caching for 256 processors. The communication overhead became predominant for 2048 processors, which led to a saturation of the gained speed-up. In the simulations we performed throughout this work, we typically used 1024 processors. With this setup the simulation of one biological second consumed approximately 3 s of computation time. This included the required computation time for the online analysis of the simulation results, which was omitted in the scalability tests.

3.1.1. Spontaneous activity of the model network

As shown above the system enters one of the two dominant stables states depending on the initial mean synaptic weight among the STN neurons. The observed spontaneous activity within the two states can be characterized in the following way: First, for low mean coupling within the STN the neurons fire in an uncorrelated, desynchronized manner. We consider this type of activity as a model for the non-pathological (healthy) state of the basal ganglia system. Second, the system is capable to fire in a strongly correlated synchronized pattern, if the mean coupling strength within the STN is strong enough. This type of activity models the pathological (“Parkinsonian”) state of the STN-GPe network.

Both states can be characterized by multiple observables, e.g., the mean membrane potentials of both populations that have been recorded with 1.0 ms time resolution. In the desynchronized state the mean membrane potentials of the STN and the GPe population showed only minor fluctuations around the resting value. The time averages over the mean membrane potentials were 〈v(t)〉_stn = −59.2 ± 0.1 mV and 〈v(t)〉_gpe = −90.2 ± 0.1 mV.

In order to measure the degree of synchronization within the simulated neuronal network, we therefore calculated the order parameter R(t) (section 2.5), which measured the phase coherence of the neurons and gave us the opportunity to detect possible cluster states. The time course of the order parameter R(t) was stable and close to zero in both neuronal populations with a time average of 〈R(t)〉_stn = 0.009 for the STN and 〈R(t)〉_gpe = 0.009 for the GPe with standard deviations smaller than 10⁻⁴. The time course of the order parameter for both sub-structures was in accordance to what could be expected from a system of noise driven and very weakly coupled oscillators.

An adequate characterization of the observed neuronal activity required the investigation of the firing frequencies of individual neurons and their distribution. Therefore, we calculated the time intervals between subsequent spikes for each neuron separately and thus obtained the distribution of interspike intervals over the two populations of neurons. Hence, this resulted in an interspike interval (ISI) histogram shown in Figure 7. The bin width of the histogram was chosen equal to 5.0 ms. In the desynchronized state the ISI histogram of the STN population showed a broadened distribution of ISIs around the median value of _stn = 261.2 ms with an average absolute deviation from the median (MD) of MD_stn = 51.1 ms. Accordingly, this value corresponded to a median firing rate of _stn = 3.8 Hz. For the GPe neurons we calculated a median spike period of _gpe = 92.1 ± 82.6 ms, which corresponded to a median spike frequency of _gpe = 10.9 Hz. The ISI distribution for the GPe was asymmetrical and is shown in Figure 7B. The asymmetry occurred due to the fact that the GPe neurons have a tendency to fire short bursts of action potentials. These bursts had relatively short ISIs resulting in a shift of the histogram to ISIs with reduced length. We preferred calculating the median over the mean value, because both distributions were asymmetrical and for a symmetrical distribution mean and median have identical values.

Figure 7

For strong mean coupling among STN neurons (w_ss > 5.0 · 10⁻³) the model system showed highly synchronized activity (R(t) ≈ 1). The synchronized state can be considered as the pathological, unhealthy state, which mimics the neuronal activity within the basal ganglia system in patients suffering from PD (Lenz et al., 1994; Plenz and Kital, 1999). The STN neurons fired collectively with a median period of _stn = 122.0 ± 8.3 ms, which corresponded to a median firing rate of = 8.2 Hz. Between these episodes of collective firing the individual STN neurons remained almost silent. The GPe neurons were driven by the STN activity and showed collective firing with the same frequency as the STN neurons. Nevertheless, individual GPe neurons emitted action potentials constantly and collective periods of silence were not observed for the GPe population.

The computation of population mean membrane potentials gives similar results. The mean membrane potentials of both populations showed high amplitude oscillations with the same frequency as the synchronized activity of the STN neurons (cf. Figure 8). The synchronized firing of the STN neurons forced periodic collective spike discharges of the GPe neurons. This collective GPe activity led to an increased inhibition of the connected STN cells. Thus, the mean membrane potential of the STN population was strongly hyperpolarized (cf. Figure 8A). This mechanism was the reason for the periods of silence in STN neurons following the episodes of collective firing. However, the inhibition from GPe neurons decayed and the mean membrane potential of the STN population recovered before the next episode of collective firing was launched. The GPe neurons, which were excited by incoming action potentials from collective firing of the STN neurons, reacted to this excitation with the emission of bursts. This bursting activity in the GPe neurons caused the very narrow peaks in the mean membrane potential of the GPe population (cf. Figure 8B).

Figure 8

A very clear indication of synchronized activity in the simulated network was furthermore provided by the time course of the order parameter R(t). For the STN sub-network R(t) was stable and close to 1 with an average value of 〈R(t)〉_stn = 0.97 indicating almost perfect in-phase synchronization. However, perfect in-phase synchronization was never achieved, because the STN sub-network was composed of non-identical individual neurons and the system always received noise inputs from surrounding structures. For the GPe population the values of the order parameter were considerably lower oscillating around a time average of 〈R(t)〉_gpe = 0.67. This finding can be explained by the fact that the weakly coupled GPe neurons fired in a desynchronized manner most of the time and are periodically synchronized by the external force of the STN. The synchronized input from the STN triggered the collective bursting of the GPe neurons and afterwards the GPe population returned to the desynchronized firing pattern.

The effect of synchronization on the distribution of ISIs is clearly visible in Figure 9. In the synchronized state the median period of spiking activity in STN neurons was shifted toward shorter durations _stn = 122.0 ± 8.3 ms. The consequence was an increased frequency _stn = 8.2 Hz and much sharper distribution of ISIs (cf. Figure 9A), than for the desynchronized network activity, shown in Figure 7A. The calculation of the median or mean did not make sense for the bi-modal distribution of ISIs in the GPe population shown in Figure 9B. The bi-modal distribution was caused by two reasons. First, the large peak for short time periods around 5 − 10 ms could be explained by the increased bursting activity of the GPe neurons. Second, the clearly visible peak at approximately 122.0 ms was caused by the driving force of connected STN neurons. Additionally, the entire ISI histogram for the GPe population indicated an increased activity distributed over all time periods. This was due to the fact that the frequency of the input from the STN to the GPe was increased in the synchronized regime.

Figure 9

3.2. Optimal stimulation amplitude

In all subsequent simulations we set the initial mean synaptic weight among STN neurons to w_ss = (18.0 ± 0.1) · 10⁻³. This corresponded to an initially strongly synchronized activity of the model system. In the first 120 s of the simulations we left the system unperturbed before we switched on the CR stimulation. As described above, we applied the randomized CR stimulation sequences (Figure 5). The contacts of the stimulation electrode were placed along the y-axis of the coordinate system. The M = 4 electrode contacts were centered within the STN volume with a distance of 2 mm between two neighboring stimulation sites. In our simulations we only varied the amplitude of the stimulation current pulses κ and the location of the electrode contacts. The remaining stimulation parameters were fixed and given by: ω = 200.0 μs, p_s = 8.0, T_p = 7.69 ms, M = 4, T = 125.0 ms, and number of OFF- and ON-cycles= 3 and 2. The first question we addressed was to look for an “optimal” stimulation amplitude. An optimal stimulation amplitude results in sufficiently strong effects, i.e., desynchronization of the target structure, while using a minimum stimulation amplitude to avoid potential side effects. Therefore, we increased the applied stimulation amplitude in a stepwise manner and investigated the desychronizing effect of the respective stimulation current amplitude κ. We performed this investigation for fixed electrode contact locations in the center of the STN volume. The value of the local order parameter r₁ obtained after stimulation application is used as a quantitative measure for the desynchronization capabilities of a particular stimulation amplitude. The order parameter was averaged in time (interval t ∈ [1140, 1200] s, i.e., the last 60 s of simulated neuronal activity). The averaging started 420 s after cessation of stimulation. Figure 10 illustrates the dependence of the local order parameter r₁ on the stimulation amplitude κ. The illustration shows that for low stimulation amplitudes a significant decrease in synchronized activity could not be achieved. For instance, with an amplitude of κ = −1.0 mA the system remained in a strongly synchronized state after stimulation, and the low standard deviation for r₁ indicated that the entire STN volume was homogeneously synchronized. This was because for these low amplitudes the stimulation signal was not strong enough, and neurons in further distance from the electrode contacts did not sense a sufficiently strong input current in order to achieve a reliable phase reset. Hence, the amount of neurons recruited by stimulation was not large enough for very low stimulation amplitudes.

Figure 10

With further increasing amplitude, the fraction of neurons sensing significant stimulation inputs linearly increased and the electrical CR stimulation unfolded its desynchronizing effect. Nevertheless, large standard deviations of r₁ for stimulation amplitudes below −3.2 mA indicated the existence of large discrepancies in the spatial distribution of the local order parameters. For illustration, the effects of CR stimulation with amplitude κ = −2.0 mA on the model system are shown in Figure 11. The mean synaptic weight was successfully decreased during stimulation, but not sufficiently enough in order to shift the system into an overall strongly desynchronized state. After the stimulation the mean synaptic weight increased again and stabilized at w_ss(t → ∞) ≈ 10.4 · 10⁻³. Accordingly, the order parameters were reduced, but both values after stimulation (〈R₁〉_stn = 0.7, 〈R₁〉_gpe = 0.63) indicated the persistence of reduced, but still highly correlated activity within the STN and GPe sub-networks. Figure 11C shows the spatial distribution of the local order parameter after CR stimulation and provides an explanation for these observations. Only the inner part of the STN volume was reliably desynchronized, while synaptic connections and thus synchronization among neurons at larger distances from the electrode could not be reduced due to an insufficient stimulation amplitude. However, these patterns of local (de-)synchronization after CR application were stable within our observation times. Above we showed that for spatially homogeneous initial conditions the system approaches one of the two dominating attractors with strongly synchronized activity (high mean coupling strength) or strongly desynchronized activity (low mean coupling strength) cf. Figure 4B. Our observations show that the coupling strength between neurons can be influenced locally by electrical CR stimulation and thus the system can be forced to approach one of the stable attractors shown in Figure 4B if the fraction of manipulated connections is large enough.

Figure 11

This figure also indicates that the synaptic connection strengths between STN neurons, which were governed by a STDP learning rule, converged either to values close to the maximum synaptic weight w^max_ss = 20.0 · 10⁻³ or to values close to the minimum of w^min_ss = 0. The minimal stimulation amplitude, which was capable to achieve a persistent desynchronization of the complete STN volume, was κ_opt = −3.3 mA. The stimulation outcome did not differ qualitatively for stimulation amplitudes above κ_opt. However, the effective desynchronization of the STN volume broke down for amplitudes above −6.0 mA. This effect could be explained by a reduced spatial specificity of the stimulation for larger amplitudes, i.e., the regions stimulated by the different electrode contacts showed high intersections among each other, which led to stimulation of the entire STN population instead of small sub-populations of STN neurons.

For intermediate amplitudes (i.e., κ ∈ [−3.3 mA, −6.0 mA]) and during active CR stimulation the STN is divided into four synchronized sub-populations. The four synchronized sub-populations origin from the four stimulation sites, i.e., the periodic stimulation at each individual stimulation site causes a local synchronization of the surrounding neurons. Through the four stimulation sites a phase-shifted stimulation pattern (i.e., the typical CR pattern) is applied, which results in an overall weakening of synaptic weights and thus decreases the synchronized activity of the whole STN population already during stimulation.

For larger amplitudes almost all STN neurons are synchronized and form one population which is firing in phase. This synchronized firing does not decrease synaptic weights.

In Figure 12 we illustrated the effects of CR stimulation at appropriate intensity, by considering the results obtained with the minimal stimulation amplitude enabling complete desynchronization, κ_opt = −3.3 mA. Before the stimulation was switched on the mean synaptic weight among STN neurons was stable and close to its maximum value. During the application of CR stimulation the mean weight w_ss monotonically decayed. After 10 min the continuous CR stimulation was switched off and without further stimulation the mean synaptic weight stabilized close to zero (Figure 9A). The order parameters for the STN and the GPe population decreased during stimulation, but saturated and fluctuated around a value of approximately R₁ = 0.35. This is due to the fact that during application of CR stimulation sub-populations of the STN are periodically synchronized before they undergo a desynchronization transient. Nevertheless, after switching stimulation off both neuronal populations fired in a strongly desynchronized manner and remained in this state.

Figure 12

3.3. Influence of electrode location

In order to investigate the influence of the stimulation electrode position within the STN volume on the desynchronizing effect of CR stimulation, we varied the electrode position, i.e., the location of the array of stimulation contacts and analyzed the stimulation outcome in comparison to an ideal reference position. The topology of the array of stimulation contacts remained constant. The electrode position with the four electrode contacts centered in the STN volume served as reference position. Starting from this reference position, we consecutively moved the electrode along the three respective coordinate axes out of the center of the target in steps of 1.0 mm. For all simulations we used the stimulation amplitude of κ = −3.3 mA and the same stimulation parameters as before. The mean local order parameter r₁ was computed over the last 60 s of simulated neuronal activity t ∈ [1140, 1200] s and after cessation of stimulation. In Figure 13 we display the dependence of the local order parameter r₁ after CR stimulation on the displacement D along the three respective coordinate axes from the reference position. With the electrode located at the reference position the desynchronization of the STN volume was most effective. With further distance from the reference position the local order parameter r₁ increased monotonically for growing displacements of the electrode irrespective of the considered coordinate axis.

Figure 13

The dependence of r₁ on the displacement along the x- and z-axis was very similar. This was caused by the fact that the ellipsoid defining the boundary of the STN volume was almost symmetrical in x- and z-direction with values for the semi-principal axes of a_stn = 2.5 mm and c_stn = 3.0 mm. The slope of the increase in r₁ was significantly lower for the displacement along the y-axis, because the volume of neurons stimulated with a sufficiently high amplitude decreased with a lower rate than for displacement of the electrode along the two other coordinate axes. Therefore, we can conclude that the effect of electrode displacement on the stimulation outcome had a more significant impact for deviations along the x- and z-axes than along the y-axis.

However, as mentioned before, the increase in stimulation amplitude is limited by several factors. First, an effective desynchronization could not be achieved for amplitudes which exceeded a certain (model dependent) threshold, as showed above. Second, larger stimulation amplitudes might induce unwanted side effects via unintentional stimulation of other brain areas or fiber tracts (Smith et al., 2012). Third, high stimulation current amplitudes might induce permanent damage to the neuronal tissue surrounding the electrode (Grill and McIntyre, 2001; Grill, 2005; Cogan, 2008).

4. Discussion

In this study we developed a large-scale model of two basal ganglia structures STN and GPe, which are hypothesized to generate rhythmic synchronized oscillations. This synchronized activity is associated with symptoms arising in PD (Lenz et al., 1994; Plenz and Kital, 1999). The described model is the continuation of our top-down approach (Tass, 1999, 2003c; Popovych et al., 2005; Tass and Majtanik, 2006; Hauptmann and Tass, 2007, 2009, 2010; Tass and Hauptmann, 2009). By gradually increasing the complexity of our model, we were able to link our findings to previous results obtained from models with a reduced level of detail. Therefore, we increased the model size, parallelized the simulations and data analysis. Thus we were able to achieve a simulation performance close to real-time. We investigated and characterized the two predominant stable synchronized and desynchronized states, which modeled the pathological and healthy state of the basal ganglia system, respectively. We showed that the transition between both predominant states was dependent on the mean coupling strength among STN neurons. This observation is in good accordance with previous modeling studies (Hauptmann and Tass, 2007, 2009; Tass and Hauptmann, 2007, 2009). In addition, with our large-scale approach we were able to reveal clinically relevant spatio-temporal effects: spatially incomplete desynchronization caused by either insufficient stimulation or by electrode displacement.

Furthermore, we investigated the effects of CR stimulation on the network activity and, in particular, on the synaptic structure of the STN sub-network. In this context we were especially interested in how the neuronal network can be shifted from the strongly synchronized to the desynchronized mode of network activity. To this end, we introduced a synaptic learning rule and a measure enabling to assess the distribution of synchronized neuronal activity within the volume of the STN. In our investigations of the desynchronizing CR stimulation technique, we observed a certain amplitude range for which desynchronization of the target structure was most effective. This finding might be useful for further experimental studies. This optimal parameter range has been reported previously for simpler models of coupled phase oscillators, FitzHugh-Nagumo spiking neurons, FitzHugh-Rinzel bursting neurons, and Hodgkin-Huxley neurons with and without STDP (Guo and Rubin, 2011; Lysyansky et al., 2011; Popovych and Tass, 2012). Additionally, in an experimental study in parkinsonian non-human primates it was shown that CR DBS at low intensity has long-lasting after effects on basal ganglia motor function (Tass et al., 2012b) and a clinical proof-of-concept trial revealed that CR DBS has lasting aftereffects in patients with Parkinsons disease as assessed by UPDRS III scores and local field potential theta and beta activity (Adamchic et al., 2014a).

In this study we expanded these observations in a more complex and microanatomically more realistic model. We incorporated more realistic coupling topologies and took the spatial arrangement of the electrode and stimulated neurons into account. However, the observed amplitude parameters highly depend on how the electric field surrounding the electrode is modeled (Richardson et al., 2003; Butson et al., 2006; Chaturvedi et al., 2010; Schmidt and van Riemen, 2012). Therefore, a generalized quantitative prediction of optimal stimulation parameters is difficult to obtain from a theoretical simulation study. Additionally, the spatial structure enabled us to vary the electrode location within the stimulation target and thus we observed a strong relation between the electrode position and the stimulation outcome. This observation substantiates previous theoretical and experimental studies, which investigated the impact of the electrode position on the stimulation outcome for HF DBS (Voges et al., 2002; Paek et al., 2011; Guo et al., 2012; Nakano et al., 2012; Reese et al., 2012).

Because of the high computational costs of our complex model, we limited our observations of CR DBS effects to timescales of several 10 of minutes. Therefore, we used a STDP learning rule which showed adequate effects on such timescales. Our model could adequately reproduce the experimentally observed long-term desynchronization effects of CR stimulation, which last on timescales from minutes to hours and days, demonstrated in slice experiments (Tass et al., 2009), a pre-clinical proof-of-concept study (Tass et al., 2012b) and a proof-of-concept trial (Adamchic et al., 2014a). Certainly, the learning rule has a significant impact on the simulation results. Therefore, more experimental data specifying the actual plasticity mechanisms in the basal ganglia and especially within the STN are required. In order to be able to simulate the long-term effects (timescales of days or weeks) of CR DBS, a modification of the STDP rule has to be considered. However, these long-term simulations necessitate a further improvement of the simulation code in terms of scalability and performance or a reduction of the model complexity. At this point less complex models might be more suitable for the investigation of long-term DBS effects, but with steadily increasing performance and computing capacity these long-term simulations will be possible in the medium-term.

We considered many physiologically relevant parameters that have been neglected in former simplified models of the STN-GPe network (e.g., spatial extent, realistic neuronal connectivity, and dimensionful neuron model) (Hauptmann et al., 2005; Feng et al., 2007; Tass and Hauptmann, 2009; Guo and Rubin, 2011; Lysyansky et al., 2011). This complex modeling approach had the advantage that characteristics observed in the model activity might be more transferable to an experimental verification. However, there are still some limitations and simplifications in our model. We used a simplified approximation for the distance dependence of the electric field (and thus the stimulation strength) surrounding the electrode. Besides the assumption that the electric field strength behaves similar to the radial component of the field for a finite line charge, this approximation was based on a homogeneous and isotropic tissue conductance. In a future refinement of the presented model the implementation of more realistic representations of the electric field (e.g., based on finite element methods McIntyre et al., 2004b) is intended (Butson and McIntyre, 2005, 2006; Butson et al., 2006; Chaturvedi et al., 2010; Schmidt and van Riemen, 2012). However, if a “realistic" representation is desired it is also necessary to consider the movements of charge carriers in the surrounding tissue (e.g., extracellular fluid). Obviously, substantiation of such a modeling approach requires more experimental data. Also, in this study we trade-off model complexity, to reduce computational costs. Accordingly, we introduced STDP only within the STN, since this is the primary stimulation target for CR DBS.

By the same token, future studies might incorporate more complex STDP rules including soft bounds (van Rossum et al., 2000). The effects of hard vs. soft bounds in a similar system were studied in detail in Popovych et al. (2013). That study provides a systematical investigation of the effects of hard and soft bounds on spontaneous dynamical states (in particular, multistability) as well as stimulation-driven dynamical states (stochastic resonance induced by uncorrelated noise) of a network. The stochastic resonance was robust with respect to these modifications. Furthermore, multistability, an essential ingredient for CR to induce long-lasting effects, as well as robust anti-kindling was observed for hard bounds (Tass and Majtanik, 2006) as well as for soft bounds (Maistrenko et al., 2007) in phase oscillator networks. Correspondingly, to demonstrate that the CR-induced anti-kindling does not depend on the particular choice of bounds in the model used in this paper, among other future developments we are planning to take into account soft bounds in a forthcoming study, too. Moreover, further studies might be used to investigate more complex phenomena in the field of STDP, for example the observation that in the presence of propagation delays in recurrent networks exhibiting population bursts, STPD rules might exert a strong decoupling force which can further the desynchronization of the network activity (Lubenov and Siapas, 2008). Also, STDP can induce oscillatory behaviors depending also on the model parameters and in particular on the shape of the considered pulse (Mikkelsen et al., 2014). Further studies might also incorporate neighboring structures which are omitted in the current approach. For example, GPe fibers and axonal pathways, which might be relevant for a detailed analysis of the effects of variations of the electrode position (Rubin et al., 2012). Additional measures, focussing on the effect of information transfer on antidromic structures might further increase the understanding of the therapeutic effect of different sites and types of stimulation (Rubin and Terman, 2004; Rubin et al., 2012). Note, optimal electrode placement for standard HF DBS might differ quite substantially from that of CR DBS. Computational results show that CR is effective no matter whether neurons are stimulated at somata or via inhibitory or excitatory synapses (Popovych and Tass, 2012). Accordingly, in our paper we concentrated on the central electrode placement, since from a clinical standpoint it might be feasible already with available depth electrodes.

Our investigations of different electrode positions showed that the placement of the stimulation electrode within the target structure is crucial for an effective desynchronization. Furthermore, the required stimulation current amplitude was dependent on the electrode location. In the context of a clinically more effective electrode placement it is still under discussion, whether the optimal stimulation site for HF DBS lies within the interior of the STN or in fiber tracts located close to the STN (Voges et al., 2002; Paek et al., 2011; Guo et al., 2012; Nakano et al., 2012; Reese et al., 2012). In order to tackle this issue utilizing a modeling approach, a future version of our model should include spatially extended multi-compartment neurons instead of pointlike model neurons. Hence, the mentioned fiber tracts should be represented by multi-compartment cable models.

Another feature of DBS relevant to possible further improvement is the shape of the stimulation pulses used. In various studies the efficacy of particular pulse shapes for neuronal stimulation was investigated (Kuncel and Grill, 2004; Kuncel et al., 2006; Butson and McIntyre, 2007; Foutz and McIntyre, 2010; Wongsarnpigoon and Grill, 2010; Hofmann et al., 2011). Apart from the stimulation pulse shapes, the other CR parameters (e.g., stimulation frequency, intra-burst frequency etc.) deserve further detailed investigations. Accordingly, our work provides a tool to probe related stimulation effects, which might be observed experimentally in future pre-clinical and clinical studies. Our model might contribute to a clinically useful computational tool for the optimization of parameters for CR DBS. In particular, the spatially extented topology of our model might be important for the development of new stimulation electrode geometries (Buhlmann et al., 2011; Martens et al., 2011).

Conflict of interest statement

Peter Tass is employed by Jülich Research Center and works as Consulting Professor at Stanford University; formerly working with ANM GmbH (Cologne, Germany), shareholder of ANM GmbH. Christian Hauptmann is employed by Jülich Research Center; formerly working with ANM GmbH (Cologne, Germany). Several patents protect electrical CR neuromodulation. The inventor of this patent is Peter Tass, and the assignee is Jülich Research Center. Martin Ebert reports no potential conflicts of interest.

References

• 1

  AdamchicI.HauptmannC.BarnikolU.PawelczykN.PopovychO.BarnikolT.et al. (2014a). Coordinated reset neuromodulation for Parkinson's disease: proof-of-concept study. Mov. Disord. 29, 1679–1684. 10.1002/mds.25923

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 2

  AdamchicI.TothT.HauptmannC.TassP. A. (2014b). Reversing pathologically increased EEG power by acoustic CR neuromodulation. Hum. Brain Mapp. 35, 2099–2118. 10.1002/hbm.22314

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 3

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

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 4

  AndersonR.FryeM.AbulseoudO.LeeK.McGillivrayJ.BerkM.et al. (2012). Deep brain stimulation for treatment-resistant depression: efficacy, safety and mechanisms of action. Neurosci. Biobehav. Rev. 36, 1920–1933. 10.1016/j.neubiorev.2012.06.001

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 5

  BaufretonJ.BevanM. D. (2008). D2-like dopamine receptor-mediated modulation of activity-dependent plasticity at GABAergic synapses in the subthalamic nucleus. J. Physiol. 586, 2121–2142. 10.1113/jphysiol.2008.151118

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 6

  BaufretonJ.KirkhamE.AthertonJ. F.MenardA.MagillP. J.BolamJ. P.et al. (2009). Sparse but selective and potent synaptic transmission from the globus pallidus to the subthalamic nucleus. J. Neurophysiol. 102, 532–545. 10.1152/jn.00305.2009

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 7

  BenabidA. L.PollakP.GervasonC.HoffmannD.GaoD. M.HommelM.et al. (1991). Long-term suppression of tremor by chronic stimulation of the ventral intermediate thalamic nucleus. Lancet337, 403–406. 10.1016/0140-6736(91)91175-T

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 8

  BevanM. D.MagillP. J.TermanD.BolamJ. P.WilsonC. J. (2002). Move to the rhythm: oscillations in the subthalamic nucleus–external globus pallidus network. Trends Neurosci. 25, 525–531. 10.1016/S0166-2236(02)02235-X

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 9

  BiG.-Q.PooM.-M. (1998). Synaptic modifications in cultured hippocampal neurons: Dependence on spike timing, synaptic strength, and postsynaptic cell type. J. Neurosci. 18, 10464–10472.

  • Pubmed Abstract
  • Google Scholar

• 10

  BiG.-Q.PooM.-M. (2001). Synaptic modification by correlated activity: Hebb's postulate revisited. Annu. Rev. Neurosci. 24, 139–166. 10.1146/annurev.neuro.24.1.139

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 11

  BrownP.OlivieroA.MazzoneP.InsolaA.TonaliP.LazzaroV. D. (2001). Dopamine dependency of oscillations between subthalamic nucleus and pallidum in Parkinson's disease. J. Neurosci. 21, 1033–1038.

  • Pubmed Abstract
  • Google Scholar

• 12

  BuhlmannJ.HofmannL.TassP. A.HauptmannC. (2011). Modeling of a segmented electrode for desynchronizing deep brain stimulation. Front. Neuroeng. 4:15. 10.3389/fneng.2011.00015

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 13

  ButsonC. R.MaksC. B.McIntyreC. C. (2006). Sources and effects of electrode impedance during deep brain stimulation. Clin. Neurophysiol. 117, 447–454. 10.1016/j.clinph.2005.10.007

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 14

  ButsonC. R.McIntyreC. C. (2005). Tissue and electrode capacitance reduce neural activation volumes during deep brain stimulation. Clin. Neurophysiol. 16, 2490–2500. 10.1016/j.clinph.2005.06.023

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 15

  ButsonC. R.McIntyreC. C. (2006). Role of electrode design on the volume of tissue activated during deep brain stimulation. J. Neural Eng. 3, 1–8. 10.1088/1741-2560/3/1/001

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 16

  ButsonC. R.McIntyreC. C. (2007). Differences among implanted pulse generator waveforms cause variations in the neural response to deep brain stimulation. Clin. Neurophysiol. 118, 1889–1894. 10.1016/j.clinph.2007.05.061

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 17

  CannonE.SilburnP.CoyneT.O'MaleyK.CrawfordJ. D.SachdevP. S. (2012). Deep brain stimulation of anteromedial globus pallidus interna for severe Tourette's syndrome. Am. J. Psychiatry169, 860–866. 10.1176/appi.ajp.2012.11101583

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 18

  ChaturvediA.ButsonC. R.LempkaS. F.CooperS. E.McIntyreC. C. (2010). Patient-specific models of deep brain stimulation: influence of field model complexity on neural activation predictions. Brain Stimul. 3, 65–77. 10.1016/j.brs.2010.01.003

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 19

  CoffeyR. J. (2009). Deep brain stimulation devices: a brief technical history and review. Arti. Organs33, 208–220. 10.1111/j.1525-1594.2008.00620.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 20

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

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 21

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

  • CrossRef
  • Google Scholar

• 22

  DayanP.AbbottL. (2000). Theoretical Neuroscience. Cambridge, MA: The MIT Press.

  • Google Scholar

• 23

  DeuschlG.RaethjenJ.BaronR.LindemannM.WilmsH.KrackP. (2000). The pathophysiology of Parkinsonian tremor: A review. J. Neurol. 247, 33–48. 10.1007/PL00007781

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 24

  DeuschlG.Schade-BrittingerC.KrackP.VolkmannJ.SchäferH.BötzelK.et al. (2006). A randomized trial of deep-brain stimulation for Parkinson's disease. N. Engl. J. Med. 355, 896–908. 10.1056/NEJMoa060281

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 25

  EpplerJ. M.HeliasM.MullerE.DiesmannM.GewaltigM.-O. (2008). PyNEST: a convenient interface to the NEST simulator. Front. Neuroinform. 2:12. 10.3389/neuro.11.012.2008

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 26

  FahnS. (2008). The history of dopamine and levodopa in the treatment of Parkinson's disease. Move. Disord. 23, 497–508. 10.1002/mds.22028

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 27

  FehlbergE. (1970). Klassische runge-kutta-formeln vierter und niedrigerer ordnung mit schrittweiten-kontrolle und ihre anwendung auf Wärmeleitungsprobleme. Computing6, 61–71. 10.1007/BF02241732

  • CrossRef
  • Google Scholar

• 28

  FeldmanD. (2000). Timing based LTP and LTD at vertical inputs to layer ii/iii pyramidal cells in rat barrel cortex. Neuron27, 45–56. 10.1016/S0896-6273(00)00008-8

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 29

  FengX.-J.Shea-BrownE.GreenwaldB.KosutR.RabitzH. (2007). Optimal deep brain stimulation of the subthalamic nucleus—a computational study. J. Comput. Neurosci. 23, 265–282. 10.1007/s10827-007-0031-0

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 30

  FoutzT. J.McIntyreC. C. (2010). Evaluation of novel stimulus waveforms for deep brain stimulation. J. Neural Eng. 7:066008. 10.1088/1741-2560/7/6/066008

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 31

  FujimotoK.KitaH. (1993). Response characteristics of subthalamic neurons to the stimulation of the sensorimotor cortex in the rat. Brain Res. 23, 185–192. 10.1016/0006-8993(93)90872-K

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 32

  GalassiM.DaviesJ.TheilerJ.GoughB.JungmanG.AlkenP.et al. (2009). GNU Scientific Library Reference Manual - Third Edition. Godalming: Network Theory Ltd.

  • Google Scholar

• 33

  GerstnerW.KempterR.van HemmenJ. L.WagnerH. (1996). A neuronal learning rule for sub-millisecond temporal coding. Nature383, 76–78. 10.1038/383076a0

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 34

  GerstnerW.KistlerW. M. (2002). Spiking Neuron Models. Cambridge: Cambridge University Press. 10.1017/CBO9780511815706

  • CrossRef
  • Google Scholar

• 35

  GewaltigM.-O.DiesmannM. (2007). NEST (NEural Simulation Tool). Scholarpedia2:1430. 10.4249/scholarpedia.1430

  • CrossRef
  • Google Scholar

• 36

  GilliesA. J.WillshawD. J. (1998). A massively connected subthalamic nucleus leads to the generation of widespread pulses. Proc. R. Soc. Lond. B265, 2101–2109. 10.1098/rspb.1998.0546

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 37

  GilliesA. J.WillshawD. J. (2004). Models of the subthalamic nucleus: the importance of intranuclear connectivity. Med. Eng. Phys. 26, 723–732. 10.1016/j.medengphy.2004.06.003

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 38

  GilliesA. J.WillshawD. J.LiZ. (2002). Subthalamic–pallidal interactions are critical in determining normal and abnormal functioning of the basal ganglia. Proc. R. Soc. Lond. B269, 545–551. 10.1098/rspb.2001.1817

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 39

  GradinaruV.MogriM.ThompsonK.HendersonJ.DeisserothK. (2009). Optical deconstruction of Parkinsonian neural circuitry. Science324, 354–359. 10.1126/science.1167093

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 40

  GrillW. M. (2005). Safety considerations for deep brain stimulation: review and analysis. Exp. Rev. Med. Dev. 2, 409–420. 10.1586/17434440.2.4.409

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 41

  GrillW. M.McIntyreC. C. (2001). Extracellular excitation of central neurons: implications for the mechanisms of deep brain stimulation. Thalamus Relat. Syst. 1, 269–277. 10.1017/S1472928801000255

  • CrossRef
  • Google Scholar

• 42

  GuoS.ZhuangP.HallettM.ZhengZ.ZhangY.LiJ.et al. (2012). Subthalamic deep brain stimulation for Parkinson's disease: correlation between locations of oscillatory activity and optimal site of stimulation. Parkinsonism Rel. Disord. 19, 109–114. 10.1016/j.parkreldis.2012.08.005

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 43

  GuoY.RubinJ. E. (2011). Multi-site stimulation of subthalamic nucleus diminishes thalamocortical relay errors in a biophysical network model. Neural Netw. 24, 602–616. 10.1016/j.neunet.2011.03.010

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 44

  HauptmannC.PopovychO.TassP. A. (2005). Effectively desynchronizing deep brain stimulation based on a coordinated delayed feedback stimulation via several sites: a computational study. Biol. Cybernet. 93, 463–470. 10.1007/s00422-005-0020-1

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 45

  HauptmannC.TassP. A. (2007). Therapeutic rewiring by means of desynchronizing brain stimulation. Biosystems89, 173–181. 10.1016/j.biosystems.2006.04.015

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 46

  HauptmannC.TassP. A. (2009). Cumulative and after-effects of short and weak coordinated reset stimulation: a modeling study. J. Neural Eng. 6, 1–13. 10.1088/1741-2560/6/1/016004

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 47

  HauptmannC.TassP. A. (2010). Restoration of segregated neuronal activity by desynchronizing stimulation in the presence of sensory input. J. Neural Eng. 7:056008. 10.1088/1741-2560/7/5/056008

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 48

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

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 49

  HofmannL.EbertM.TassP. A.HauptmannC. (2011). Modified pulse shapes for effective neural stimulation. Front. Neuroeng. 4:9. 10.3389/fneng.2011.00009

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 50

  HolgadoA. J. N.TerryJ. R.BogaczR. (2010). Conditions for the generation of beta oscillations in the subthalamic nucleus-globus pallidus network. J. Neurosci. 30, 12340–12352. 10.1523/JNEUROSCI.0817-10.2010

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 51

  HumphriesM. D.GurneyK. N. (2012). Network effects of subthalamic deep brain stimulation drive a unique mixture of responses in basal ganglia output. Eur. J. Neurosci. 36, 2240–2251. 10.1111/j.1460-9568.2012.08085.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 52

  KitaH.ChangH.KitaiS. T. (1983). The morphology of intracellularly labeled rat subthalamic neurons: a light microscopic analysis. J. Comp. Neurol. 215, 245–257. 10.1002/cne.902150302

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 53

  KitaH.TachibanaY.NambuA.ChikenS. (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. 10.1523/JNEUROSCI.1719-05.2005

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 54

  KühnA. A.KempfF.BrückeC.DoyleL. G.Martinez-TorresI.PogosyanA.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. 10.1523/JNEUROSCI.0282-08.2008

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 55

  KuncelA. M.CooperS. E.WolgamuthB. R.ClydeM. A.SnyderS. A.MontgomeryE. B.Jr.et al. (2006). Clinical response to varying the stimulus parameters in deep brain stimulation for essential tremor. Move. Disord. 21, 1920–1928. 10.1002/mds.21087

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 56

  KuncelA. M.GrillW. M. (2004). Selection of stimulus parameters for deep brain stimulation. Clin. Neurophysiol. 115, 2431–2441. 10.1016/j.clinph.2004.05.031

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 57

  KuramotoY. (1984). Chemical Oscillations, Waves, and Turbulence. Berlin: Springer. 10.1007/978-3-642-69689-3

  • CrossRef
  • Google Scholar

• 58

  LehnertzK.BialonskiS.HorstmannM.-T.KrugD.RothkegelA.StaniekM.et al. (2009). Synchronization phenomena in human epileptic brain networks. J. Neurosci. Methods183, 42–48. 10.1016/j.jneumeth.2009.05.015

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 59

  LenzF.KwanH.MartinR.TaskerR.DostrovskyJ.LenzY. (1994). Single unit analysis of the human ventral thalamic nuclear group - Tremor-related activity in functionally identified cells. Brain117, 531–543. 10.1093/brain/117.3.531

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 60

  LevesqueJ.-C.ParentA. (2005). GABAergic interneurons in human subthalamic nucleus. Move. Disord. 20, 574–584. 10.1002/mds.20374

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 61

  LimousinP.SpeelmanJ. D.GielenF.JanssensM. (1999). Multicentre european study of thalamic stimulation in Parkinsonian and essential tremor. J. Neurol. Neurosurg. Psychiatry66, 289–296. 10.1136/jnnp.66.3.289

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 62

  LubenovE. V.SiapasA. G. (2008). Decoupling through synchrony in neuronal circuits with propagation delays. Neuron58, 118–131. 10.1016/j.neuron.2008.01.036

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 63

  LysyanskyB.PopovychO. V.TassP. A. (2011). Desynchronizing anti-resonance effect of m:n ON–OFF coordinated reset stimulation. J. Neural Eng. 8:036019. 10.1088/1741-2560/8/3/036019

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 64

  MaistrenkoY. L.LysyanskyB.HauptmannC.BurylkoO.TassP. A. (2007). Multistability in the Kuramoto model with synaptic plasticity. Phys. Rev. E75:066207. 10.1103/PhysRevE.75.066207

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 65

  MalletN.PogosyanA.MártonL. F.BolamJ. P.BrownP.MagillP. J. (2008). Parkinsonian beta oscillations in the external globus pallidus and their relationship with subthalamic nucleus activity. J. Neurosci. 28, 14245–14258. 10.1523/JNEUROSCI.4199-08.2008

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 66

  MarkramH.LübkeJ.FrotscherM.SakmannB. (1997). Regulation of synaptic efficacy by coincidence of postsynaptic APs and EPSPs. Science275, 213–215. 10.1126/science.275.5297.213

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 67

  MartensH. C. F.ToaderE.DecréM. M. J.AndersonD. J.VetterR.KipkeD. R.et al. (2011). Spatial steering of deep brain stimulation volumes using a novel lead design. Clin. Neurophysiol. 122, 558–566. 10.1016/j.clinph.2010.07.026

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 68

  McIntyreC. C.GrillW. M.ShermanD. L.ThakorN. V. (2004a). Cellular effects of deep brain stimulation: model-based analysis of activation and inhibition. J. Neurophysiol. 91, 1457–1469. 10.1152/jn.00989.2003

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 69

  McIntyreC. C.MoriS.ShermanD. L.ThakorN. V.VitekJ. L. (2004b). Electric field and stimulating influence generated by deep brain stimulation of the subthalamic nucleus. Clin. Neurophysiol. 115, 589–595. 10.1016/j.clinph.2003.10.033

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 70

  McIntyreC. C.SavastaM.Kerkerian-Le GoffL.VitekJ. L. (2004c). Uncovering the mechanism(s) of action of deep brain stimulation: activation, inhibition, or both. Clin. Neurophysiol. 115, 1239–1248. 10.1016/j.clinph.2003.12.024

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 71

  MikkelsenK.ImparatoA.TorciniA. (2014). Sisyphus effect in pulse-coupled excitatory neural networks with spike-timing-dependent plasticity. Phys. Rev. E89:062701. 10.1103/PhysRevE.89.062701

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 72

  MiocinovicS.ParentM.ButsonC. R.HahnP. J.RussoG. S.VitekJ. L.et al. (2006). Computational analysis of subthalamic nucleus and lenticular fasciculus activation during therapeutic deep brain stimulation. J. Neurophysiol. 96, 1569–1580. 10.1152/jn.00305.2006

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 73

  MontgomeryE. B.Jr.GaleJ. T. (2008). Mechanisms of action of deep brain stimulation (DBS). Neurosci. Biobehav. Rev. 32, 388–407. 10.1016/j.neubiorev.2007.06.003

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 74

  MoroE. E.EsselinkR. J.XieJ.HommelM.BenabidA. L.PollakP. (2002). The impact on Parkinson's disease of electrical parameter settings in STN stimulation. Neurology59, 706–713. 10.1212/WNL.59.5.706

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 75

  MorrisonA.MehringC.GeiselT.AertsenA.DiesmannM. (2005). Advancing the boundaries of high-connectivity network simulation with distributed computing. Neural Comput. 17, 1776–1801. 10.1162/0899766054026648

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 76

  NakanoN.TanedaM.WatanabeA.KatoA. (2012). Computed three-dimensional atlas of subthalamic nucleus and its adjacent structures for deep brain stimulation in Parkinson's disease. ISRN Neurol. 2012, 1–13. 10.5402/2012/592678

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 77

  NiniA.FeingoldA.SlovinH.BergmanH. (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.

  • Pubmed Abstract
  • Google Scholar

• 78

  NishibayashiH.OguraM.KakishitaK.TanakaS.TachibanaY.NambuA.et al. (2011). Cortically evoked responses of human pallidal neurons recorded during stereotactic neurosurgery. Move. Disord. 26, 469–476. 10.1002/mds.23502

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 79

  PaekS. H.LeeJ.-Y.KimH.-J.KangD.LimY. H.KimM. R.et al. (2011). Electrode position and the clinical outcome after bilateral subthalamic nucleus stimulation. J. Korean Med. Sci. 26, 1344–1355. 10.3346/jkms.2011.26.10.1344

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 80

  PinskyP. F.RinzelJ. (1995). Synchrony measures for biological neural networks. Biol. Cybernet. 73, 129–137. 10.1007/BF00204051

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 81

  PlenzD.KitalS. T. (1999). A basal ganglia pacemaker formed by the subthalamic nucleus and external globus pallidus. Nature400, 677–682. 10.1038/23281

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 82

  PlesserH.EpplerJ.MorrisonA.DiesmannM.GewaltigM.-O. (2007). Efficient parallel simulation of large-scale neuronal networks on clusters of multiprocessor computers, in Euro-Par 2007 Parallel Processing, Volume 4641 of Lecture Notes in Computer Science, eds KermarrecA.-M.BougéL.PriolT. (Berlin: Springer), 672–681.

  • Google Scholar

• 83

  PlesserH. E.AustvollK. (2009). Specification and generation of structured neuronal network models with the NEST topology module. BMC Neurosci. 10, 1–2. 10.1186/1471-2202-10-S1-P56

  • CrossRef
  • Google Scholar

• 84

  PopovychO. V.HauptmannC.TassP. A. (2005). Effective desynchronization by nonlinear delayed feedback. Phys. Rev. Lett. 94:164102. 10.1103/PhysRevLett.94.164102

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 85

  PopovychO. V.TassP. A. (2012). Desynchronizing electrical and sensory coordinated reset neuromodulation. Front. Hum. Neurosci. 6:58. 10.3389/fnhum.2012.00058

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 86

  PopovychO. V.YanchukS.TassP. A. (2013). Self-organized noise resistance of oscillatory neural networks with spike timing-dependent plasticity. Sci. Rep. 11:2926. 10.1038/srep02926

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 87

  PyragasK.NovicenkoV.TassP. (2013). Mechanism of suppression of sustained neuronal spiking under high-frequency stimulation. Biol. Cybernet. 107, 669–684. 10.1007/s00422-013-0567-1

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 88

  RanckJ. B.Jr. (1975). Which elements are excited in electrical stimulation of mammalian central nervous system: a review. Brain Res. 98, 417–440. 10.1016/0006-8993(75)90364-9

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 89

  ReeseR.PinskerM. O.HerzogJ.WodargF.SteigerwaldF.Pötter-NergerM.et al. (2012). The atypical subthalamic nucleus–an anatomical variant relevant for stereotactic targeting. Move. Disord. 27, 544–548. 10.1002/mds.24902

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 90

  RichardsonK. A.GluckmanB. J.WeinsteinS. L.GloschC. E.MoonJ. B.GwinnR. P.et al. (2003). In vivo modulation of hippocampal epileptiform activity with radial electric fields. Epilepsia44, 768–777. 10.1046/j.1528-1157.2003.35402.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 91

  Rodriguez-OrozM. C.ObesoJ. A.LangA. E.HouetoJ.-L.PollakP.RehncronaS.et al. (2005). Bilateral deep brain stimulation in Parkinson's disease: a multicentre study with 4 years follow-up. Brain128(Pt 10), 2240–2249. 10.1093/brain/awh571

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 92

  RubinJ. E.McIntyreC. C.TurnerR. S.WichmannT. (2012). Basal ganglia activity patterns in parkinsonism and computational modeling of their downstream effects. Eur. J. Neurosci. 36, 2213–2228. 10.1111/j.1460-9568.2012.08108.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 93

  RubinJ. E.TermanD. (2004). High frequency stimulation of the subthalamic nucleus eliminates pathological thalamic rhythmicity in a computational model. J. Comput. Neurosci. 16, 211–235. 10.1023/B:JCNS.0000025686.47117.67

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 94

  SadekA. R.MagillP. J.BolamJ. P. (2007). A single-cell analysis of intrinsic connectivity in the rat globus pallidus. J. Neurosci. 27, 6352–6362. 10.1523/JNEUROSCI.0953-07.2007

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 95

  SchmidtC.van RiemenU. (2012). Modeling the field distribution in deep brain stimulation: the influence of anisotropy of brain tissue. IEEE Trans. Biomed. Eng. 59, 1583–1592. 10.1109/TBME.2012.2189885

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 96

  ShenK.-Z.JohnsonS. W. (2006). Subthalamic stimulation evokes complex EPSCs in the rat substantia nigra pars reticulata in vitro. J. Physiol. 573, 697–709. 10.1113/jphysiol.2006.110031

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 97

  ShenK.-Z.ZhuZ.-T.MunhallA.JohnsonS. W. (2003). Synaptic plasticity in rat subthalamic nucleus induced by high-frequency stimulation. Synapse50, 314–319. 10.1002/syn.10274

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 98

  ShinkE.BevanM.BolamJ.SmithY. (1996). The subthalamic nucleus and the external pallidum: two tightly interconnected structures that control the output of the basal ganglia in the monkey. Neuroscience73, 335–357. 10.1016/0306-4522(96)00022-X

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 99

  SilchenkoA.AdamchicI.HauptmannC.TassP. A. (2013). Impact of acoustic coordinated reset neuromodulation on effective connectivity in a neural network of phantom sound. Neuroimage77, 133–147. 10.1016/j.neuroimage.2013.03.013

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 100

  SmirnovD. A.BarnikolU. B.BarnikolT. T.BezruchkoB. P.HauptmannC.BührleC.et al. (2008). The generation of Parkinsonian tremor as revealed by directional coupling analysis. Eur. Lett. 83:20003. 10.1209/0295-5075/83/20003

  • CrossRef
  • Google Scholar

• 101

  SmithY.WichmannT.FactorS. A.DeLongM. R. (2012). Parkinson's disease therapeutics: new developments and challenges since the introduction of levodopa. Neuropsychopharmacology37, 213–246. 10.1038/npp.2011.212

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 102

  SoR. Q.KentA. R.GrillW. M. (2012). Relative contributions of local cell and passing fiber activation and silencing to changes in thalamic fidelity during deep brain stimulation and lesioning: a computational modeling study. J. Comput. Neurosci. 32, 499–519. 10.1007/s10827-011-0366-4

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 103

  SteigerwaldF.PötterM.HerzogJ.PinskerM.KopperF.MehdornH.et al. (2008). Neuronal activity of the human subthalamic nucleus in the Parkinsonian and non-Parkinsonian state. J. Neurophysiol. 100, 2515–2524. 10.1152/jn.90574.2008

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 104

  TassP. A. (1999). Phase Resetting in Medicine and Biology - Stochastic Modelling and Data Analysis. Berlin: Springer. 10.1007/978-3-540-38161-7

  • CrossRef
  • Google Scholar

• 105

  TassP. A. (2003a). Desynchronization by means of a coordinated reset of neural sub-populations. Prog. Theor. Phys. Suppl. 150, 281–296. 10.1143/PTPS.150.281

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 106

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

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 107

  TassP. A. (2003c). A model of desynchronizing deep brain stimulation with a demand-controlled coordinated reset of neural subpopulations. Biol. Cybernet. 89, 81–88. 10.1007/s00422-003-0425-7

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 108

  TassP. A.AdamchicI.FreundH.-J.von StackelbergT.HauptmannC. (2012a). Counteracting tinnitus by acoustic coordinated reset neuromodulation. Restor. Neurol. Neurosci. 30, 137–159. 10.3233/RNN-2012-110218

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 109

  TassP. A.QinL.HauptmannC.DoveroS.BezardE.BoraudT.et al. (2012b). Coordinated reset has sustained after-effects in Parkinsonian monkeys. Ann. Neurol. 72, 816–820. 10.1002/ana.23663

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 110

  TassP. A.HauptmannC. (2007). Therapeutic modulation of synaptic connectivity with desynchronizing brain stimulation. Int. J. Psychophysiol. 64, 53–61. 10.1016/j.ijpsycho.2006.07.013

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 111

  TassP. A.HauptmannC. (2009). Anti-kindling achieved by stimulation targeting slow synaptic dynamics. Restor. Neurol. Neurosci. 27, 589–609. 10.3233/RNN-2009-0484

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 112

  TassP. A.MajtanikM. (2006). Long-term anti-kindling effects of desynchronizing brain stimulation: a theoretical study. Biol. Cybernet. 94, 58–66. 10.1007/s00422-005-0028-6

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 113

  TassP. A.PopovychO. (2012). Unlearning tinnitus-related cerebral synchrony with acoustic coordinated reset stimulation: theoretical concept and modelling. Biol. Cybernet. 106, 27–36. 10.1007/s00422-012-0479-5

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 114

  TassP. A.SilchenkoA. N.HauptmannC.BarnikolU. B.SpeckmannE.-J. (2009). Long-lasting desynchronization in rat hippocampal slice induced by coordinated reset stimulation. Phys. Rev. E80:011902. 10.1103/PhysRevE.80.011902

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 115

  TassP. A.SmirnovD. A.KaravaevA.BarnikolU. B.BarnikolT. T.AdamchicI.et al. (2010). The causal relationship between subcortical local field potential oscillations and Parkinsonian resting tremor. J. Neural Eng. 7:016009. 10.1088/1741-2560/7/1/016009

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 116

  TermanD.RubinJ.YewA.WilsonC. (2002). Activity patterns in a model for the subthalamopallidal network of the basal ganglia. J. Neurosci. 22, 2963–2976.

  • Pubmed Abstract
  • Google Scholar

• 117

  van RossumM. V. W.BiG. Q.TurrigianoG. G. (2000). Stable hebbian learning from spike timing-dependent plasticity. J. Neurosci. 20, 8812–8821.

  • Pubmed Abstract
  • Google Scholar

• 118

  VogesJ.VolkmannJ.AllertN.LehrkeR.KoulousakisA.FreundH.-J.et al. (2002). Bilateral high-frequency stimulation in the subthalamic nucleus for the treatment of Parkinson disease: correlation of therapeutic effect with anatomical electrode position. J. Neurosurg. 96, 269–279. 10.3171/jns.2002.96.2.0269

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 119

  VolkmannJ. (2004). Deep brain stimulation for the treatment of Parkinson's disease. J. Clin. Neurophysiol. 21, 6–17. 10.1097/00004691-200401000-00003

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 120

  VolkmannJ.HerzogJ.KopperF.DeuschlG. (2002). Introduction to the programming of deep brain stimulators. Move. Disord. 17, 181–187. 10.1002/mds.10162

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 121

  WichmannT.DeLongM. R.GuridiJ.ObesoJ. A. (2011). Milestones in research on the pathophysiology of Parkinson's disease. Move. Disord. 26, 1032–1041. 10.1002/mds.23695

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 122

  WilkinsonB.AllenM. (2004). Parallel Programming: Techniques and Applications using Networked Workstations and Parallel Computers. Upper Saddle River, NJ: Prentice Hall.

  • Google Scholar

• 123

  WongsarnpigoonA.GrillW. M. (2010). Energy-efficient waveform shapes for neural stimulation revealed with a genetic algorithm. J. Neural Eng. 7:046009. 10.1088/1741-2560/7/4/046009

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 124

  YelnikJ.PercheronG.FrancoisC. (1984). A golgi analysis of the primate globus pallidus. ii. Quantitative morphology and spatial orientation of dendritic arborizations. J. Comp. Neurol. 227, 200–213. 10.1002/cne.902270206

  • Pubmed Abstract
  • CrossRef
  • Google Scholar