# Dynamical Analyses on Beta Oscillations in a STN-GPE-GPI Model of Parkinson’s Disease

Authors: Yuanhong Bi; Quansheng Liu; Jingyi Zhao; Wuritu Yang
Journal: Complexity
 (2020)

---

## Abstract

Exploring the behaviors of beta oscillations in the basal ganglia is helpful to understand the mechanism of Parkinson’s disease. Studies have shown that the external and internal segments (GPe, GPi) of the globus pallidus receive different intensities of signals from the striatum in Parkinson’s disease and play different roles in the production of beta oscillations, but the relevant mechanism still remains unclear. Based on a model of the subthalamic nucleus (STN) and globus pallidus (GP), we propose an extended STN-GPe-GPi model and analyze the dynamical behaviors of beta oscillations in this model. The stability condition is obtained through theoretical analyses, and the generation of beta oscillations by the inputs from the cortex and striatum is further considered. The influence of some parameters related to GPi on its firing rate oscillations is discussed. The results obtained in this paper are expected to play a guiding role in the medical treatment of Parkinson’s disease.


---

## Body

## 1. Introduction

Parkinson’s disease is a chronic neurodegenerative disease with the symptoms of involuntary tremor of the hand and head, muscle rigidity, slow movement, and imbalance of posture [1]. The main pathological causes of Parkinson’s disease ascribe to the loss of dopamine neurons in the basal ganglia [2], which consists of the striatum, globus pallidus (GP), subthalamic nucleus (STN), compacta (SNc), and reticular (SNr) structure of the substantia nigra [3–5]. The loss of dopamine neurons causes beta oscillations with frequencies ranging from 13 Hz to 30 Hz in the basal ganglia [6, 7]. Therefore, it is necessary for understanding the mechanism of Parkinson’s disease to analyze the conditions of beta oscillations in the basal ganglia [8].Many researches explored the origin of beta oscillations in the basal ganglia [9–12]. Van Albada et al. believed that oscillations originated from the cortical-thalamic loop and then spread to the basal ganglia with the development of the disease [13]. Holgado et al. found that the STN-GP loop in the basal ganglia plays an important role in generating oscillations, which are related to connection weight and synaptic transmission time between STN and GP [14]. Furthermore, a model with two STN and one GP populations are considered to get the stability boundary of oscillations [15, 16]. However, GP population has not been divided into the external and internal segment (GPe and GPi) in the above models, where GPi is the main output structure of the basal ganglia and also used to treat dystonia by deep brain electrical stimulation in medicine [5]. Actually, GPi and GPe are affected by different intensities from striatum to result in Parkinson’s disease. Therefore, it is necessary to add both GPe and GPi into the basal ganglia network.Based on the above considerations, we introduce both GPe and GPi in the above model proposed by Holgado et al. as a new STN-GPe-GPi model. The mechanism of generating beta oscillations for the new model is explored through theoretical analyses and numerical simulation. The model of STN-GPe-GPi loop is given in Section2. Section 3 shows the results. Stability analyses and bifurcation for this model are given through theory analyses. Also, the effect of inputs from the cortex and striatum and some parameters related to GPi on beta oscillations is discussed in Section 4. Finally, the conclusion is given in Section 5.


## 2. Model

The model of STN- GPe-GPi loop，which is an extended STN-GP model, is considered here in order to understand the mechanism of Parkinson’s disease. Actually, GP can be divided into two parts, GPe and GPi. They can receive different excitatory inputs from the cortex in Parkinson’s disease, and GPi as the output part of the basal ganglia mainly affects the activity of neurons in thalamic and cortical areas. GPe and GPi receive inhibitory and excitatory stimulations from striatum and STN, respectively, but GPe sends inhibitory signal to GPi and STN. Besides, STN also receives excitatory stimulation from the cortex (see Figure1).Figure 1
Schematic diagram of STN-GPe-GPi. Each rectangular represents neural population, arrows indicate excitatory inputs, lines ending with black dots mean inhibitory inputs, and the striatum and cortex represent their inputs, respectively.The following firing rate equations are used to describe the dynamical behavior of the model, as shown in Figure1 [17–18]：(1)τGG1′t=FG−WGGG2t−TGG+WSGISt−TSGI−WXGStr−G1t,(2)τGG2′t=FGWSGESt−TSGE−WXGStr−G2t,(3)τSS′t=FS−WGSG2t−TGS+WCSCtx−St,where G1t, G2t, and St represent the firing rate of GPi，GPe, and STN; τG and τS are time constants of GP and STN; W and T are synaptic connection weight and the delay of signal transmission between neural populations, respectively, Wij and Tij are connection weight and time delay between the neural populations i and j, WSGE, WSGI, TSGE, and TSGI denote connection weights and the time delays from STN to GPe and GPi, respectively, Str and Ctx are input constants from the striatum and cortex. FS and FG are activation functions of STN and GP, which are given by the following formulas [14]:(4)FSx=MS1+MS−BS/BSexp−4x/MS,FGx=MG1+MG−BG/BGexp−4x/MG,where Mi is the maximum firing rate of neuron population i and Bi is the firing rate of neuron population i with no input. The activation functions FSx and FGx with their derivatives are shown in Figures 2(a) and 2(b), respectively.Figure 2
(a) Curves of the activation functionsFSx and FGx. (b) Derivatives of activation functions FSx and FGx.
(a)(b)The parameters and their source are given in Table1, and the connection weights between neuron populations in healthy and disease states are given in Table 2. Figure 3 shows time series of the firing rate of healthy and disease states, which oscillates for Parkinson’s disease and reaches the steady state in the healthy state.Table 1
Parameter values in the model and their sources.
ParameterValueReferenceT6 ms[19, 20]τS6 ms[21–23]τG14 ms[24]Ctx27 spk/s[25]Str2 spk/s[26]MS300 spk/s[27]BS17 spk/s[27]MG400 spk/s[20, 28]BG75 spk/s[28, 29]Table 2
Connection weights between neurons in different states.
ParameterHealthy stateParkinson’s diseaseWGS1.1210.7WSGE19.020.0WSGI19.020.0WGG6.6012.3WCS2.429.2WGG15.1139.4Figure 3
Time series of systems (1)–(3): (a) Parkinson’s disease and (b) healthy state.
(a)(b)


## 3. Results and Discussion

### 3.1. Stability Analyses for the STN-GPE-GPI Model

In this section, we obtain the following equations (5)–(7) by simplifying equations (1)–(3) with Fx=x without inputs from the striatum and cortex and the identical τ and T:(5)τG1′t=−WGGG2t−T+WSGISt−T−G1t,(6)τG2′t=WSGESt−T−G2t,(7)τS′t=−WGSG2t−T−St.Equations (5)–(7) are given in matrix form as follows:(8)G1′tG2′tS′t+0WGGτ−WSGIτ00−WSGEτ0WGSτ0G1t−TG2t−TSt−T+1τ0001τ0001τG1tG2tSt=0.Let(9)A=0WGGτ−WSGIτ00−WSGEτ0WGSτ0,B=1τ0001τ0001τ.Using Laplace transform [30],(10)Lft−T=e−sTFs,Lft=Fs,Lf′t=sFs−F0.Equation (8) can be given as follows:(11)sG1sG2sSs−G10G20S0+Ae−sTG1sG2sSs+BG1sG2sSs=0.Without loss of generality, let(12)G10G20S0=0.Then,(13)sI+Ae−sT+B=0,where S is the eigenvalue of the characteristic equation and I is the unit matrix. Taking the determinant on both sides of matrix (13),(14)detsI+Ae−sT+B=0,then(15)s+1τ3+WGSWSGEτ2s+1τe−2sT=s+1τs+1τ2+WGSWSGEτ2e−2sT=0.s≠−1/τ, so(16)s+1τ2+WGSWSGEτ2e−2sT=0.Lets=iλ and e−2sT be expanded by Euler transformation, and equation (16) is changed into(17)−λ2+1τ2+2iλτ+WGSWSGEτ2cos2λT−isin2λT=0.For simplification, letτ=1 and T^=T/τ，equation (17) can be written as follows：(18)−λ2+1+2iλ+WGSWSGEcos2λT^−isin2λT^=0.Let the real part and the imaginary part be zero, respectively:(19)−λ2+1+WGSWSGEcos2λT^=0,(20)2λ−WGSWSGEsin2λT^=0.Add the squares of (19) and (20) to get the following equation:(21)λ=WGSWSGE−1.Bring equation (21) into equation (19), and we get stability boundary of the linear system：(22)T^=12WGSWSGE−1arccos1−2WGSWSGE.So, the linear model equations (5)–(7) oscillate at the following condition:(23)Tτ>12WGSWSGE−1arccos1−2WGSWSGE.


### 3.2. Stability Analysis for the Nonlinear Model

For the nonlinear model with activation functions, we linearize the activation functions at the steady stateG1∗,G2∗,S∗ and get the characteristic equation (24) in a matrix form based on the steps in Section 3.1.(24)detsI+A1e−sT^+B=0,where(25)A1=0−FG1∗′WGGτFG1∗′WSGIτ00FG2∗′WSGEτ0−FS∗′WGSτ0,B=1τ0001τ0001τ,FG1∗′=FG1′−WGGG2∗+WSGIS∗−WXGStr,FG2∗′=FG2′WSGES∗−WXGStr,FS∗′=FS∗−WGSG2∗+WCSCtx.Therefore, the stability boundary of the nonlinear model is as follows:(26)T^=12FS∗′WGSFG2∗′WSGE−1arccos1−2FS∗′WGSFG2∗′WSGE.The oscillation condition is(27)Tτ>12FS∗′WGSFG2∗′WSGE−1arccos1−2FS∗′WGSFG2∗′WSGE.We analyze the stability of the model shown in Figure1 and get the stability conditions equations (22) and (26) for the linear and nonlinear model, respectively. Next, we draw the stability boundary curves of equations (22) and (26) and explore the effect of the inputs from the striatum and cortex on oscillation through numerical simulation.


### 3.3. Numerical Simulation of Stability Conditions

Figures4(a) and 4(b) describe the stability boundary curves of linear and nonlinear models based on equations (22) and (26), respectively. The decreasing boundary curve is infinitely close to x-axis with the increase of WGSWSGE⟶∞. The system oscillates for parameter values above the curve while it is stable for ones below the curve. However, the nonlinear system oscillate for a larger weight WGSWSGE than the one of linear system at the same T/τ.Figure 4
(a) Stability boundary curve for the linear system; (b) stability boundary for the nonlinear system; (c) the effect of cortical input on oscillation; (d) the effect of the striatum input on the oscillation. Ctx and Str denote the inputs from the cortex and striatum.
(a)(b)(c)(d)Figures4(c) and 4(d), respectively, show the influence of the inputs from the cortex and striatum on the stability of the nonlinear system. The boundary curve of the cortex is in the shape of “U,” where the system will oscillate in the area of above “U” and it will reach a stable steady state in the area below “U.” While the boundary curve of the striatum decreases to x-axis with the increase of the input of the striatum.


## 3.1. Stability Analyses for the STN-GPE-GPI Model

In this section, we obtain the following equations (5)–(7) by simplifying equations (1)–(3) with Fx=x without inputs from the striatum and cortex and the identical τ and T:(5)τG1′t=−WGGG2t−T+WSGISt−T−G1t,(6)τG2′t=WSGESt−T−G2t,(7)τS′t=−WGSG2t−T−St.Equations (5)–(7) are given in matrix form as follows:(8)G1′tG2′tS′t+0WGGτ−WSGIτ00−WSGEτ0WGSτ0G1t−TG2t−TSt−T+1τ0001τ0001τG1tG2tSt=0.Let(9)A=0WGGτ−WSGIτ00−WSGEτ0WGSτ0,B=1τ0001τ0001τ.Using Laplace transform [30],(10)Lft−T=e−sTFs,Lft=Fs,Lf′t=sFs−F0.Equation (8) can be given as follows:(11)sG1sG2sSs−G10G20S0+Ae−sTG1sG2sSs+BG1sG2sSs=0.Without loss of generality, let(12)G10G20S0=0.Then,(13)sI+Ae−sT+B=0,where S is the eigenvalue of the characteristic equation and I is the unit matrix. Taking the determinant on both sides of matrix (13),(14)detsI+Ae−sT+B=0,then(15)s+1τ3+WGSWSGEτ2s+1τe−2sT=s+1τs+1τ2+WGSWSGEτ2e−2sT=0.s≠−1/τ, so(16)s+1τ2+WGSWSGEτ2e−2sT=0.Lets=iλ and e−2sT be expanded by Euler transformation, and equation (16) is changed into(17)−λ2+1τ2+2iλτ+WGSWSGEτ2cos2λT−isin2λT=0.For simplification, letτ=1 and T^=T/τ，equation (17) can be written as follows：(18)−λ2+1+2iλ+WGSWSGEcos2λT^−isin2λT^=0.Let the real part and the imaginary part be zero, respectively:(19)−λ2+1+WGSWSGEcos2λT^=0,(20)2λ−WGSWSGEsin2λT^=0.Add the squares of (19) and (20) to get the following equation:(21)λ=WGSWSGE−1.Bring equation (21) into equation (19), and we get stability boundary of the linear system：(22)T^=12WGSWSGE−1arccos1−2WGSWSGE.So, the linear model equations (5)–(7) oscillate at the following condition:(23)Tτ>12WGSWSGE−1arccos1−2WGSWSGE.


## 3.2. Stability Analysis for the Nonlinear Model

For the nonlinear model with activation functions, we linearize the activation functions at the steady stateG1∗,G2∗,S∗ and get the characteristic equation (24) in a matrix form based on the steps in Section 3.1.(24)detsI+A1e−sT^+B=0,where(25)A1=0−FG1∗′WGGτFG1∗′WSGIτ00FG2∗′WSGEτ0−FS∗′WGSτ0,B=1τ0001τ0001τ,FG1∗′=FG1′−WGGG2∗+WSGIS∗−WXGStr,FG2∗′=FG2′WSGES∗−WXGStr,FS∗′=FS∗−WGSG2∗+WCSCtx.Therefore, the stability boundary of the nonlinear model is as follows:(26)T^=12FS∗′WGSFG2∗′WSGE−1arccos1−2FS∗′WGSFG2∗′WSGE.The oscillation condition is(27)Tτ>12FS∗′WGSFG2∗′WSGE−1arccos1−2FS∗′WGSFG2∗′WSGE.We analyze the stability of the model shown in Figure1 and get the stability conditions equations (22) and (26) for the linear and nonlinear model, respectively. Next, we draw the stability boundary curves of equations (22) and (26) and explore the effect of the inputs from the striatum and cortex on oscillation through numerical simulation.


## 3.3. Numerical Simulation of Stability Conditions

Figures4(a) and 4(b) describe the stability boundary curves of linear and nonlinear models based on equations (22) and (26), respectively. The decreasing boundary curve is infinitely close to x-axis with the increase of WGSWSGE⟶∞. The system oscillates for parameter values above the curve while it is stable for ones below the curve. However, the nonlinear system oscillate for a larger weight WGSWSGE than the one of linear system at the same T/τ.Figure 4
(a) Stability boundary curve for the linear system; (b) stability boundary for the nonlinear system; (c) the effect of cortical input on oscillation; (d) the effect of the striatum input on the oscillation. Ctx and Str denote the inputs from the cortex and striatum.
(a)(b)(c)(d)Figures4(c) and 4(d), respectively, show the influence of the inputs from the cortex and striatum on the stability of the nonlinear system. The boundary curve of the cortex is in the shape of “U,” where the system will oscillate in the area of above “U” and it will reach a stable steady state in the area below “U.” While the boundary curve of the striatum decreases to x-axis with the increase of the input of the striatum.


## 4. The Effect of Parameters Related to GPI on Frequency and Amplitude of Nonlinear Model Oscillation

In this section, the effect of three groups of parameters related to GPi,WGG−TGG, WSGI−TSGI, and WSGE−WSGI, on oscillation frequency (Figures 5 (a1)–(a3)) and amplitude (Figures 5 (b1)–(b3)) are also considered. Besides, time series of the GPI firing rate are given in Figure 6 for typical parameter values in each group of parameters in order to clearly see the influence of these parameters on the frequency and amplitude of GPi oscillation. We set the frequency of the stable steady state to be 0. As can be seen from Figures 5 (a1) and (b1), the plane is divided into two parts, where the firing rate of GPi reaches a steady state for parameters taken in the blue area of the upper half while it will oscillate for parameters taken in the red area of the lower half. Frequency of oscillation with 14 Hz almost is not affected by WGG and TGG, while smaller WGG and larger TGG increase the amplitude of oscillation to 180. Figures 5 (a2) and (b2) shows that the firing rate of GPi always oscillates with beta frequency band of about 14 Hz and lower amplitudes for WSGI<20 otherwise alpha oscillations with frequency of 8 Hz–12 Hz and higher amplitude for WSGI>20. Besides, Figure 5 (a3) shows that the firing rate of GPi reaches the steady state for smaller WSGE (WSGE<3.0) while it oscillates for large WSGE (WSGE>3.0), where the frequency and amplitude increase for smaller WSGE and larger WSGI.Figure 5
The effect ofWGG−TGG, WSGI−TSGI, and WSGE−WSGI on the frequency of the GPi firing rate (a1)–(a3) and the amplitude of GPi (b1)–(b3). Other parameters are taken from Parkinson’s state in Table 2.Figure 6
Typical time series of the firing rate of GPi for some parameters: (a1)–(a3) WGG−TGG, (b1)–(b3) WSGI−TSGI, and (c1)–(c3) WSGE−WSGI.Furthermore, time series of the GPi firing rate are given in Figure6 for three pairs of parameters from each group of parameter plane in Figure 5 to clearly see the influence of these parameters on the frequency and amplitude of the GPi firing rate. It can be seen from Figures 6 (a1)–(a3) that the firing rate of GPi oscillates for smaller WGG and larger TGG. According to Figures 6 (b1)–(b3), the firing rate of GPi always oscillates with smaller amplitude and larger frequency for smaller WSGI and TSGI while the case is opposite for larger WSGI and TSGI. Figures 6 (c1)–(c3) show that larger WSGE makes the firing rate of GPi oscillate and larger WSGI increases the amplitude of oscillation.To sum up, the connection weights between neurons have a great influence on oscillation of the GPi firing rate, where smallerWGG and larger WSGE, WSGI easily induce the oscillation of the GPi firing rate. However, transmission delays between neurons affect the amplitude of the GPi firing rate, where larger TGG and smaller TSGI will increase the amplitude.


## 5. Conclusions

Analyzing the conditions of beta oscillation in the basal ganglia is helpful to understand the mechanism of Parkinson’s disease. In this paper, we analyze the conditions of beta oscillation in the STN-GPe-GPi model for different cases. First, stability analyses give stability boundary condition equations (22) and (26) for the linear system and nonlinear system, respectively, which are shown in Figures 4(a) and 4(b). The nonlinear system oscillates for larger connection weight than the one of the linear system due to activation function. In addition, we consider the influence of the cortex and striatum as external input of the STN-GPe-GPi loop on the stability boundary of the nonlinear system. As can be seen from Figures 4(c) and 4(d), the nonlinear system is in the state of oscillation for 100<ctx<400 or 5<str<25. Comparing with the results in [15], the stability boundary system in this model moves to the right and other results almost are consistent with the ones in [15]. Furthermore, it can be seen from numerical simulation that the influence of the connection weights and delays is related to GPi on its oscillation. Smaller connection weight from GPe to GPi (WGG) and larger one from STN to GPe (WSGE) make the system oscillate easily, regardless of connection weight (WSGI) and transmission delay (TSGI) from STN to GPi. We hope that the results may provide guidance for the therapy of reducing pathological oscillations of PD, especially for the operation of the region related to GPi. However, it is necessary to consider a more complete neural network related to Parkinson’s disease and explore the conditions of beta oscillation in response to different time delay, noise, and temperature [31]. Furthermore, we will investigate the pathogenesis of Parkinson’s disease from the perspective of systems’ biology in the future [32].



---

*Source: 4376279-2020-09-09.xml*