Abstract
A large-scale computational model of the basal ganglia network and thalamus is proposed to describe movement disorders and treatment effects of deep brain stimulation (DBS). The model of this complex network considers three areas of the basal ganglia region: the subthalamic nucleus (STN) as target area of DBS, the globus pallidus, both pars externa and pars interna (GPe-GPi), and the thalamus. Parkinsonian conditions are simulated by assuming reduced dopaminergic input and corresponding pronounced inhibitory or disinhibited projections to GPe and GPi. Macroscopic quantities are derived which correlate closely to thalamic responses and hence motor programme fidelity. It can be demonstrated that depending on different levels of striatal projections to the GPe and GPi, the dynamics of these macroscopic quantities (synchronisation index, mean synaptic activity and response efficacy) switch from normal to Parkinsonian conditions. Simulating DBS of the STN affects the dynamics of the entire network, increasing the thalamic activity to levels close to normal, while differing from both normal and Parkinsonian dynamics. Using the mentioned macroscopic quantities, the model proposes optimal DBS frequency ranges above 130 Hz.
Similar content being viewed by others
Avoid common mistakes on your manuscript.
1 Introduction
1.1 Basal ganglia connectivity
Parkinson’s disease (PD) and dystonia, including different types, belong to the most common movement disorders and hence pose a considerable health burden (Chesselet and Delfs 1996; Defazio 2010; de Lau and Breteler 2006). It is generally accepted that they arise from a dysfunction of the basal ganglia (BG), shown as simplified circuitry in Fig. 1 and result in hypokinetic or hyperkinetic symptoms, depending on which part of the circuitry is affected. In this circuitry, cortical glutamatergic projections are thought to activate GABAergic medium spiny neurons (MSN) and interneurons of corpus striatum (striatum in Fig. 1). From the medium spiny, inhibitory neurons, the so-called direct pathway, project to globus pallidus pars interna (GPi). The inhibition of GPi leads to activation of thalamus (THA) (via disinhibition); one can thus speculate that the thalamus faithfully responds to initial cortical signals (to initiate movement). In the so-called indirect pathway, the activation of striatum inhibits the globus pallidus pars externa (GPe) which projects to the subthalamic nucleus (STN), enhancing its activity. An increased STN activation, in turn, will lead to an increment of GPi activity, resulting in an inhibition of the thalamus (and hence reduction of locomotive activity Calabresi et al. 2014), see also Fig. 1. Parkinson’s disease constitutes a paradigmatic hypokinetic syndrome, which results from the degeneration of the substantia nigra (SN in Fig. 1). It is characterised by rigidity, tremor and hypokinesia, i.e. the inability to start movements fluently (DeLong 1990; Maiti et al. 2017). Looking at the circuitry, the motor symptoms can well be explained by a loss of both activating and inhibitory projections of the SN to the striatum, leading to disinhibition in the so-called indirect (striatum–GPe–STN–GPi–THA) and over-activation in the so-called direct pathways (striatum–GPi–THA) (Calabresi et al. 2014). Hyperkinetic syndromes, characterised by involuntary movements or muscle contractions, in turn, are generally thought to originate from functional or structural damage or degeneration of striatum (dystonia or choreatic syndromes) and of STN; ballistic syndromes. Again, the simplified circuitry depicted in Fig. 1 can serve to explain the functional outcome of e.g. striatal over-activation speculated to result in a shift of balance toward the so-called direct pathway in dystonias (Wichmann and Dostrovsky 2011).
1.2 Role of activity patterns
While the simplified circuitry suggests an explanation for the emergence of hyper- or hypokinetic syndromes, in fact it does not take into consideration the patterning of activity although information in the nervous system is actually conveyed by spatio-temporal activity patterns. Pattern propagation, however, will typically depend on nonlinear couplings between the interacting compartments of the system (Bevan et al. 2002). Hence, one cannot assume that inhibitory or excitatory activity will straight-forwardly propagate across the network. Importantly, both in PD and in dystonia, changes in oscillatory activity patterns in the basal ganglia and the cortex seem to be markers of the diseases (Eusebio and Brown 2007). Thus, in PD, synchronised beta-band activity in both cortex and STN seems to be associated with hypokinesia (Crowell et al. 2012; Kühn et al. 2008). This prominent beta-band is speculated to be caused by the synchronised network activity between the different basal ganglia nuclei (Schmidt et al. 2013). Indeed, a shift of network behaviour from autonomous oscillations of STN and GPe neurons, to synchronous abnormal low-frequency bursting, is observed (Bevan et al. 2002). These findings provide further support for the view that the basal ganglia use both the pattern and the rate of neuronal activity to encode information. Synchronous abnormal patterns of activity of this local circuit should also be reflected in similar changes regarding to the activity of the basal ganglia output and the thalamic activity.
1.3 Deep brain stimulation
Deep brain stimulation (DBS) has been deemed to be the most important innovation in movement disorder therapy and has revolutionised treatment for PD, dystonia and essential tremor patients, first having been approved for the latter condition by the FDA (Food and Drug Administration, USA) in 1997 (Krack et al. 2019). Moreover, DBS is currently being introduced also for the therapy of mental disorders such as depression and obsessive-compulsive disorder (Holtzheimer and Mayberg 2011). DBS improves levodopa-related motor complications in PD and often, motor symptoms in dystonia (Deuschl et al. 2006; Vidailhet et al. 2013). For PD, we know that the efficiency of the treatment depends strongly on the frequency of the stimulus; high-frequency stimulation (HFS), \(f>90\) Hz (McConnell et al. 2012) improves motor symptoms, while low frequencies are ineffective or worsen the motor disorders (Koeglsperger et al. 2019). Unfortunately, HFS may worsen frontal functions such as verbal fluency. In contrast to motor functions, very low frequency (10Hz) on STN DBS has significantly better results in verbal fluency (Wojtecki et al. 2006). Until now, the main mechanism of DBS remains elusive (Krack et al. 2019; Ashkan et al. 2017; Udupa and Chen 2015). The main hypotheses put forward so far are (Rubin and Terman 2004; Guo et al. 2008) (a) an inactivation of the target nuclei, (b) changes in transmitter release, (c) neuroprotective and electrostatic effects (i.e. structural plasticity), and network effects resulting in firing pattern alterations. In addition, DBS does not cause a total silencing of neuronal activity of the target nuclei, but rather complex responses. These include sequences of prolonged activation or activity reduction (Luo and Kiss 2016), possibly because of a dissociation between somatal and axonal activation (Holsheimer et al. 2000). One of the most attractive explanations at the moment is the disruption of hypersynchronised oscillations (Kühn et al. 2008). Indeed, recordings from animal models suggest that DBS in the STN results in more periodic and regular firing at higher frequencies in the thalamus (Xu et al. 2008). Thalamic neurons receive inhibitory signals from GPi, an output nucleus of the basal ganglia which provides inhibitory input to the thalamus. Thus GPi activity affects its thalamic targets, generating a possible pathway for STN-DBS, to modify basal ganglia–thalamocortical activity. As a consequence, STN-DBS alters Parkinsonian GPi activity which might improve thalamo-cortical fidelity (Guo et al. 2008; So et al. 2012; Santaniello et al. 2015).
1.4 Predicting network dynamics under DBS using computational approaches
From a computational perspective, one of the main obstacles to explain or predict DBS effects on network activity is the lack of a coherent framework which could bridge the different scales (Deco et al. 2008; Siettos and Starke 2016) of network models ranging from microscopic (cellular activity) to macroscopic (symptom) (Pavlides et al. 2012, 2015). An intermediate level, the mesoscopic level, is related to the dynamics of specific networks of neurons in different nuclei of the basal ganglia. These mesoscopic networks constitute the bridge between micro- and macroscales. Examples of these network dynamics and variations of activity pattern are confirmed in a number of animal model studies in Parkinson models. For instance, DBS in the STN leaves firing rates unaltered in GPe and increases the firing rate in the thalamus (i.e. the nucleus the GPi projects on). While single unit recordings thus indicate changes in firing rates in projection areas of the STN during DBS, one overarching motive of all observations is that DBS effectively changes firing patterns (McConnell et al. 2012; Xu et al. 2008; So et al. 2012, 2017; Dorval et al. 2010) as the essential element of its therapeutic success.
More specifically, in McConnell et al. (2012), DBS leads to reduced low-frequency neuronal oscillations in GPe and SNr, increased neuronal oscillations at the stimulation frequency, and increased phase locking with the stimulus pulses. Moreover, coherence within and across the GPe and SNr during HFS was reduced in the band of pathological low-frequency oscillations and increased in the stimulation frequency band. These findings provide evidence that effective high-frequency DBS suppresses low-frequency network oscillations and entrains neurons in the basal ganglia. Therefore, these results support the hypothesis that the effectiveness of HFS stems from its ability to override pathological firing patterns in the basal ganglia by inducing a new regularized pattern of synchronous neuronal activity.
In this spirit we propose a computational large-scale biophysical model related to the Parkinson disease (PD) and DBS treatment. Based on the work presented in Terman et al. (2002); Rubin and Terman (2004) and using complex network theory (Bassett and Bullmore 2006; Watts and Strogatz 1998; Bullmore and Sporns 2009), three areas of the BG: the globus pallidus (partes externa/interna) (GPe-GPi), the subthalamic nucleus (STN), and the thalamus, are modelled. We show, in accordance with the dopaminergic dysfunction during Parkinson’s disease, that different levels of striatal inhibition to BG areas change the behaviour from “normal” to “Parkinsonian”, i.e. switching from faithful transfer of information through the network to a state of disturbed information transmission. The model can also reproduce the action of DBS in the STN and illustrates how high-frequency stimulation (HFS) influences the whole network with respect to its computational features. Specifically, during DBS conditions, the model reveals a de-synchronisation or declustering of GPe and GPi activity which is projected to the thalamus. Defining and quantifying the response efficacy of thalamic activation during DBS, we deduce ranges of stimulation frequencies optimal for therapeutic success. Our model uses a significantly large number of neurons and connections (approx. two orders of magnitude larger compared to Terman et al. 2002; Rubin and Terman 2004; So et al. 2012) which makes the model more realistic and allows us to describe the behaviour of the neural network macroscopically. Although the detailed description is on the microscale (level of neurons), a macroscopic analysis is of interest, specifically the mean activity of interconnected neurons. The results of the macroscopic analysis suggest that a strong nonlinear response is obtained in the basal ganglia network, similar to resonating mechanical systems e.g. (Asadi et al. 2021; Bureau et al. 2014, 2013; Schilder et al. 2015). The following analysis proposes optimal frequencies above 130Hz, suggesting the investigation of DBS treatment with frequencies beyond 130Hz in animal experiments.
2 Mathematical modelling of basal ganglia and thalamic neurons
Our network model includes in addition to 3 areas of the basal ganglia: the subthalamic nucleus (STN), the globus pallidus internal (GPi) and external (GPe) also a part of the thalamus (THA), see Fig. 1. In this section we formulate the mathematical description for the neurons in each area of the basal ganglia (BG) and the thalamus.
In the model, the STN plays a key role as DBS is applied there to treat Parkinson’s disease. The dynamics of each STN, GPe, and GPi neuron are governed by a Hodgkin–Huxley formalism, and the current balance equation for the membrane potential reads (Terman et al. 2002; Bevan and Wilson 1999; Popovych and Tass 2019):
where C is the membrane capacity, \(V_i\) is the membrane potential of the i-th neuron, \(x_i\) denotes the gating variables n, h, r, and \([\text {Ca}^{2+}]_i\) is the intracellular concentration of calcium. For all basal ganglia areas the currents are described below: The leak currents \(I_{\text {LEAK}}=g_{\text {LEAK}}(V_i-E_{\text {LEAK}}),\) the potassium calcium and sodium currents are given by \(I_{\text {K}}=g_{\text {K}}n^4(V_i-E_{\text {K}})\), \(I_{\text {Ca}}=g_{\text {Ca}}s^2_{\infty }(V_i-E_{\text {Ca}})\) and \(I_{\text {Na}}=g_{\text {Na}}m^3_{\infty }h(V_i-E_{\text {Na}})\), while the low-threshold T-type calcium current for STN is \(I_{\text {T}}=g_{\text {T}}a^3_{\infty }b^2_{\infty }(V_i-E_{\text {Ca}})\). In the case of GPe, GPi neurons, the low-threshold calcium current has the form \(I_{\text {T}}=g_{\text {T}}a^3_{\infty }r(V_i-E_{\text {K}})\), reducing the bursting activity of the GPe relative to STN neurons. The current underlying the after-hyperpolarizing potential has the form \(I_{\text {AHP}}=g_{\text {AHP}}([\text {Ca}^{2+}]/(k_1+[\text {Ca}^{2+}])(V_i-E_{\text {K}})\).
The current \(I_{\text {DBS}}\) in eq. (1) models the deep brain stimulation of STN neurons and is set to the value 0 in the absence of DBS. The current \(I_{\text {app}}\) is applied to the STN, GPe, and GPi, but with different physiological meaning. In the case of STN neurons, \(I_{\text {app}}\) simulates the afferent synaptic input from the cortex (Terman et al. 2002), while in the case of GPe and GPi the \(I_{\text {app}}\) represents the incoming signal from the striatum with different levels of inhibition of GPe and GPi, respectively.
In eq. (2) the equilibrium state is \(x_{\infty }=x_{\infty }(V_i)=1/(1+e^{-(V_i-\theta _x)/\sigma _x})\) for \( x=n,m,h,a,r,s\), while for the equilibrium state of T-type current the following form is used: \(b_{\infty }(V_i)=1/(1+e^{(r_i-\theta _b)/\sigma _b})-1/(1+e^{-\theta _b/\sigma _b})\). The voltage-dependent timescale \(\tau _{x}\) has the form \(\tau _{x}(V_i)=[\tau _{x0}+\tau _{x1}/(1+e^{-(V_i-\theta _{\tau x})/\sigma _{\tau x}})]/ A_x\), for the STN neurons and \(\tau _{x}(V_i)=\tau \) for GPe and GPi neurons (Terman et al. 2002).
Figure 2 depicts the dynamics of one (uncoupled) STN neuron, firing at a frequency of 3Hz, while in all simulations we use \(I_{\text {app}}=4\)pA resulting in an STN activity of 6–7Hz (Bevan and Wilson 1999) (all values of parameters are given in table 1, see also Terman et al. 2002). When a negative current is applied for a short time, the neuron is hyperpolarized accordingly. Due to the presence of hyperpolarization-activated currents (HCN currents), a rebound burst occurs after the current injection.
According to experimental findings (Plenz and Kital 1999; Kita and Kitai 1991), GP neurons show similar ionic currents with respect to the STN cells, yet with different parameters. In the following we explain that the resulting dynamical properties are reproduced by our mathematical model. The main attributes are a spontaneous firing activity at a frequency of \(\approx \) 30Hz (Cooper and Stanford 2000; Kita and Kitai 1991), (see also (Terman et al. 2002) and the references there in), and a rebound response to a hyperpolarizing current (Terman et al. 2002; Kita and Kitai 1991; Cooper and Stanford 2000). Simulations of one GPe/GPi neuron are shown in Fig. 3, where the GP neuron fires at a frequency of 30Hz at rest shown in subfigure 3a. Figure 3b depicts the dynamics for small negative \(I_{\text {app}}\) resulting in intermittent bursting activity. A small depolarizing current, in turn, results in fast tonic discharges at nearly twice resting discharge frequency 3c. The entire tuning curve depicting responses to different levels of current injections is depicted in 3d.
2.1 Description of the basal ganglia synaptic connectivity
The coupling between the neurons is described by the synaptic current \(I_{\text {syn}}\). In the model, GPe and GPi neurons are connected through a Watts and Strogatz (WS) small-world topology (Watts and Strogatz 1998; Bullmore and Sporns 2009; Stam and Reijneveld 2007; Gafarov 2016; De Santos-Sierra et al. 2014; Netoff et al. 2004; Bertalan et al. 2017; Fang et al. 2017), where neurons nearby are densely connected; additionally, a small number of large distance connections exist. Following the experimental findings (Gouty-Colomer et al. 2018) which suggest sparse connections in the STN area, we chose a modified small-world topology with sparse internal connectivity. The detailed description of the network structure is given in sec. 2.3.
The small-world network is considered in the synaptic currents defined by the activation variable s, which are given by Laing and Chow (2002), Ermentrout and Terman (2012), Compte et al. (2000):
where H(V) is a smooth approximation of a step function, i.e. \(H(V)=1/(1+e^{-(V-\theta _x)/\sigma _x}\)).
The excitatory and inhibitory synaptic currents for the i-th neuron are given, respectively, by
with \(E_{\text {Glu}}=-10mV\), and
with \(E_{\text {GABA}}=-70mV\), where \(A_{ij}\) has the value 1 or 0, depending on whether neurons i and j are connected or not. The summation is taken over all presynaptic neurons.
In the case of STN neurons the \(I_{\text {syn}}\) current is given by the summation \(I_{\text {syn}}=I_{\text {STST}}+I_{\text {GPST}}\) and indicates the internal excitation between the STN neurons and the incoming inhibition from the GPe neurons, respectively. The excitatory glutaminergic connections within the STN are expressed by \(I_{\text {STST}}\) which follows eq. (5), while the inhibitory current \(I_{\text {GPST}}\) is given by eq. (6) and expresses the inhibition from the GPe area.
The synaptic current \(I_{\text {syn}}\) for the GPe region is defined by \(I_{\text {syn}}=I_{\text {GPeGPe}}+I_{\text {STGPe}}\), where the first term \(I_{\text {GPeGPe}}\) express the intra-layer inhibitory interaction of GPe neurons (i.e. follows eq. (6)), while \(I_{\text {STGPe}}\) describes excitation from STN neurons. For the GPi region the current \(I_{\text {syn}}\) is given by \(I_{\text {syn}}=I_{\text {GPiGPi}}+I_{\text {GPeGPi}}+I_{\text {STGPi}}\), where the first two terms \(I_{\text {GPiGPi}}\) and \(I_{\text {GPeGPi}}\) are inhibitory connections, connections from GPi to itself and from GPe to GPi, respectively, while \(I_{\text {STGPi}}\) describes excitations from STN neurons. The values of the parameters are given in table 1, see also tables 1, 2 of Terman et al. (2002).
2.2 Modelling and simulations of neurons in the thalamus
Modelling the basal ganglia network, the crucial behaviour of the model is the response of thalamic neurons to synaptic input from GPi neurons. The following section hence addresses this question. Here, a movement programme-associated sensory-motor cortex input to the thalamus is simulated by a repetitive periodic synaptic activation. This is modelled by 5 pA, 5 ms current injections, at 40 Hz (as depicted in Fig. 4), of the form Rubin and Terman (2004): \( I_{\text {SM}}=A_{\text {SM}}H(\sin (2\pi t/T_{\text {SM}})\cdot (1-H(\sin (2\pi (t+\delta _{\text {SM}})/T_{\text {SM}}))\), where H is the Heaviside function. The GPi neuronal input (resulting, in turn, from basal ganglia network activity) and the sensory-motor cortex input, determine the thalamic firing activity.
The mathematical description of thalamic neurons is given in the following equation
where C is the membrane capacity and \(V_i\) is the membrane potential of the i-th neuron. The leak current has the form \(I_{\text {LEAK}}=g_{\text {LEAK}}(V_i-E_{\text {LEAK}}),\) the other ionic currents, i.e. potassium and sodium, are given by \(I_{\text {K}}=g_{\text {K}}(0.75(1-h))^4(V_i-E_{\text {K}})\) and \(I_{\text {Na}}=g_{\text {Na}}m^3_{\infty }h(V_i-E_{\text {Na}})\), while the low-threshold T-type calcium current is described by \(I_{\text {T}}=g_{\text {T}}p^2_{\infty }r(V_i-E_{\text {T}})\). The gating variables h, r follow the differential equation with first- order kinetics as in eq. (2)
The equilibrium function has the form \(x_{\infty }(V_i)=1/(1+e^{-(V_i-\theta _x)/\sigma _x})\) for \(x=r,h,m,p\), and the voltage-dependent timescale for the gating variable r is \(\tau _{r}=28+(1+e^{-(V_i+25)/10.5})\), while for h is defined as \(\tau _{h}=1/(a_h+b_h)\) with \(a_h=0.128e^{-(V_i+46)/18}\) and \(b_{h}=4/(1+e^{-(V_i+23)/5})\). The current \(I_{\text {GPTH}}\) represents the inhibition of the thalamus by the GPi and has the form of eq. (6). For a detailed description of thalamic neurons and for the arithmetic values of parameters, see Rubin and Terman (2004) and table 2.
As Fig. 4 shows, within an isolated cortico-thalamic interaction, the thalamus follows cortical input absolutely faithfully. The next important question is how the basal ganglia network input will modify this strong thalamocortical interplay.
2.3 Network structure
For modelling and analysis \(N_{\text {STN}}=500\) STN neurons, \(N_{\text {GPe}}=500\) GPe neurons, \(N_{\text {GPi}}=500\) GPi neurons were used, while the thalamus was represented by \(N_{\text {THA}}=200\) neurons. Following (Watts and Strogatz 1998; Bassett and Bullmore 2006; Bullmore and Sporns 2009; Stam and Reijneveld 2007; Spiliotis and Siettos 2011; Shefi et al. 2002), the GPe/GPi layers were modelled as separate small-world networks, i.e. the connections of neurons follow a small-world topology. In such small-world complex networks (Newman 2003), not only does each neuron (node) in the network interact with its k nearest neighbours, but there are also a few randomly chosen remote connections (Watts and Strogatz 1998). For the construction of nuclei networks (GPe and GPi) the value \(k=20\) was used, while remote connections were created according to Watts and Strogatz (1998) with \(p=0.005\). Figure 5 shows a characteristic snapshot of the network. For the STN, we chose a modified approach with reduced internal connectivity, as it is already mentioned in sec. 2.1, matching experimental findings which suggest sparse connections (Gouty-Colomer et al. 2018; Ammari et al. 2010). Specifically, in these experimental findings (Gouty-Colomer et al. 2018), 80% of STN neurons do not have any connections and the remaining 20% of the STN neurons form connections with each other. For these connecting STN neurons, measuring the distance of axonal endings with collaterals (within STN) reveals that roughly 30% of all synapses lie within a 200 \(\mu \)m radius and another 45% within the 200-400 \(\mu \)m radius. The other 20% are contacts which occur farther away, i.e. \(>500\) \(\mu \)m (Gouty-Colomer et al. 2018). In this sense, the 20% of neurons, which form STN connections, have both local and remote connections similar to small-world property (i.e. 20% of neurons only showing an average of 25 connections each, see Gouty-Colomer et al. 2018). Finally, the case of denser connectivity between STN neuron is also studied in section 4.5 and the results are shown in Fig. 10.
The coupling between the nuclei is modelled in the following way: Each GPe neuron is connected to one STN neuron and vice versa. Furthermore, each STN neuron projects to one GPi neuron. In addition, the GPe also sends inhibitory signals to the GPi according to a small-world connectivity (one GPe neuron inhibits locally GPi neighbours but also inhibits a few randomly selected remote GPi neurons; in total, one GPe is connected on average with 20 GPi neurons), see Fig. 5. Each thalamic neuron, in turn, is receiving inhibitory inputs from three GPi neurons. Since the main interest of this study is the impact of the basal ganglia output activity on thalamic response to sensorimotor signals (shaping the thalamocortical interplay), this model does not include intrathalamic connections. The spatial arrangement in Fig. 5 is made according to the indices of the neurons; their detailed position in space is not relevant for the model (we use a coordinate system only for visualisation). The ring satisfies periodic boundaries conditions.
3 Macroscopic description of basal ganglia dynamics
3.1 Synchronisation analysis of basal ganglia dynamics
In the following, a nonlinear dynamical system of the form
is considered, where f is a nonlinear vector field. For a periodic (or oscillatory) solution the property
holds with period T, where x(t) defines a periodic orbit in the phase space. If the periodic orbit shows normal hyperbolicity (Izhikevich and Kuramoto 2004), the system of eqs. (9) can be transformed into an angle or phase equation for the angle \(\theta \), i.e.
where \(\omega =2\pi /T\) is the natural frequency of the oscillator.
Coupled oscillators described as system of differential equations (9) are of particular interest in neuroscience. Under certain conditions, identical or nearly identical interacting oscillators converge to a common frequency and synchronise (Strogatz 2001).
The weak coupling of n oscillators (Izhikevich and Kuramoto 2004; Kuelbs et al. 2020) is formulated by
where \(\epsilon \ll 1\) represents the coupling strength. The summation is over the coupling between j and i, namely over the adjacent of the i-th oscillator (connections which project to the i-th oscillator). Similarly, there is a transformation which allows to express the dynamics in phase variables, i.e.
where \(\theta _i \in [0,2\pi ]\). In the presence of external stimulation (e.g. DBS), eq. (12) can be generalised (Monga et al. 2019; Kuelbs et al. 2020), by including a time-periodic function on the rhs. The weakly coupled system in terms of the phase variable \(\theta \), i.e. eq. (14), shows a phase-locked solution if there is a constant integer matrix \(K_{n-1,n}\), such that \(K\cdot \theta =c\) with \(c_k>0\) (Izhikevich and Kuramoto 2004). Furthermore, the coupled oscillators are synchronised (in-phase) when
Following the approach of Kuramoto (Arenas et al. 2008; Kuramoto 1984b, a; Bertalan et al. 2017) which is applied when oscillators are near a supercritical Andronov–Hopf bifurcation, in case of fully connected network eq. (14), we obtain
while with arbitrary connectivity i.e. a complex network topology
where the coefficient \(A_{ij}\in \{0,1\}\) is derived from the network adjacency matrix (Bertalan et al. 2017). The emerging macroscopic dynamics (in terms of phase) can be obtained by taking the mean value of all phase populations (in exponential form, \(e^{i\theta }\)). This defines the synchronisation index r, i.e. (Kuramoto 1984b; Strogatz 2001; Kuelbs et al. 2020)
This index is used in the following sections as a measure describing the level of synchronised patterns within the GPi.
The synchronisation index acts as macroscopic variable (order parameter) with range \(r\in [0,1]\). In the case of perfect synchronisation, i.e. eq. (15), the index can be written as
while the case of \(r \rightarrow 0\) corresponds to incoherent phase dynamics. The phase \(\theta _k(t)\) used in eq. (18) of the k-th neuron can be approximated linearly according to the following equation
where \(t_n\) corresponds to n-th firing time of the k-th neuron and \(t\in [t_n, t_{n+1}]\). Another option to compute the phase without assuming a linear dependence of the angle \(\theta \) on time, see eq. (20) is the Hilbert transform, applied by Gabor in Gabor (1946). For a given function \(x=x(t)\), the Hilbert transform is defined as
with P denoting the Cauchy principal value. The complex signal \(z=z(t)=x(t)+iX(t)=A(t)e^{i\theta (t)}\) is defined in order to extract the phase, where \(A, \theta \) is the amplitude and the phase of the complex signal z, respectively. The instantaneous phase is defined as:
The synchronisation index eq. (18) will be computed in the next section to characterise the network activity for the cases of normal, Parkinsonian, and DBS treatment.
An alternative macroscopic quantity which can be used to estimate neuronal activity is the mean synaptic activity index l. It is defined by
where \(s_i\) is the synaptic variable as it is defined from eq. (4). Coupling theoretical modelling tasks with experimental data, the synchronisation property in essence reflects local field potentials (LFP). The flow of the extracellular current generating the LFP is represented by the summed up postsynaptic potentials from local cell groups (similar to eq. (23), see Buzsáki 2004; Popovych and Tass 2019; Manos et al. 2018; Popovych and Tass 2018). In the next section this mean synaptic activity will be used to measure the transfer of neuronal activity (current flow information) from the basal ganglia output, i.e. from GPi to thalamus.
4 Simulating different functional states of the basal ganglia network
Three cases of dynamic network behaviours were studied. In all cases the thalamus receives a periodic sensorimotor input (simulated by periodic 5 pA current injections Rubin and Terman 2004) which represents the signal for the initiation of movement. In addition, a continuous input current to STN is applied in order to obtain rhythmic activity (simulating afferent synaptic input from cortex Terman et al. 2002). During the first, normal case, see Fig. 1b, basal ganglia network activity allows a relatively faithful response of the thalamic neurons, closely following the sensorimotor cortex input. The second case considered is the Parkinsonian state, depicted schematically in Fig. (1c), with resulting overall increase in inhibitory projections to the thalamus. The third case considers the situation of therapeutic intervention with deep brain stimulation (DBS) to STN, which is simulated as a high-frequency current injection into STN neurons, resulting in switching the entire network dynamics and restoring functional thalamic output close to normal behaviour.
In order to capture the transitions and to measure qualitatively the impact of the network dynamics, we additionally defined macroscopic variables or observables such as the synchronisation index and the average synaptic activity, focusing on GPi as output region of the basal ganglia network.
All simulations are made with MATLAB 2020a using the solver ode23, a adaptive time step integration Runge–Kutta scheme. The order of the method is three, and the default relative and absolute tolerances in MATLAB \([10^{-6},10^{-6}]\) are used.
4.1 Modelling the normal state
Parameters were tuned to simulate normal-healthy conditions as shown in Fig. 1b. Inhibitory signals from the striatum to the GPe are described by the current \(I_{\text {app}}=5 \text {pA}\), whereas from striatum to GPi the current has the value \(I_{\text {app}}=4\)pA. Figure 6 depicts the time-dependent activity of all nuclei under normal conditions. The STN neurons fire irregularly at around 4–7 Hz, i.e. close to values reported in the literature (cf. Bevan and Wilson (1999) Rubin et al. 2012), see Fig. 6a. The GPe activity is plotted in Fig. 6b. It is characterised by high-frequency irregular bursting firing with individual clustering of action potential series (which is not visible on population scale, but only on the level of a single neuron as shown in the right column). Overall, there is little correlation between STN and GPe activity. Similar dynamics appear in GPi, which, in turn, also shows high-frequency firing. Looking more closely, there is, however, an underlying weak rhythmicity at approx. 12 Hz, i.e. in the \(\beta \)-band, see Fig. 6, and also 12(a). The resulting time-dependent activity plot of all neurons in the thalamus is shown in Fig. 6d, demonstrating that under these conditions, the thalamus faithfully responds to sensorimotor input (Fig. 6e), without however being in a perfect phase-locked state as in the case of single neuron, Fig. 4. Accordingly, the thalamic response efficacy R (a macroscopic variable which measures the response of the thalamic neurons to sensorimotor input) is approximately 0.5 under normal conditions. The response efficacy R ranges in the interval [0,1], while the value 1 corresponds to an activation of the whole thalamus (exact definition of R, see sec. 6 and Fig. 14). In conclusion, the GABAergic synaptic output from the GPi to the thalamus converts thalamic responses (otherwise phase locked to the cortical input) to more loosely firing of the thalamic neurons, without subduing it altogether (thalamic activity \(R=0.5\)) as in the Parkinsonian state (see Fig. 7).
4.2 Modelling the Parkinsonian state
As outlined above, in Parkinson’s disease, a degeneration of nigrostriatal dopaminergic neurons, accompanied by a reduction in the number of dendritic spines of striatal medium spiny neurons (Gagnon et al. 2017; Fiore et al. 2016), leads to a loss of dopamine in the striatum. The resulting overall reduction of D1/D2 receptor-mediated activity affects the direct/indirect pathway functionality. In the direct pathway, the (D1) receptors malfunction results in disfacilitation of striatal projection neurons and as a consequence a reduced inhibition of the GPi neurons. Thus, the disinhibited GPi increases its neuronal activity, sending higher levels of inhibition to the thalamus.
In the indirect pathway, loss of D2-receptor activation will disinhibit striatal projection neurons. These, in turn, now decrease the activity of the GPe, to which they project. By inhibiting the GPe, the activation of STN is increased, and the overactive STN will enhance the neuronal activity in the GPi which again leads to even more pronounced thalamic inhibition, see Fig. 1c.
Consistent with this concept, the model assumes a decrease in the level of inhibition of GPi neurons \(I_{\text {app}}\). This is modelled by increasing depolarising current from 4 to 8 pA (corresponding to disinhibition). At the same time, the model assumes an increase in the level of inhibition of GPe neurons (therefore, the depolarising current is decreased from 5 to 3pA). Figure 7 shows the dynamics of the network under Parkinsonian conditions. The time-dependent activity of all neurons in STN is changed compared to the normal state: Firing becomes more regular and occurs at higher frequency (approx. 11 Hz). For the GPe, Fig. 7b, the consequence is a lower burst frequency, almost following the STN activation (see insets). Figure 7c shows that in GPi, the altered activity is translated to high-frequency bursting, with an underlying enhanced rhythmicity in \(\beta \)-activity (13-15 Hz) compared to the normal condition, see also Fig. 12b. As a result, inhibition to the thalamus strongly increases, and the thalamus is no longer able to transmit signals in response to sensorimotor input, as its firing becomes very sparse.
4.3 Modelling the effects of DBS
The network structure and the conductances were kept invariant with respect to Parkinsonian conditions. DBS treatment was simulated by a high-frequency current of 184 Hz, applied to all STN neurons (this value resulted after the analysis in section 6 and is suggested as candidate for an optimal DBS frequency in experiments). The high-frequency current is modelled as periodic short pulses of the form (Rubin and Terman 2004)
In all three basal ganglia nuclei and thalamus, STN-DBS induces dramatic alternations in firing dynamics. Regarding the STN, neurons follow the strong DBS signal and fire tonically at stimulation frequency, see Fig. 8a. In GPe, this results in regular bursting of neurons at about 90 Hz, i.e. at 3x higher frequency than normal and 4x higher frequency than under Parkinsonian conditions, respectively. In the GPi, in turn, the increased activity in GPe nearly normalises the firing to slightly irregular firing at 30 Hz (compare insets in this figure and 6) and thus strongly reduces the firing frequency compared to the Parkinsonian state (which ranged around 50–80 Hz) see Fig. 8. In summary, there is an overall loss of synchronisation between STN, GPE, and GPi. Simultaneously the STN and GPe nuclei now fire tonically. With the reduction of the high-frequency tonic firing of the GPi, the thalamus is disinhibited and resumes firing, following the sensorimotor input relatively closely.
4.4 STN-GPi interaction: synchronisation changes and local travelling waves formation
Different spatio-temporal patterns of activation in STN on the one hand, and GPi on the other, can be observed under normal, Parkinsonian, and DBS conditions. Figure 9 depicts the whole neuronal activity of STN and GPi regions in the case of sparse connectivity between STN neurons. In normal case (a) and (b), both areas STN and GPi show irregular patterns. By contrast, in the Parkinsonian case (c) and (d), the regions show higher synchronous activation. Figure 9c shows sparse, but synchronous firing in STN in the Parkinsonian case, with local travelling waves. These waves are restricted to a few neurons, demonstrating clustered organisation of neuronal activity. An interpretation for this is that in the Parkinsonian state, the reduced inhibition from GPe to the STN leads to a higher spatial synchronisation between these nuclei, resulting in a different pattern formation with respect to normal case 9a, b.
Which consequence does this have for the activity in the GPi?
This nucleus receives both input from the STN (indirect pathway) and from the striatum (direct pathway), and hence one would expect two competing activity patterns. Under Parkinsonian conditions, in turn, the inhibitory activity of the direct pathway is reduced. Hence the activity in the GPi is more pronounced, which shows series of quasi-synchronous activations of GPi neurons, with local travelling waves forming small clusters.
Under DBS conditions, the massive synchronisation in the STN (abolishing all travelling waves and imposing a 184 Hz rhythm on nearly all neurons) is mirrored in a closely matching synchronised activity pattern in the GPi, but at lower frequency of around 60Hz, again overriding all other activity. This presumably results in a nearly tonic activation of inhibitory neurons, to which the thalamic neurons desensitise, as will be discussed below.
4.5 Effects of dense STN connectivity
Since the degree of connectivity among STN neurons is still under debate (Amadeus Steiner et al. 2019), in a next step also a higher connectivity was considered. The increment of connectivity, i.e. higher number of connections between STN neurons, leads to different spatio-temporal pattern. Figure 10 shows the normal, Parkinsonian, and DBS states. The qualitative difference here is the strong connectivity within STN. Each neuron now has mutual connections with other neurons, with a mean number of connections equal to 20. The STN structure is described as a small-world topology. Under normal and Parkinsonian conditions, the basal ganglia network depicts strongly correlating travelling waves, which propagate and are more variable in the normal state. The travelling wave represents a propagation of similar activity levels along one direction, i.e. as activity peaks in neighbouring neurons in a co-moving frame, preserving the activity shape. Also with this much higher intra-STN connectivity, the Parkinsonian state is characterised by increased synchronisation in GPi. Under DBS, rhythmic activity of STN in turn is reflected in less synchronisation in GPI. Thus, the model is rather robust with respect to GPi output, within a wide range of STN connectivity.
5 Synchronisation and synaptic activity indices characterise transitions of network dynamics
In order to shed light on the transitions between the different states, and to examine the appearance of distinct patterns in our model, we analysed both the level of synchronisation within the GPi (computing the synchronisation index r) and the levels of synaptic GABAergic projection activity from the GPi to the thalamus (defining GPi activity index l).
The theory of phase synchronisation (Izhikevich and Kuramoto 2004; Ermentrout and Terman 2012; Tass 1999; Pikovsky et al. 2001; Kuramoto 1984b) described in section 3.1, allows us to characterise and analyse different attributes and the dynamics which results from the mathematical model. The main observables are the synchronisation index r defined in eq. (18) and the mean GPi synaptic activity index l defined in eq. (23) and depicted in Figs. 11 and 12. We further performed a Fourier analysis of the GPi mean synaptic activity index l, in order to measure the dominant frequencies (rhythms), and we determined its slope of the power spectrum to estimate the distribution of frequencies in the power spectrum (see Fig. 13).
5.1 Macroscopic dynamics in the normal state
Considering the healthy, normal state, as Fig. 11a shows, synchronisation within the GPi is generally low (around 0.5). The synchronisation index shows oscillatory behaviour (at around 12 Hz). The local maxima of the synchronisation index thus correspond to increased synchronous activation of GPI neurons in the low \(\beta \) range. The GPi activity results in an overall periodic, but low-magnitude inhibitory drive to the thalamus, see 12a. Corresponding to this behaviour, the mean synaptic GPi activity phasically oscillates around 0.2 (even if the microscopic behaviour of neurons forms bursting clusters) see, Fig. 12a. This rhythm lies within the lower frequencies (i.e. \(\alpha \) and lower \(\beta \) band), see Fig. 13a. The spectrum analysis reveals a first peak at 7 Hz and a higher peak at around 11 Hz, i.e. activity in the lower \(\beta \) band. This rhythmic peak is centred relatively narrowly around this dominating frequencies, since the slope of linear approximation is relatively steep (around \(-3.6\) Hz/sec), meaning that the higher frequencies do not contribute in the magnitude of the power spectrum.
5.2 Macroscopic dynamics in the Parkinsonian state
Under Parkinsonian conditions, the GPi neurons are generally much more active than in the normal state. The synchronisation is very high (close to 1). This index reflects prolonged intervals of high activity, interrupted only briefly by small decrements in synchronisation levels. This indicates that neurons are synchronised during the prolonged and accentuated bursting in GPi in the Parkinsonian state. As a result, the synaptic, inhibitory projection to the thalamus remains periodic, but predominantly strong, and is again only slightly reduced during short burst intervals, see Fig. 12b. The mean GPi synaptic activity (fluctuating around 0.7) is commensurate with this high synchronisation, which is supported by prolonged bursting of GPi neurons, only periodically interrupted by brief periods of reduced activity. The Fourier analysis reveals a dominant \(\alpha \) /low\(\beta \) activity with a strong peak at 10Hz, and a secondary at 18Hz, i.e. higher \(\beta \) rhythm than the normal conditions. Importantly, the frequency spectrum of this bursting is much broader, as indicated by the corresponding linear approximation: The slope in Fig. 13b is increasing to \(-2.5\) Hz/sec, which reflects that the higher frequencies contribute more (with respect to the normal case) to the power spectrum.
5.3 Macroscopic dynamics during DBS
Simulating DBS, the macroscopic activities (as the order parameters r and l show) change dramatically (depending on a specific range of frequencies), introducing a state dissimilar to both normal or Parkinsonian states. The synchronisation index in Fig. 11c within the GPi shows a dynamic development over time, starting at values of \(\approx \) 0.9, and then decreasing to \(\approx \) 0.4, during the first 500 ms of tonic irregular firing in the GPi. After this period, the system adapts and oscillates at high frequencies with values depending on the applied frequency (e.g. between 0.8 and close to 1 in case of 184 Hz, between 0.6 and 0.8 in the case of 210 Hz). Similar dynamics emerge regarding the mean GPi synaptic activity, with a first transient period of 600 ms where the synaptic drive is decreased (to \(\approx \) 0.2), and a second period (\(t>600\) ms) where high synaptic activity is maintained (with fluctuations around a mean value which depend on the applied frequency (e.g. \(\approx \) 0.6 in case of 184 Hz, \(\approx \) 0.5 in case of 210 Hz) i.e. higher than under normal conditions and lower than under Parkinsonian. Surprisingly, this would appear to impose a strong inhibitory drive to the thalamus, see Fig. 12c. This, however, differs from the Parkinsonian condition as it is tonic high-frequency and not repeating bursting blocks appearing at \(\beta \) frequency. This qualitative change in rhythmicity is crucial: While a relatively high inhibitory tonic drive seems to suggest strong thalamic inactivation, the contrary is the case. Due to the fact that the input to the thalamus is high-frequency tonic (peak power at 184 Hz), the GABAergic synapses are deactivated (see eqs. (4)). Indeed, the frequency analysis of the GPi activity shows a main peak seeming to resonate with the external DBS stimulus (around 60Hz) and the corresponding harmonics. Importantly, the peaks in the lower (normal state) and higher (Parkinsonian state) \(\beta \) band disappear, while the slope changes back to a value around −3.1 Hz/sec, close to the normal case. Although this slope is close to normal conditions, it now reflects a narrow high-frequency band and not unclustered activity.
Taken together, the findings underscore the importance of synchronous regular and brief, clustered burst firing in the GPi for successful inhibition in the thalamus, which normally takes place at \(\beta \) frequency. Although DBS thus does not restore normal rhythmicity in the GPi, the thalamic activity is disinhibited by the loss of synchronisation and clustering within the GPi.
6 DBS efficiency depends on stimulation frequency
As the previous considerations show, the macroscopic order parameters r and l, which are defined in eq. (18) and (23), respectively, allow us to describe the effects of DBS stimulation and to compare the different states (normal, Parkinsonian, and DBS).
One critical parameter quantifying the effect of DBS simulation is the thalamic response with respect to the sensorimotor cortical signals. Faithful thalamic activation will send strong excitatory signals to cortex, see Fig. 1a, alleviating in this way the Parkinsonian symptoms (Guo et al. 2008). Naturally, under normal conditions the thalamic response should neither be completely uncoupled (with respect to the sensorimotor cortical signal), as under Parkinsonian conditions (see Fig. 7d), nor completely coupled as in an isolated cortico-thalamic system (Fig. 4). Indeed, under normal conditions, as shown in Fig. 6d, the thalamus follows the cortical input closely, but not in absolute synchrony.
In order to quantify the thalamic response to cortical input, under normal conditions or DBS with various stimulation frequencies, we define the response efficacy R of thalamic neurons as a mean value of the fraction of activated thalamic neurons (per stimulus) during simulated cortical activation. This simulation, as already discussed in Fig. 4, comprised sensorimotor current injections of length \(\delta \) (5 ms) defined by the interval \([t,t+\delta ]\), with frequency of 40 Hz. The mathematical formulation of response efficacy R is thus defined as:
where \(a_i(t)\) is the proportion of activated neurons, i.e. the number of activated neurons within the time interval of \([t,t+2\delta ]\) divided by the total number of thalamic neurons \(N_{\text {THA}}=200\). The summation is taken over the number of intervals \(N_{\text {int}}\), for times \(500<t<1500\) ms, (i.e. \(N_{\text {int}}=100\)). Under normal conditions, as expected, R is approximately 0.5, i.e. suggesting a coupling to cortical input at a value offering the broadest dynamic range.
Under DBS at various frequencies, R shows a nonlinear behaviour, starting at 0 for low DBS frequencies (around 50 Hz) to values \(>0\) with small peaks for frequencies around 70 and 130 Hz, see Fig 14. Very prominent peaks are found around 184 Hz, 210 Hz, and 244 Hz. Thus, there are three dominant optima at where the DBS maximises the thalamic response close to normal values of R. As Fig. 8 thus shows, the frequency of DBS is critical for thalamic firing outcome. At low frequencies (i.e. 50–150 Hz in Fig. 14), only few frequency bands can be found at which \(R > 0\), while at this low range, the effects on thalamic activation are only transient, e.g. 50 Hz and 150 Hz (see Fig. 14). Exceptions are the small peaks at 70 and 130 Hz; Fig. 14; reasonably stable responses can thus be seen at 130 Hz. Only from 160 Hz onward, stable thalamic firing can be achieved by DBS. Interestingly, this is very similar to experimental findings in hemi-Parkinsonian rats: low-frequency DBS stimulation up to 75 Hz actually results only in transient effects, while permanent reductions of circling behaviour (the Parkinsonian equivalent in this animal model) were only achieved at DBS frequencies > 130 Hz (So et al. 2017).
Computing the Shannon entropy for the macroscopic variables r, l we confirm the optimal DBS frequency. As suggested by Dorval et al. (2008), Arle et al. (2018), Deco et al. (2012), an optimised DBS frequency can be achieved by regularisation of the whole basal ganglia activity (Dorval et al. 2008), i.e. minimisation of the entropy and thus a more ordered state. The Shannon entropy E is defined by
where \(x_i\) expresses the macroscopic variable of interest (here \(x=r, l\), as defined in equation (18) and (23)). The computation of the entropy with respect to the frequency was performed for the first 1500 ms including also the transient period of the first 500ms to be able to compare the initial response in experiments once the DBS is switched on (So et al. 2017; Mottaghi et al. 2020). The results are depicted in Fig. 15. From this figure it can be derived that the optimal frequencies for DBS are at 184 Hz and 210 Hz, where the entropy of the macroscopic variables (order parameters) r, l reach a local minima (in accordance with the thalamic response efficacy R, cf. Fig. 14, circle and square markers). A similar minima can be found when increasing f further at 244 Hz, which would require more energy for the DBS stimulation, likely raising the possibility of energy-dependent side effects.
7 Discussion and conclusions
A large-scale basal ganglia computational model has been used to study network dynamics in movement disorders. The network model consists of 1700 interactive neurons with approximately 30000 connections (representing the microscopic level). The connectivity within the areas follows small-world topology where nearby neurons are densely connected; additionally, a small number of long-distance connections exist. The resulting structure combines properties of both regular and random topology, e.g. high clustering and at the same time short path length. Small-world structure enhances signal-propagation speed and synchronizability among the areas of the neuronal network.
For different values of parameters, the model reproduces phase transitions for normal, Parkinsonian, and DBS, in the emergent dynamics which are captured with suitable macroscopic order parameters for example the change from oscillatory to bursting behaviour of mean synaptic activity in Fig. 12a, b. Importantly, transitions produced by the model are consistent with the physiology and experimental observations of aberrant functionality of direct and indirect pathways.
Experimental findings of animal studies regarding the impact of STN-DBS on GPe, GPi, and thalamic firing could be approximated qualitatively. Specifically, in the transition from normal to Parkinsonian state shown in Figs. 6 and 7, the model alters the dynamics of GPe and GPi neurons due to varying levels of striatal inhibition. Remarkably, similar alternations in GPe/GPi dynamics were observed in Galvan and Wichmann (2008) (Fig. 2 therein) in monkeys and mice treated with MPTP (methyl-4-phenyl-1,2,3,6-tetrahydropyridine) (Galvan and Wichmann 2008; Ogawa et al. 1985). Specifically, in GPe single-cell recordings of monkey, the firing becomes sparser. This is indeed the case also in our model, see Fig. 7b. Concerning the GPi single cell activity, the firing rate increases strongly in the Parkinsonian state which is also reflected in our model, see Fig. 7c. In the STN area, the monkey single cell activity depicts bursting behaviour, while in our model burst firing only occurs sparsely, or not at all. However, the overall frequency increases similar to monkeys MPTP.
Furthermore, under DBS, in the model, the GPi and GPe fire tonically at high frequency, following STN firing at stimulation frequency (184 Hz), abolishing synchronised \(\beta \) activity. This corresponds to the findings of McConnell et al. (2012), who showed that STN stimulation results in a sharp power peak at the same frequency also in GPe, and with findings of Wang et al. (2018) where DBS disrupts pallidal beta oscillations and the cortical coherence in Parkinson disease.
Regarding the \(\beta \)-band hypersynchrony hypothesis of Parkinson’s disease (Kühn et al. 2008), the model thus importantly shows that DBS actually abolishes \(\beta \)-band synchrony, which also in the model is prevalent in GP without stimulation. Hence, the model faithfully replicates what is known from clinical and animal model studies both on the cellular and on the network activity pattern levels. The model produces clusters of local travelling waves, more pronounced in the Parkinsonian state. Remarkably, in the experimental findings of Cagnan et al. (2015), it has been found that the analysis of local field potential recordings from the subthalamic nucleus and globus pallidus of patients with Parkinson’s disease shows beta-band propagation waves within the globus pallidus.
Further, the model predicts that a faster firing GPi projections under DBS will result in regular firing of the thalamus, which in turn will essentially follow the sensorimotor input—while under PD conditions, in fact thalamic firing was very sparse and irregular, and did not faithfully mirror cortical signals. These results are supported by studies in a Parkinsonian animal model in rhesus monkeys (Xu et al. 2008), where STN-DBS (as in Fig. 8) produces a change in the pattern and periodicity of neuronal activity in the basal ganglia thalamic network, resulting in a regular, higher-frequency firing pattern in the thalamus.
Moreover, macroscopic properties can be derived from our model. Since these essentially mirror local field potential activity, these properties can be used to test predictive modelling in future studies. For example, the model predicts a power law power spectrum, for the mean synaptic activity (see Fig. 13) with variable critical exponent a. Similar power law is reported in He et al. (2010) in the dynamics of visual cortex, hippocampus, and cerebellum with similar changes in the critical exponent a. They explained these differences due to the different areas of activation during the experiments. A power law behaviour (in basal ganglia areas and cortex) with deviations also is reported in West et al. (2018). In Huang et al. (2020) using local field potentials for 12 Parkinson’s disease patients, a power law activity in STN and cortex is reported. The STN activity shows strong variations of the exponent a during awake and loss of consciences state. In our case, this exponent is the slope of the linear approximation in Fig. 13. Our analysis shows that the slope can be distinguished between normal Parkinsonian and DBS state. Further, the model also allows to extract variables such as a macroscopic synchronisation and GPi synaptic activity indices which allow to predict the transitions during DBS application, which might pave the way for feedback control of DBS in the future (Popovych and Tass 2018, 2019; Manos et al. 2018).
Beyond this, the detailed analysis of the macroscopic parameters and derived values such as entropy also allow for an optimisation of DBS. As an outcome measure, the response efficacy of thalamic neurons reflects the degree of thalamic activity correlating to cortical input, which under normal conditions is \(\approx \) 0.5. Critically, raising stimulation frequency beginning with 50 Hz, the first prominent entropy local minima of the two mesoscopic parameters precisely coincide with the first peak of the response efficacy R close to normal values (\(\approx \) 0.23) at 184 Hz and more efficient frequencies over the 200Hz. These frequencies are not identical to the most commonly used in clinical settings i.e. 130 Hz, while in this frequency range, the model predicts sub-therapeutic action (see Fig. 14). In addition, one can argue that our modelling study does not fully replicate clinical findings, which show that at frequencies \(>185\) Hz, side effects like dyskinesias seem to become more problematic (Dayal et al. 2017; Karl et al. 2019). Contrariwise, the match between the experimental findings in hemi-Parkinsonian rats (So et al. 2012, 2017; Mottaghi et al. 2020) and the response efficacy R is very interesting. In So et al. (2012) the authors describe a sustained depression of circling to values smaller than 5 turns per minute only at frequencies > 75 Hz. This incidentally matches well with our findings, which show peaks in the response efficacy R at 70 Hz and 130 Hz, and then at peaks of 180, 210, and 245 Hz, in close approximation to the effective frequencies reported in So et al. (2012, 2017), i.e. 185 and 260 Hz. Remarkably, similar frequency-dependent behaviour is reported in the recent study of Mottaghi et al. (2020), where strong impact of stimulation frequency on the induced rotation (caused by DBS) was found in the range of 250Hz.
Our mathematical model extends previous computational work (So et al. 2012; Dorval et al. 2010) by using a larger number of neurons for each basal ganglia areas and complex small-world connectivity. Additionally, we propose different approximations in the DBS frequency. The frequency analysis of Dorval et al. (2010) differs from our investigation since they focus on perturbations around 130 Hz (Fig. 1 and 2 in Dorval et al. 2010). In So et al. (2012) the activation patterns of local cells and fibres passage are studied with respect to the fidelity of thalamus. Our model suggests that the basal ganglia network behaviour to DBS stimulates frequency has a strong nonlinear response similar to resonating mechanical systems with optimal frequencies above 130Hz, suggesting the investigation of DBS treatment beyond the 130Hz. In summary, the mathematical model and the analysis is considered a powerful tool to explore the effects of parameter-dependent changes of DBS and to optimise the medical treatment.
8 Outlook
One major future topic will be the study of network topological variations and the impact on the emergent dynamics. Thus, the functional effects of neuroanatomical changes observed in Parkinson’s disease (Prakash et al. 2016), such as massive decreases of dendritic length of medium spiny neurons (Stephens et al. 2005), could be modelled. Considering movement disorders beyond Parkinson’s disease, the model could also be extended to further elements of the basal ganglia, to gauge the effect of alterations in cortico-striatal communication such as those observed in dystonic hamsters (Köhling et al. 2004), where alternations in excitability in dystonic tissue were described to be related to both changes in intrinsic neuronal properties and presynaptic release probability at glutamatergic synapses. Furthermore, our analysis can be extended beyond the initial interval of the first 1.5 sec (which contains also transient effects). The analysis of the long-time behaviour under DBS including the effects of neuromodulators and the changes in synapses functionality (network plasticity) (Marschler et al. 2014; Morrison et al. 2008; Droste et al. 2013), is an interesting subject for future studies.
In future work, it will be important also to investigate the parameter dependence of network dynamics, using numerical bifurcation tools for multiscale problems (Spiliotis and Siettos 2011; Marschler et al. 2014; Moon et al. 2015; Schmidt et al. 2018), for an in-depth understanding of the functional network changes occurring in movement disorders.
References
Amadeus Steiner L, Barreda Tomás F, Planert H, Alle H, Vida I, Geiger J (2019) Connectivity and dynamics underlying synaptic control of the subthalamic nucleus. J Neurosci 39(13):2470–2481
Ammari R, Lopez C, Bioulac B, Garcia L, Hammond C (2010) Subthalamic nucleus evokes similar long lasting glutamatergic excitations in pallidal, entopeduncular and nigral neurons in the basal ganglia slice. Neuroscience 166(3):808–818
Arenas A, Díaz-Guilera A, Kurths J, Moreno Y, Zhou C (2008) Synchronization in complex networks. Phys Rep 469(3):93–153
Arle J, Mei L, Carlson K, Shils J (2018) Theoretical effect of dbs on axonal fibers of passage: Firing rates, entropy, and information content. Stereot Funct Neurosurg 96(1):1–12
Asadi K, Yeom J, Cho H (2021) Strong internal resonance in a nonlinear, asymmetric microbeam resonator. Microsys Nanoeng7(1)
Ashkan K, Rogers P, Bergman H, Ughratdar I (2017) Insights into the mechanisms of deep brain stimulation. Nature Rev Neurol 13(9):548–554
Bassett D, Bullmore E (2006) Small-world brain networks. Neuroscientist 12(6):512–523
Bertalan T, Wu Y, Laing C, Gear C, Kevrekidis I (2017) Coarse-grained descriptions of dynamics for networks with both intrinsic and structural heterogeneities. Front Comput Neurosci 11
Bevan M, Magill P, Terman D, Bolam J, Wilson C (2002) Move to the rhythm: oscillations in the subthalamic nucleus-external globus pallidus network. Trends Neurosci 25(10):525–531
Bevan M, Wilson C (1999) Mechanisms underlying spontaneous oscillation and rhythmic firing in rat subthalamic neurons. J Neurosci 19(17):7617–7628
Bullmore E, Sporns O (2009) Complex brain networks: graph theoretical analysis of structural and functional systems. Nat Rev Neurosci 10(3):186–198
Bureau E, Schilder F, Elmegråd M, Santos I, Thomsen J, Starke J (2014) Experimental bifurcation analysis of an impact oscillator-determining stability. J Sound Vib 333(21):5464–5474
Bureau E, Schilder F, Ferreira Santos I, Juel Thomsen J, Starke J (2013) Experimental bifurcation analysis of an impact oscillator—tuning a non-invasive control scheme. J Sound Vib 332(22):5883–5897
Buzsáki G (2004) Large-scale recording of neuronal ensembles. Nature Neurosci 7(5):446–451
Cagnan H, Duff E, Brown P (2015) The relative phases of basal ganglia activities dynamically shape effective connectivity in parkinson‘s disease. Brain 138(6):1667–1678
Calabresi P, Picconi B, Tozzi A, Ghiglieri V, Di Filippo M (2014) Direct and indirect pathways of basal ganglia: a critical reappraisal. Nature Neurosci 17(8):1022–1030
Chesselet MF, Delfs J (1996) Basal ganglia and movement disorders: an update. Trends Neurosci 19(10):417–422
Compte A, Brunel N, Goldman-Rakic P, Wang XJ (2000) Synaptic mechanisms and network dynamics underlying spatial working memory in a cortical network model. Cerebral Cortex 10(9):910–923
Cooper A, Stanford I (2000) Electrophysiological and morphological characteristics of three subtypes of rat globus pallidus neurone in vitro. J Physiol 527(2):291–304
Crowell A, Ryapolova-Webb E, Ostrem J, Galifianakis N, Shimamoto S, Lim D, Starr P (2012) Oscillations in sensorimotor cortex in movement disorders: an electrocorticography study. Brain 135(2):615–630
Dayal V, Limousin P, Foltynie T (2017) Subthalamic nucleus deep brain stimulation in parkinson‘s disease: the effect of varying stimulation parameters. J Parkinson‘s Dis 7(2):235–245
De Santos-Sierra D, Sendiña-Nadal I, Leyva I, Almendral J, Anava S, Ayali A, Papo D, Boccaletti S (2014) Emergence of small-world anatomical networks in self-organizing clustered neuronal cultures. PLoS ONE 9(1)
Deco G, Jirsa V, Robinson P, Breakspear M, Friston K (2008) The dynamic brain: from spiking neurons to neural masses and cortical fields. PLoS Computat Biol 4(8)
Deco G, Senden M, Jirsa V (2012) How anatomy shapes dynamics: a semi-analytical study of the brain at rest by a simple spin model. Front Comput Neurosci 6:1–7
Defazio G (2010) The epidemiology of primary dystonia: current evidence and perspectives. Eu J Neurol 17(SUPPL. 1):9–14
DeLong M (1990) Primate models of movement disorders of basal ganglia origin. Trends Neurosci 13(7):281–285
Deuschl G, Schade-Brittinger C, Krack P, Volkmann J, Schäfer H, Bötzel K, Daniels C, Deutschländer A, Dillmann U, Eisner W, Gruber D, Hamel W, Herzog J, Hilker R, Klebe S, Kloß M, Koy J, Krause M, Kupsch A, Lorenz D, Lorenzl S, Mehdorn H, Moringlane J, Oertel W, Pinsker M, Reichmann H, Reuß A, Schneider GH, Schnitzler A, Steude U, Sturm V, Timmermann L, Tronnier V, Trottenberg T, Wojtecki L, Wolf E, Poewe W, Voges J (2006) A randomized trial of deep-brain stimulation for Parkinson‘s disease. New Engl J Med 355(9):896–908
Dorval A, Kuncel A, Birdno M, Turner D, Grill W (2010) Deep brain stimulation alleviates parkinsonian bradykinesia by regularizing pallidal activity. J Neurophysiol 104(2):911–921
Dorval A, Russo G, Hashimoto T, Xu W, Grill W, Vitek J (2008) Deep brain stimulation reduces neuronal entropy in the mptp-primate model of parkinson‘s disease. J Neurophysiol 100(5):2807–2818
Droste F, Schwalger T, Lindner B (2013) Interplay of two signals in a neuron with heterogeneous synaptic short-term plasticity. Front Comput Neurosci 7:86
Ermentrout B, Terman D (2012) Neural networks as spatio-temporal pattern-forming systems. Springer, New York
Eusebio A, Brown P (2007) Oscillatory activity in the basal ganglia. Parkinson Related Disorders 13(SUPPL. 3):S434–S436
Fang J, Chen H, Cao Z, Jiang Y, Ma L, Ma H, Feng T (2017) Impaired brain network architecture in newly diagnosed parkinson‘s disease based on graph theoretical analysis. Neurosci Lett 657:151–158
Fiore V, Rigoli F, Stenner MP, Zaehle T, Hirth F, Heinze HJ, Dolan R (2016) Changing pattern in the basal ganglia: motor switching under reduced dopaminergic drive. Sci Rep 6
Gabor D (1946) Theory of communication. J IEEE London 93:429–457
Gafarov F (2016) Emergence of the small-world architecture in neural networks by activity dependent growth. Phys A Stat Mech Appl 461:409–418
Gagnon D, Petryszyn S, Sanchez M, Bories C, Beaulieu J, De Koninck Y, Parent A, Parent M (2017) Striatal neurons expressing d1 and d2 receptors are morphologically distinct and differently affected by dopamine denervation in mice. Sci Rep 7
Galvan A, Wichmann T (2008) Pathophysiology of parkinsonism. Clin Neurophysiol 119(7):1459–1474
Gouty-Colomer LA, Michel F, Baude A, Lopez-Pauchet C, Dufour A, Cossart R, Hammond C (2018) Mouse subthalamic nucleus neurons with local axon collaterals. J Compar Neurol 526(2):275–284
Guo Y, Rubin J, McIntyre C, Vitek J, Terman D (2008) Thalamocortical relay fidelity varies across subthalamic nucleus deep brain stimulation protocols in a data-driven computational model. J Neurophysiol 99(3):1477–1492
He BJ, Zempel JM, Snyder AZ, Raichle ME (2010) The temporal structures and functional significance of scale-free brain activity. Neuron 66(3):353–369
Holsheimer J, Dijkstra E, Demeulemeester H, Nuttin B (2000) Chronaxie calculated from current-duration and voltage-duration data. J Neurosci Methods 97(1):45–50
Holtzheimer P, Mayberg H (2011) Deep brain stimulation for psychiatric disorders. Annu Rev Neurosci 34:289–307
Huang Y, Hu K, Green A, Ma X, Gillies M, Wang S, Fitzgerald J, Pan Y, Martin S, Huang P, Zhan S, Li D, Tan H, Aziz T, Sun B (2020) Dynamic changes in rhythmic and arrhythmic neural signatures in the subthalamic nucleus induced by anaesthesia and tracheal intubation. Br J Anaesthesia 125(1):67–76
Izhikevich E, Kuramoto Y (2004) Weakly coupled oscillators. Encyclop Math Phys Five Volume Set, pp 448–453
Karl J, Ouyang B, Verhagen Metman L (2019) A novel dual-frequency deep brain stimulation paradigm for parkinson‘s disease. Neurol Therapy 8(2):483–489
Kita H, Kitai S (1991) Intracellular study of rat globus pallidus neurons: membrane properties and responses to neostriatal, subthalamic and nigral stimulation. Brain Res 564(2):296–305
Koeglsperger T, Palleis C, Hell F, Mehrkens J, Bötzel K (2019) Deep brain stimulation programming for movement disorders: Current concepts and evidence-based strategies. Front Neurol 10
Krack P, Volkmann J, Tinkhauser G, Deuschl G (2019) Deep brain stimulation in movement disorders: from experimental surgery to evidence-based therapy. Movement Disorders 34(12):1795–1810
Kuelbs D, Dunefsky J, Monga B, Moehlis J (2020) Analysis of neural clusters due to deep brain stimulation pulses. Biolo Cyber 114(6):589–607
Kuramoto Y (1984) Chemical oscillations, waves, and turbulence. Springer, New York
Kuramoto Y (1984) Cooperative dynamics of oscillator community—a study based on lattice of rings . Prog Theor Phys 79
Köhling R, Koch UR, Hamann M, Richter A (2004) Increased excitability in cortico-striatal synaptic pathway in a model of paroxysmal dystonia. Neurobiol Dis 16(1):236–245
Kühn A, Kempf F, Brücke C, Doyle L, Martinez-Torres I, Pogosyan A, Trottenberg T, Kupsch A, Schneider GH, Hariz M, Vandenberghe W, Nuttin B, Brown P (2008) High-frequency stimulation of the subthalamic nucleus suppresses oscillatory \(\beta \) activity in patients with parkinson‘s disease in parallel with improvement in motor performance. J Neurosci 28(24):6165–6173
Laing C, Chow C (2002) A spiking neuron model for binocular rivalry. J Comput Neurosci 12(1):39–53
de Lau L, Breteler M (2006) Epidemiology of parkinson‘s disease. Lancet Neurol 5(6):525–535
Luo F, Kiss Z (2016) Cholinergic mechanisms of high-frequency stimulation in entopeduncular nucleus. J Neurophysiol 115(1):60–67
Maiti P, Manna J, Dunbar G, Maiti P, Dunbar G (2017) Current understanding of the molecular mechanisms in parkinson’s disease: Targets for potential treatments. Transl Neurodegeneration 6(1)
Manos T, Zeitler M, Tass P (2018) Short-term dosage regimen for stimulation-induced long-lasting desynchronization. Front Physiol 9(APR)
Marschler C, Faust-Ellsässer C, Starke J., Van Hemmen J (2014) Bifurcation of learning and structure formation in neuronal maps. EPL 108(4)
McConnell G, So R, Hilliard J, Lopomo P, Grill W (2012) Effective deep brain stimulation suppresses low-frequency network oscillations in the basal ganglia by regularizing neural firing patterns. J Neurosci 32(45):15657–15668
Monga B, Wilson D, Matchen T, Moehlis J (2019) Phase reduction and phase-based optimal control for biological systems: a tutorial. Biol Cyber 113(1–2):11–46
Moon S, Cook K, Rajendran K, Kevrekidis I, Cisternas J, Laing C (2015) Coarse-grained clustering dynamics of heterogeneously coupled neurons. J Math Neurosci 5(1):1–20
Morrison A, Diesmann M, Gerstner W (2008) Phenomenological models of synaptic plasticity based on spike timing. Biol Cyber 98(6):459–478
Mottaghi S, Buchholz O, Hofmann U (2020) Systematic evaluation of dbs parameters in the hemi-parkinsonian rat model. Front Neurosci 14
Netoff T, Clewley R, Arno S, Keck T, White J (2004) Epilepsy in small-world networks. J Neurosci 24(37):8075–8083
Newman M (2003) The structure and function of complex networks. SIAM Rev 45(2):167–256
Ogawa N, Hirose Y, Ohara S, Ono T, Watanabe Y (1985) A simple quantitative bradykinesia test in mptp-treated mice. Res Commun Chem Pathol Pharmacol 50(3):435–441
Pavlides A, Hogan S, Bogacz R (2015) Computational models describing possible mechanisms for generation of excessive beta oscillations in parkinson‘s disease. PLoS Comput Biol 11(12):e1004609
Pavlides A, John Hogan S, Bogacz R (2012) Improved conditions for the generation of beta oscillations in the subthalamic nucleus-globus pallidus network. Eur J Neurosci 36(2):2229–2239
Pikovsky A, Rosenblum M, Kurths J (2001) Synchronization: a Universal Concept in Nonlinear Sciences. Cambridge University Press, Cambridge
Plenz D, Kital S (1999) A basal ganglia pacemaker formed by the subthalamic nucleus and external globus pallidus. Nature 400(6745):677–682
Popovych O, Tass P (2018) Multisite delayed feedback for electrical brain stimulation. Front Physiol 9(FEB)
Popovych O, Tass P (2019) Adaptive delivery of continuous and delayed feedback deep brain stimulation—a computational study. Sci Rep 9(1)
Prakash K, Bannur B, Chavan M, Saniya K, Kumar S, Rajagopalan A (2016) Neuroanatomical changes in parkinson‘s disease in relation to cognition: an update. J Adv Pharm Technol Res 7(4):123–126
Rubin J, Terman D (2004) High frequency stimulation of the subthalamic nucleus eliminates pathological thalamic rhythmicity in a computational model. J Comput Neurosci 16(3):211–235
Santaniello S, McCarthy M, Montgomery EB, J, Gale J, Kopell N, Sarma S (2015) Therapeutic mechanisms of high-frequency stimulation in parkinson’s disease and neural restoration via loop-based reinforcement. In: Proceedings of the National Academy of Sciences of the United States of America 112(6):E586–E595
Schilder F, Bureau E, Santos I, Thomsen J, Starke J (2015) Experimental bifurcation analysis—continuation for noise-contaminated zero problems. J Sound Vib 358:251–266
Schmidt H, Avitabile D, Montbrió E, Roxin A (2018) Network mechanisms underlying the role of oscillations in cognitive tasks. PLoS Comput Biol14(9)
Schwab B, Heida T, Zhao Y, Marani E, van Gils S, van Wezel R (2013) Synchrony in parkinson’s disease: Importance of intrinsic properties of the external globus pallidus. Front Syst Neurosci 7
Shefi O, Golding I, Segev R, Ben-Jacob E, Ayali A (2002) Morphological characterization of in vitro neuronal networks. Phys Review E Stat Phys Plasmas Fluids and Related Interdiscip Top 66(2)
Siettos C, Starke J (2016) Multiscale modeling of brain dynamics: from single neurons and networks to mathematical tools. Wiley Interdisc Rev Syst Biol Med 8(5):438–458
So R, Kent A, Grill W (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(3):499–519
So R, McConnell G, August A, Grill W (2012) Characterizing effects of subthalamic nucleus deep brain stimulation on methamphetamine-induced circling behavior in hemi-parkinsonian rats. IEEE Trans Neural Syst Rehabil Eng 20(5):626–635
So R, McConnell G, Grill W (2017) Frequency-dependent, transient effects of subthalamic nucleus deep brain stimulation on methamphetamine-induced circling and neuronal activity in the hemiparkinsonian rat. Behav Brain Res 320:119–127
Spiliotis K, Siettos C (2011) A timestepper-based approach for the coarse-grained analysis of microscopic neuronal simulators on networks: Bifurcation and rare-events micro- to macro-computations. Neurocomputing 74(17):3576–3589
Stam C, Reijneveld J (2007) Graph theoretical analysis of complex networks in the brain. Nonlinear Biomed Phys 1
Stephens B, Mueller A, Shering A, Hood S, Taggart P, Arbuthnott G, Bell J, Kilford L, Kingsbury A, Daniel S, Ingham C (2005) Evidence of a breakdown of corticostriatal connections in parkinson‘s disease. Neuroscience 132(3):741–754
Strogatz S (2001) Exploring complex networks. Nature 419(1):268–276
Tass P (1999) Phase resetting in medicine and biology—stochastic modelling and data analysis. r. Springer, Berlin
Terman D, Rubin J, Yew A, Wilson C (2002) Activity patterns in a model for the subthalamopallidal network of the basal ganglia. J Neurosci 22(7):2963–2976
Udupa K, Chen R (2015) The mechanisms of action of deep brain stimulation and ideas for the future development. Prog Neurobiol 133:27–49
Vidailhet M, Jutras MF, Grabli D, Roze E (2013) Deep brain stimulation for dystonia. J Neurol Neurosurg Psychiatry 84(9):1029–1042
Wang D, de Hemptinne C, Miocinovic S, Ostrem J, Galifianakis N, Luciano M, Starr P (2018) Pallidal deep-brain stimulation disrupts pallidal beta oscillations and coherence with primary motor cortex in parkinson‘s disease. J Neurosci 38(19):4556–4568
Watts D, Strogatz S (1998) Collective dynamics of “small-world” networks. Nature 393(6684):440–442
West T, Berthouze L, Halliday D, Litvak V, Sharott A, Magill P, Farmer S (2018) Propagation of beta/gamma rhythms in the cortico-basal ganglia circuits of the parkinsonian rat. J Neurophysiol 119(5):1608–1628
Wichmann T, Dostrovsky J (2011) Pathological basal ganglia activity in movement disorders. Neuroscience 198:232–244
Wojtecki L, Timmermann L, Jörgens S, Südmeyer M, Maarouf M, Treuer H, Gross J, Lehrke R, Koulousakis A, Voges J, Sturm V, Schnitzler A (2006) Frequency-dependent reciprocal modulation of verbal fluency and motor functions in subthalamic deep brain stimulation. Arch Neurol 63(9):1273-1276
Xu W, Russo G, Hashimoto T, Zhang J, Vitek J (2008) Subthalamic nucleus stimulation modulates thalamic neuronal activity. J Neurosci 28(46):11916–11924
Acknowledgements
KS would like to acknowledge Georgios Georgiou and Constantinos Siettos for their help, discussions, and encouragement during the early stages of this research. All authors thank Susann Dittmer for her help drawing Fig. 1.
Open Access
This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
Funding
Open Access funding enabled and organized by Projekt DEAL. The authors thank the DFG for support through the Collaborative Research Center CRC 1270 (Deutsche Forschungsgemeinschaft, Grant/ Award Number: SFB 1270/1–299150580).
Author information
Authors and Affiliations
Corresponding author
Additional information
Communicated by Benjamin Lindner.
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Spiliotis, K., Starke, J., Franz, D. et al. Deep brain stimulation for movement disorder treatment: exploring frequency-dependent efficacy in a computational network model. Biol Cybern 116, 93–116 (2022). https://doi.org/10.1007/s00422-021-00909-2
Received:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s00422-021-00909-2