Effect of Stimulation of Subthalamic Nucleus onbeta Oscillations and Thalamus Tremor Activity in aComputational Model of Parkinson’s Disease  


 Parkinson’s disease (PD) is associated with abnormal b band oscillations (13-30 Hz) in the cortico-basal ganglia circuits.Abnormally increased striato-pallidal inhibition and strengthening the synaptic coupling between subthalamic nucleus (STN)and globus pallidus externa (GPe), due to the loss of dopamine, are accounted as the potential sources of b oscillations in thebasal ganglia. Deep brain stimulation (DBS) of the basal ganglia subregions is known as a way to reduce the pathological boscillations and motor deficits related to PD. Despite the success of the DBS, its underlying mechanism is poorly understoodand, there is controversy about the inhibitory or excitatory role of the DBS in the literature. Here, we utilized a computationalnetwork model of basal ganglia which consists STN, GPe, globus pallidus interna (GPi), and thalamus neuronal population.This model can capture healthy and pathological b oscillations as what has been observed in experimental studies. Using thismodel, we investigated the effect of DBS to understand whether its effect is excitatory or inhibitory. Our results show that theexcitatory DBS (EDBS) is able to quench the pathological synchrony and b oscillations, while, applying inhibitory DBS (IDBS)failed to quench the PD signs. In addition, the EDBS ameliorated the thalamic activity related to tremor in the model, while,the IDBS outperformed. However, with the help of the model results, we conclude that the effect of the DBS on its target isexcitatory

Reduction of firing rate of the stimulated neuronal area has been observed in human 38,39 and monkeys 40 with PD which remarks the inhibitory role of the DBS. Moreover, several mechanisms suggested to explain the inhibitory role of the DBS such as depolarization block 36,41 , inactivation of voltage-gated currents [42][43][44] , and activation of inhibitory afferents 38,40,[45][46][47][48][49][50] . In a study has been suggested the effect of DBS is inhibitory 45 . its authors showed the GPi neurons were inhibited during stimulation, while the neurons were inhibited using GABA blocker. On the other hand, the excitatory role of the DBS has also been suggested by several studies 48,51 . Applying DBS on internal segment of globus pallidus (GPi) reduces firing rates of thalamic neurons which are inhibited by the GPi 52 . Also, applying the DBS on STN neurons (the excitatory neuronal population) increases firing rate of GPi, globus pallidus externa (GPe), and substantia nigra pars reticulata (SNr) of human and animal with PD [53][54][55] which support the excitatory effect of the DBS.
where I L , I Na , I K , I Ca , I T , and I AHP are the leak, sodium, potassium, high threshold calcium, low threshold calcium, and after hyper polarization currents, respectively. I app is the external current applied to the neurons (i.e. the DBS current). I pre−>post is synaptic current from the presynaptic to the postsynaptic neuron. X represents gating channels such as potassium channels (n), opening (m) and closing (h) sodium channels, and low threshold calcium channels (r). The τ X (v) in equation 2 is defined as follows: While, in the GPe and GPi neurons the τ X (v) is constant and equal to τ r . The ionic currents used in equation 1 were computed as follows: In the equations 4 to 8 the X = n, h is the ionic gating channel variables (h for closing sodium channel and n for potassium). In these equations the X ∞ = m, a, r or s is the steady state of the ionic gating channels (m for opening sodium channel, a for T-type and s for L-type calcium channel) and is computed by the equation 9.

3/23
But, the function b ∞ (r) used in 8 is computed with different equation: The after hyper-polarization (AHP) current used in equation 1 (I AHP ) is where the [Ca] is the intra-cellular calcium concentration: The parameters and their values of STN, GPe, and GPi neurons are presented in table 1 to 3.
Membrane potential of thalamus neurons in the network model is computed using the following differential equations: where I L , I Na , I K , and I T are the leak, sodium, potassium, and low threshold calcium currents, respectively. I GPi−>T h is the synaptic current from a GPi neuron to a thalamus neuron in the network model. The I SMC represents cortico-thalamic sensorimotor pulses applied to the thalamus neurons. The equations 16 to 19 compute the ionic currents used in equation 13.
and the functions used in equations 14 to 19 are computed as follows: The parameters and their values of thalamus neurons are presented in table 4. In addition to modification of the network structure we modified the network model parameters compared to the Terman et al. (2004) as follows: for the STN neurons g Na was decreased from 37.5 to 30 nS/µm 2 . The g K was decreased from 45 to 40 nS/µm 2 . The value of φ was taken from Park et al. (2011; i.e. φ n = φ h = 5, φ r = 2). The value of ε was considered to be 3 × 10 −5 ms −1 . The g GPe−>ST N and I app of the STN in the healthy state is 2.2 nS/µm 2 and 8.4 pA/µm 2 , respectively. In the PD state of the network model these two parameters were changed to 7 nS/µm 266-69 and 3 pA/µm 270 , respectively.  To simulate GPe neurons in the network model, we set φ h = 0.135, φ n = 0.165, φ r = 1, and ε = 0.0055 (similar to Park et al., 2011). The g GPe−>GPe , g ST N−>GPe , and I app of GPe in the healthy state were 0.01, 0.01 nS/µm 2 , and 5.9 pA/µm 2 , respectively. To simulate the PD state of the network model, these parameters were changed to 0.9, 0.55 nS/µm 266-69 , and 0.5 pA/µm 2 , respectively. Note that in the PD state of the network model I app of the GPe and GPi decreases leading to less activity of the GPe

6/23
and GPi neurons due to the increasing striatal inhibition (explained in 56 ) in the network model. Parameters of the GPi neurons are similar to the GPe neurons with the difference that for the GPi neurons, φ h = 0.1 and φ n = 0.135. The g GPe−>GPi , g ST N−>GPi , and I app of GPi in the healthy state were 0.01, 0.005 nS/µm 2 , and 7.7 pA/µm 2 . To simulate the PD state of the network model, these parameters similar to GPe neurons were changed to 1.9, 1.1 nS/µm 2 , and 4 pA/µm 2 , respectively. Parameters of the thalamic neurons are the same as in 18 with the difference that in our network model g GPi−>T h was 0.05 nS/µm 2 . These modifications moved our network model activity more close to the experimental results.
The synaptic model used here is of conductance based type similar to the model used in 17,18,61,62,71,72 . The synaptic currents used in equations 1 and 13 are computed as follows: where the j is the index of presynaptic neuron. The parameter s in equation 28 is where the H ∞ (v) as follows: population firing rate To compute the time resolved population firing rate for each neuronal population in the network model we used 10 milliseconds sliding window and shifted with steps of 1 millisecond over the entire simulation time while for each step we counted the number of spikes for all neurons in the population and converted it to spikes per second (i.e. to Hz).

Sensorimotor and DBS pulses
The sensorimotor and DBS pulses are simulated using the following equation: where A, f , t, and δ are pulse amplitude, frequency, time (in ms), and pulse duration (in ms), respectively. The H(.) is the Heaviside function. To find an appropriate value for the A for the excitatory and inhibitory DBS, we computed the STN mean firing rate (over 20 trials) by varying A = 50 pA/µm 2 to 300 pA/µm 2 ( Figure 1-C), then we chose the constraint the firing rate of the STN in excitatory DBS (EDBS) and inhibitory DBS (IDBS) to be equal in the cases. Therefore, to simulate EDBS, we set A = 126.57 pA/µm 2 and 147.36 pA/µm 2 , f = 150 Hz, and δ = 0.1 ms. Parameter settings for IDBS were the same as EDBS, with the difference that A = −126.57 pA/µm 2 and −147.36 pA/µm 2 . The DBS parameters are adapted from 10,54,[73][74][75][76][77] . For sensorimotor pulses (i.e. cortico-thalamic input), A = 4.5 pA/µm 2 , f = 20 Hz, and δ = 5 ms. To generate irregular pulses we used Poisson process in equation 31.

Thalamus fidelity
Thalamus neurons in the network model show four types of responses to the cortico-thalamic input pulse: 1) Correct spike: a single thalamic spike in response to a cortico-thalamic input pulse. 2) Missed spike: refers to the case when there is no thalamic spike in response to a cortico-thalamic input spike. 3) Extra spike: refers to the case when a thalamic neuron shows more than one spike in response to a cortico-thalamic input pulse. 4) Undesired spike: occurs when a thalamic neuron spikes while there is no cortico-thalamic input spike. According to the four response types of the thalamus neurons to a single cortical pulse, the thalamus fidelity is computed using the following equation: where the N M , N E , and N U are the number of missed spikes, extra spikes, and undesired spikes, respectively. The N exp is the number of expected thalamic spikes due to cortico-thalamic pulses. Since each cortical pulse is given to all thalamus neurons in the network model, we expect to observe that each thalamus neuron, relaying the cortico-thalamic pulse, emits a spike in response to the cortico-thalamic pulse. Therefore, the number of expected thalamic spikes in response to the cortical inputs (i.e. N exp ) equals to the number of thalamus neurons multiplied by the number of cortico-thalamic pulses 18, 78 .

Synchrony index
We used Fano Factor (FF) to measure synchronous spiking activity for each neuronal population in the network model. To compute FF we used the following equation: where the Var(.) and E(.) are the variance and mean of the population firing rate (PFR), respectively. Higher FF values represent more synchrony in the spiking activity of a neuronal population in the network model 56,79 .

Mean power spectral density
The power spectral density (PSD) of the population firing rates was computed using the Welch's method in python 2.7 (i.e. using scipy.signal.welch python package; 80 ). The sampling rate and the segment length were set to 1000 Hz and 1000 data points, respectively. Other parameters required for the scipy.signal.welch function were set to the predefined default values (see https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.signal. welch.html). The mean power spectral density was computed by averaging over 50 simulations.

Oscillation index
Oscillation index was computed by dividing the area under the curve of a PSD in the β frequency range (i.e. between 13 to 30 Hz) by the area under the curve for the whole frequency range (i.e. from 1 to 500 Hz).

Tremor frequency
We assumed that the extra spikes of the thalamus neurons during the resting state of the network model (i.e. when there is no cortico-thalamic sensorimotor pulses) are related to limb tremor. So, to compute the tremor frequency, we measured the mean firing rate of thalamus neurons (in Hz), over the whole simulation period, during the resting state.

Simulation
The simulations were implemented in python 2.7. All differential equations were solved using odeint from SciPy library 80 with 0.05 ms time resolution (see https://docs.scipy.org/doc/scipy-0.18.1/reference/generated/ scipy.integrate.odeint.html). To reduce the simulation time, we performed parallel programming using python message passing interface (MPI) in a cluster computing with 30 core processors (Intel 3.2 GHz). To avoid the initial transient network model responses, we did not consider the first 250 ms of each simulation in our analysis.

The network model captures features of the healthy and PD BG
The activity of different BG regions in the healthy state is non-oscillatory and desynchronized 3, 81-83 . This feature is captured in our network model. Similar to the experimental results, the STN spiking activity in the healthy network model is asynchronous irregular (Figure 2A, top panel; the same for the GPe and GPi spiking activity, data not shown). This is also reflected in the STN population activity in the healthy state of the network model ( Figure 2B, top panel). The STN population activity in the healthy network model is non-oscillatory. This leads to a flat PSD of the STN population firing rate in the healthy state ( Figure 2C). Altogether, these results indicate that the healthy activity of the network model is irregular, inline with the experimental studies showing that the BG activity in the healthy state is non-oscillatory. The STN, GPe, and GPi mean population firing rates, in the healthy state of the network model, are 19.4 ± 1.1 Hz, 45.47 ± 1.2 Hz and 56.52 ± 2 Hz, respectively which match the previously reported experimental values 83,84 .
In the PD state, the network model neurons show synchronized bursts of spiking activities in the β frequency range (see the middle panels in Figures 2A and 2B). This also matches experimental studies which indicate synchronized β band oscillatory spiking activities in the BG as a hallmark of PD 7-10, 84-98 .
To bring the network model from the healthy to the PD state, we followed three steps. First, the I app applied to the GPe and to the GPi neurons (Equation 1) was decreased from 5.9 pA/µm 2 and 7.7 pA/µm 2 (healthy state) to 0.5 pA/µm 2 and 4 pA/µm 2 (PD state), respectively. This leads to reduction in the activity of the GPe neurons (39.1 ± 0.8 Hz; independent two-tailed t-test, p < 0.001) in the PD state of the network model, compared to the healthy state ( Figure 2D). This is inline with the experimental studies indicating that the GPe firing activity decreases during PD 83,99,100 . However, despite decreasing the I app applied to the GPi neurons, the GPi firing rate increases (64.9 ± 0.86 Hz; independent two-tailed t-test, p < 0.001) compared to the healthy state ( Figure 2D; 99,101,102 ). The reason is that the lower GPe activity during the PD state of the network 8/23 model disinhibits the GPi neurons. Thereby, the GPi population firing rate in the PD network model increases compared to the healthy state ( Figure 2D).
Second, the I app applied to the STN neurons, representing cortico-subthalamic input, in the network model (Equation 1) was decreased from 8.4 pA/µm 2 (healthy state) to 3 pA/µm 2 (PD state). Such a change in the network model is inline with the experimental studies showing that the cortical activity decreases in PD 70 which can lead to less cortico-subthalamic drive, due to direct cortico-subthalamic connectivity [103][104][105] . Note that despite decreasing the I app of the STN neurons in the PD state of the network model, the STN activity increases (27.9 ± 1.5 Hz; independent two-tailed t-test, p < 0.001) compared to the healthy state ( Figure 2D; 84,85 ). The reason for this is STN disinhibition due to reduction in the activity of the GPe units in the PD state of the network model ( Figure 2D Third, the synaptic connectivity in the subthalamo-pallidal circuit (STN to GPe and GPe to STN synapses) were strengthened in the PD state, compared to the healthy state (see materials and methods). Such subthalamo-pallidal synaptic strengthening is inline with the experimental studies showing that both STN to GPe and GPe to STN synapses are strengthened during PD [66][67][68][69] .
Applying these three changes brings the network model from the healthy state to the PD state. STN neurons of the PD network model show synchronized bursts of spiking activity in the β frequency range (Figure 2A; the GPe and GPi neurons show the same behavior). The STN PD-like β band oscillations are also observed in the STN population firing rate (Figures 2B, 2C) as well as in the GPe and GPi population activities (Figures 2E, 2F). In addition, we measured the burst rate of GPi neurons and we found a significant increase in PD state compared to the healthy (p < 0.001, independent two-tailed t-test; Figure 4A) which matches with the previous experimental 54 and computational studies 61 .
To check whether the thalamus neurons are able to process the cortical inputs (i.e. cortico-thalamic motor commands), we measured thalamus fidelity in the healthy and PD states of the network model. Our results indicate that thalamus fidelity (see materials and methods) to the cortical motor commands decreases in the PD state of the network model compared to the 9/23 healthy state ( Figure 2G). This result is inline with the previous computational studies 18,59,60,78,106 indicating reduction of the thalamus fidelity during PD.
All in all, our network model can reproduce features of the experimental data for both healthy and PD state of the BG. Mainly, the network model show asynchronous irregular spiking activity in the healthy state and synchronous β band oscillatory activity in the PD state. In addition, thalamus fidelity decreases in the PD state of the network model compared to the healthy state.

Effects of EDBS and IDBS on the network model
The STN high frequency DBS has therapeutic effects on PD signs such as reducing the pathological β oscillations in the cortico-BG loop 33,36,74,[107][108][109] , and improving PD-related motor symptoms [27][28][29][30][31][32][110][111][112] . However, the mechanism(s) underlying STN DBS is yet unknown. To understand whether the STN DBS therapeutic effects are due to excitation of STN or inhibition of STN, we applied EDBS (i.e. excitatory DBS) and IDBS (i.e. inhibitory DBS) to the STN neurons in the network model and investigated the effect(s) of each DBS type. We tested which DBS type can quench the PD-like β oscillations in the network model. Furthermore, we also tested which DBS type can improve the thalamus fidelity to the cortical motor commands. Here, we used 126.57pA and 147.36pA current amplitudes for both types of DBS. The STN mean firing rate is approximately the same in these two amplitudes (Figure 2 and 3D; not significant, independent two-tailed t-test; see materials and methods). Due to the similarity of results we reported the DBSs with 147.36pA current except for the Figure 5.

EDBS effects
Applying the EDBS (see materials and methods) to the STN neurons in the network model quenches the PD-like β oscillations in the STN (and in the other regions included in the network model; Figures 2A-C, bottom panels). We applied two EDBS types to the PD network model (i.e. regular and irregular EDBS; regular DBS is rectangular pulses with constant phase and the irregular DBS is similar to the regular with random phases; see materials and methods for more details on each EDBS type). To quantify the performance of regular and irregular EDBS on the PD-like β oscillations in the network model, we measured the Fano factor and the oscillation index of the STN, GPe, and GPi population firing rates while the STN was stimulated. For regular EDBS, both Fano factor and oscillation index dramatically decreased compared to the PD state (P < 0.001, independent two-tailed t-test), for all regions measured (i.e. STN, GPe, and GPi; Figures 2E, F). We also compared the performance of the regular and irregular EDBS on improving the thalamus fidelity in the network model. We found that both regular and irregular EDBS increased the thalamus fidelity (almost to the healthy level; Figure 2G). In addition, we found that the burst rate of the GPi neurons significantly decreased (P < 0.001, independent two-tailed t-test) in both types of EDBS state compared to the PD and also to the healthy ( Figure 4A) which matches the previous experimental 54 and computational studies 61 .
All in all, this result indicates that applying regular/irregular EDBS to the STN neurons in the network model, quenches the PD-like β oscillations and improves the thalamus fidelity to the cortical motor command.

IDBS effects
Thus far, we showed that whether regular or irregular high frequency EDBS is able to quench PD-like β oscillations and improve the thalamus fidelity in the network model. Next, we tested whether high frequency inhibition of the STN neurons (IDBS) has the same effects as the EDBS in the network model. To do that, we stimulated the STN neurons in the pathological network model by high frequency inhibitory pulses (see materials and methods) and measured the PSD, Fano factor, and oscillation index of the STN population firing rate as well as the thalamus fidelity in the network model ( Figure 3). Our simulation results show that regular and irregular IDBS (unlike the EDBS) fail to quench PD-like β oscillations and to improve the thalamus fidelity ( Figure 3). Although the burst rate of the GPi neurons in regular and irregular IDBS slightly decreases (p < 0.001, independent two-tailed t-test) compared to the PD state, the GPi burst rate were still above the healthy state when the network model exposed to regular and irregular IDBS ( Figure 4A). However, comparing the performance of EDBS and IDBS in the network model reveals that EDBS outperforms IDBS in improving PD signs (i.e. quenching PD-like β oscillations and improving the thalamus fidelity) in the network model. In the following, we explain why IDBS fails to improve the PD signs in the network model.

IDBS can not treat the pathological STN rebound bursting activity in the network model
We showed that only EDBS (and not IDBS) is able to quench the PD-like β oscillations in the network model (Compare Figure 2 with Figure 3). The reason why IDBS fails to quench the PD-like β oscillations in the network model is rebound bursting activity of the STN neurons due to IDBS. In line with the experimental studies 87,113,114 , the STN neurons in the network model show rebound bursts of spiking activity when the strong inhibition finished. As the IDBS is a barrage of inhibitory inputs to the STN neurons in the network model (see materials and methods), the STN neurons react to it by rebound bursts of spiking activities ( Figure 5C). STN rebound bursts lead to increase in the GPe spiking activity through subthalamo-pallidal excitatory pathway. Then, the higher GPe activity gives rise to inhibition of the STN neurons which resulted in rebound bursting due the T-type calcium current ( Figure 5C). Such a mechanism retains the PD-like β oscillations in the network model during IDBS. Therefore, the reason why IDBS fails to quench the PD-like β oscillations in the network model is propagation of the STN rebound bursts through the subthalamo-pallidal loop.
During EDBS, two behaviors can occur depending on the current amplitude. First, for 147.36pA, the STN neurons in the network model do not show rebound bursting activity (Figures 2, and 5A). The reason is that the STN spiking activity is driven by the EDBS pulses, counteracts with the pallido-subthalamic inhibition, thereby no STN rebound bursting activity can happen due to the non-sufficient T-type calcium current. Second, for 126.57pA, the EDBS extends the burst durations and the inter-burst intervals which resulting in quenching the PD-like β oscillations ( Figure 5B).

Tremor signs in healthy and PD state of the network model
So far, using our simulation results, we showed that the STN EDBS (and not IDBS) can quench pathological β oscillations emerging in the PD state of the network model. Next, we investigated the effect(s) of each stimulation scenario (i.e. EDBS or IDBS) on the resting state activity of the healthy and the PD network model. To simulate the resting state of the network model, we excluded the cortico-thalamic sensorimotor pulses (i.e. by setting the I SMC in the equation 13 to zero for both healthy and PD network model simulations; see materials and methods). Our simulation results indicate that the thalamus neurons in the healthy network model, do not show any spiking activity during the resting state ( Figure 6A). However, in the PD network model, the thalamus neurons show the mean spiking activity of 8.24 ± 1 Hz during the resting state ( Figure 7C). We considered this as resting state tremor-like activity in the network model because of three reasons. First, such thalamus activity is specific to the PD state of the network model (compare Figures 6A and 7C). Second, according to the experimental studies [115][116][117] the tremor is driven by the pathological thalamus activity in PD. Third, the resting state pathological thalamus activity in the network model approximately is 8 Hz which matches the tremor frequency reported in the experimental studies 63, 64, 88, 118-121 . In the following, we explain why such tremor-like thalamus spiking activity emerges only in the PD network model (and not in  Figure 1A; note that the cortico-thalamic input is set to zero during the resting state). Thereby, the resting state thalamus activity in the network model is driven by inhibitory inputs from the GPi. Our resting state simulation results show that the healthy GPi spiking activity is tonic ( Figure 6B, bottom panel). Such tonic GPi spiking activity or weak bursting activity (burst with lower number of spikes; Figure 6B, two bottom panels; and also see Figure4B and C) is not strong enough to sufficiently flow the T-type calcium current of the thalamus neurons (see the T-Current in Figure 5B). This leads to no thalamus spiking activity during the resting state and in the healthy condition of the network model ( Figure 5B, top panel).

STN EDBS improves the thalamus tremor-like activity in the network model
Next, we investigated the effect(s) of EDBS and IDBS on the tremor-like thalamus activity in the network model. Mainly, we tested which DBS type (i.e. EDBS or IDBS) can reduce the resting state tremor-like thalamus spiking activity in the PD network model. Our simulation results show that applying EDBS to the STN neurons in the network model leads to a better

13/23
performance as compared to applying the IDBS (Figure 6). In other words, STN EDBS reduces the tremor-like thalamus spiking activity nearly to the healthy state while STN IDBS fails to do so ( Figure 7A-C ). In the following, we explain why STN EDBS outperforms STN IDBS in improving the tremor-like thalamus activity in the network model.
In line with previous experimental and computational studies 54,61 , GPi neurons in the PD state of our network model switch between long lasting hyperactivity and no activity states when the STN is stimulated by the EDBS (Figure 7D, bottom panel and also see the burst duration and number of spikes per burst in Figure 4B and C). During the GPi hyperactivity periods, thalamus neurons are strongly inhibited and thereby do not show tremor-like spiking activity anymore ( Figure 7D). However, a single thalamus spike occurs at the end of the GPi activity ( Figure 7D, top panel; also see Figure 7A). This single thalamus spike is due to after-hyperpolarization increase in the T-type calcium current of the thalamus neuron receiving the GPi strong inhibitory input ( Figure 7D).
However, when the STN is exposed to the IDBS, the GPi neurons in the network model show regular oscillatory bursts of spiking activity in ( Figure 7E, bottom panel). These strong GPi bursts (similar property to PD state; see Figure 4B and C) cause sufficient flow of the T-type calcium current of the thalamus neurons in the network model (see the T current fluctuations in Figure 7E) which consequently, leads to tremor-like thalamus spiking activity.
Altogether, our simulation results indicate that the EDBS outperforms the IDBS, not only in quenching the PD-like β oscillations in the network model (Figure 2), but also in reducing the tremor-like spiking activity of the thalamus neurons nearly to the healthy state ( Figure 7C). Tremor-like frequency of thalamus in the healthy and PD states and when the STN is exposed to regular/irregular EDBS and regular/irregular IDBS (error bars show standard deviation). (D) From top to bottom, membrane potential, low threshold T-type calcium current, synaptic inhibitory input from the connected GPi and membrane potential of the corresponding GPi neuron (see Materials and methods) when STN is exposed to EDBS in the network model. (E) The same as D, when STN is exposed to the IDBS.

14/23
Discussion Permanent and excessive β oscillations in BG is the hallmark of the PD [11][12][13] . The DBS quenches β oscillations and improves motor symptoms related to PD 27-33, 36, 74, 107-112 . But, the mechanism of the DBS is poorly understood. In the present study, we investigated the generation of the pathological β oscillations in the BG by neural modeling and the effects of the DBS in different scenarios. We showed that the subthalamo-pallidal circuit is the potential source of the generation of the β oscillations based on the computational network model. Our findings confirm the rebound burst activity of the STN neurons is the key reason to generation of the β oscillations as what has been suggested in 17,23 . Then we investigated the role of the DBS on quenching the β oscillations. We applied inhibitory and excitatory DBS on STN neurons to compare their effect with those observed in experimental studies. Thereby, we tuned current amplitude of each DBS types to have equal STN firing rates. The results show that the EDBS can counteract the pallido-subthalamic inhibition and therefore, the rebound burst of the STN neurons were controlled ( Figure 5). Also, EDBS can extend the duration of STN bursts which results in quenching the β oscillations ( Figure 5). While, the IDBS is not able to quench the β oscillations due to its ineffectiveness in the suppression of rebound burst activity of STN neurons and in the extension of burst duration. Our results suggest that effect of DBS on STN neurons is excitatory. Our simulation results also help to clarify the relationship between tremor and the activity of the basal ganglia and thalamus. We showed that the postinhibitory rebound activity of the thalamus neurons due to the strong inhibition of the GPi neurons is the potential reason of the resting state tremor. Like other signs of PD, the tremor activity of thalamus was quenched by EDBS while the IDBS was worsened it.
In the following, we discuss about the network model, the generation of β oscillations, the tremor, and the role of DBS in improving the PD signs in detail.

Network model
Our BG model has been created using Hodgkin-Huxley type neuron model that can generate neural behavior in detail (i.e. ion currents). Previous computational BG studies 17, 18, 57-62, 72, 122 except 24, 123 which were created by Hodgkin-Huxley type neurons, were not reported β oscillations in PD state. In these studies, one of changes to set the model in PD state is strengthening subthalamo-pallidal synapses, based on the experimental results [66][67][68][69] . The PD β oscillations also was reported by the model proposed in 56 which was created using LIF neuron model. This model moves to PD state only by increasing striato-pallidal spiking activity (without changing any synaptic strength). They show the PD β oscillations might arise from the network effects, without taking into account the details of neuronal dynamics. While, 23 used integrate-and-fire-or-burst neuron in their network model for subthalamo-pallidal circuit. Their model moves to PD state by changing in synaptic strength and input spikes in subthalamo-pallidal circuit and represents the PD β oscillations. In this model, the rebound burst activity of the STN and GPe neurons played an important role in the generation of the β oscillations. The PD β oscillations has also reported in the firing rate models 20,[124][125][126] . These models cannot be used for single neuron study in the network. By contrast, the role of the calcium bursting can not be demonstrated in these models. However, the modified BG model in present study was be able to generate the β oscillations in PD state. Also, the other properties of BG model (i.e. firing rates and changes in firing rates when moving from the healthy to the PD state) were matched the experimental observations. these modifications were made by changing in some intrinsic and connectivity neuronal properties (see materials and methods).

Generation of beta oscillations
Computational and experimental studies have implicated the subthalamo-pallidal circuit is the potential source of the β oscillations in PD 7,8,[14][15][16][17][18][19][20][21][22][23][24]78 . Moreover, the induction of β oscillations from cortex into the BG has been claimed by 6,25,26 . However, more detailed models for the generation of the β oscillations, the rebound burst activity of the STN neurons due to the inhibition of the GPe neurons in the BG has been considered 15,16,114 . Somehow, mentioned hypothesize is confirmed when the motor symptom of patients with PD suppressed after receiving T-type calcium channel blocker such as Zonisamide 127 . Also, the computational studies such as 17,18,23,24,123 confirmed this hypothesis. While, the other studies demonstrate that the excessive inhibition on inhibitory population (such as striato-pallidal inhibition) and/or excessive excitation on excitatory population (such as cortico-subthalamic) result in β oscillations 20,56,72,128,129 . Our findings confirm the role of post-inhibitory rebound bursts of STN neurons in generation of the β oscillations.

The role of the DBS is excitatory
The DBS improves PD-related motor symptoms [27][28][29][30][31][32][110][111][112] and PD-related neuronal behavior in BG 10-13, 33, 36, 74, 107-109 . In our network model, we quantified the PD-related signs by β oscillations, thalamus fidelity, and tremor frequency (see below and materials and methods). Here, we applied two scenarios of the DBS in our network model to see which satisfy our expectation about the amelioration of PD signs with DBS. Our findings suggest that the EDBS improves the PD neuronal behavior and the motor symptoms while the IDBS worsened them. Although the previous studies 38-41, 43-47, 49, 50, 130 indicate that the role of the DBS on its target is inhibitory, in 48,51,[53][54][55] has been demonstrated the opposite role in accordance with our 15/23 findings. On the other hand, the computational study in 56 has claimed the excitation of excitatory and/or inhibition of inhibitory populations leads the network to oscillations (in this case is β oscillations). Therefore, the pathological β oscillations have been quenched using high frequency inhibitory spike trains on STN, that is not consistent with EDBS in our network. still, there is opposite evidence for this study's results: the initiation of movement accompanied with increasing cortical activity [131][132][133][134] and the STN neurons are excited by cortex [103][104][105] which quenches briefly the β oscillations in patients with PD 8,135,136 and in a computational model proposed in 24 . Moreover, the neuronal bursting activity were not investigated in that study, while, in 137 a computational model was proposed that the neuronal bursting activities have been investigated. This study showed the variety of STN burst range can affect on β oscillations. In detail, the low burst rate of STN neurons quenched the β oscillations while the high STN burst rate generated the β oscillations with a little shift in frequency. Indeed, the intrinsic behavior of STN neurons in non-pathological state are bursty and the higher burst rate of BG nuclei in PD condition are also shown ( 138 ; and see also Figure 4 for GPi) which is consistent with the results that reported in 137 . The consistency of this study and our network model is justified with considering that 1-the STN burst rate in our network model corresponds to high STN burst rate of the model in 137 and 2-with considering the STN burst rate when it exposes to EDBS corresponds to low STN burst rate in that study. Consequently, it seems to the therapeutic effects of DBS acts by excitation of the STN.
In addition, in 37,139,140 it is shown that the DBS inhibits the STN neurons while it induces spike in their axons. In this hypothesis, despite the inhibition of STN neurons, their post-synaptic neurons (GPe and GPi) receive glutamate. So, inducing excitatory pulses in STN neurons in the network model matches with this hypothesis.
In addition, our results suggested that the post inhibitory rebound burst of the STN neurons is the main cause of the β oscillations. The EDBS by counteracting the pallidal inhibition in the network model eliminates the cause of the rebound bursting of the STN neurons. While, the IDBS helps the pallidal inhibition that resulted in more rebound burst activity of the STN neurons in the network model. Altogether, based on our network model results we conclude that the role of the DBS on its target is excitatory (and not inhibitory).

Tremor
The tremor in patients with PD occurs when they are at rest. The frequency of the resting tremor in PD state is reported 3 to 8 Hz 64, 88, 118-121 and it correlates with thalamus neuronal activity [115][116][117] . Previous computational studies 17,18,57 demonstrated tremor activity by single neuron spiking activity in the tremor frequency in BG populations such as STN. While, in 88 the correlation between high frequency activity of STN and tremor has been shown. In a computational study the tremor activity in the PD state has been represented by synaptic input from the GPi to the thalamus neurons in tremor frequency 62 . However, the resting state was not simulated in mentioned computational studies. In our network model, we simulated resting state with interrupting sensorimotor command pulses to the thalamus neurons. Our results, show the extra thalamic spikes in tremor frequency in PD state (and not in healthy state). This result was consistent with previous experimental tremor-related thalamic activity observations. The main reason of extra spiking activity of thalamus in PD state is post inhibitory rebound spiking activity due to the strong inhibitory input from the GPi. The T-type calcium current of the thalamus neurons plays an important role in the rebound spiking activity.
However, the tremor-related thalamus activities were quenched when the network model exposed to the EDBS, while the tremor activities were increased when the network model exposed to the IDBS. This finding, again, confirmed that the role of the DBS is excitatory. Note that, we assigned the extra spiking activity of the thalamus neurons in the network model during the PD state to the involuntary tremor activity of the limbs, which was consistent with tremor frequency. By applying the EDBS on STN neurons, the tremor frequency was quenched, while in 76 has been reported by applying STN DBS the frequency of limb activity was increased and the amplitude of the limb movement was suppressed. In fact, the relation between the amplitude of the limb tremor and thalamus activity was not clarified in our network model.

Conclusion
We utilized a computational model of the BG to investigate the underlying mechanisms of the DBS. With the help of the model we concluded that first, the rebound burst activity related to T-type calcium current of the neurons has a key role for the generation of β oscillations. Second, we found that the role of the STN DBS is excitatory (and not inhibitory), because, the excitation of the STN neurons suppressed their rebound burst. Third, again, the rebound burst of thalamus neurons gave rise to the generation of resting tremor. Next, with exposing the STN neuron to high frequency excitation the thalamus tremor activity has quenched, while, with exposing the SNT neurons to high frequency inhibition it has been worsened. In summary, based on our model, we conclude that the role of the DBS is excitatory on its target.