The mechanism of Parkinson oscillation in the cortex: Possible evidence in a feedback model projecting from the globus pallidus to the cortex

1.   Introduction

Parkinson's disease is the second most common neurodegenerative brain disease and has a higher incidence rate among the elderly ^[1,2]. Parkinson's disease places a great burden on society and the economy, and its research has attracted more and more attention from society. Postural imbalance, rigidity, tremors, etc. are typical symptoms of Parkinson's patients ^[3,4]. The loss of neurons in the substantia nigra pars compacta, which causes a large loss of dopamine in the striatum, is considered to be the main cause of these symptoms ^[5,6]. The precise pathological mechanism of dopamine loss is still unclear, and many different views have been formed. Moreover, the excessive beta frequency synchronous oscillation observed in electroencephalogram is generally considered to be closely related to the emergence and disappearance of Parkinson's symptoms ^[7,8,9].

In the past, many brain network computing models have been developed to explore the origin mechanism of beta oscillation ^{[10,11,12,13,14,15,16,17,18]}. Especially, many studies have shown that beta oscillations in the basal ganglia may be induced by anomalous interaction of circuits that consist of the subthalamic nucleus (STN) and the globus pallidus pars externa (GPe) ^{[19,20,21,22,23,24]}. The STN contains a large number of glutamatergic neurons that exert excitatory inputs to the GPe; in turn, the GPe is mainly composed of GABAergic neurons, which give inhibitory feedback to the STN. Together, they form an inhibitory-excitatory coupling loop. Nevado Holgado et al. pointed out that the STN-GPe network can produce beta oscillation independently, and coupling weights and delays in the loop have great influence on oscillation ^[19]. Recently, Hu et al. explored the bidirectional Hopf bifurcation mechanism in an STN-GPe model, which can uniformly explain the origin of the observed oscillation ^[24].

The cortex has a close input-output relationship with the basal ganglia. The cortex exerts strong excitatory projections to the basal ganglia, which is necessary for the existence of Parkinson's oscillation ^[25,26]. In fact, the beta oscillation may also originate in the cortex and spread to other brain regions ^[27,28]. The loss of dopamine in the striatum can disrupt cortical function, which may be the cause of Parkinson's symptoms ^[29]. The function of the primary motor cortex is damaged in Parkinson's disease, which may be an effective goal for neurostimulation in alleviating motor dysfunction ^[30]. In anatomical structure, the cortex is mainly composed of long-range excitatory pyramidal neurons (CEX) and short-range inhibitory interneurons (CIN) ^[27,28]; they also form excitatory-inhibitory circuits, the structure of which allows easy generation of oscillations independently ^[31]. Pavlides et al. proposed a STN-GPe-CEX-CIN mean-field network ^[32]. In this model, the CEX-CIN loop exerts excitatory inputs to the STN. They found that beta oscillations can occur in the cortex, and the STN-GPe loop resonated with inputs from the cortex. Recently, Chen et al. obtained theoretical critical conditions for Hopf bifurcation of beta oscillation in a similar model ^[33]. We note that the mechanism of beta oscillation in the cortex is rarely discussed in previous models.

In most previous models, the cortex had strong excitability projection to the subthalamic nucleus, and the effect of cortical outputs on beta oscillations has been widely studied ^{[20,27,28,34,35,36]}. In fact, there is a significant inhibitory feedback effect from basal ganglia to cortex ^[32], which also has a great influence on the origin of Parkinson's oscillation, but there are few relevant theoretical models. Recently, experiments have shown that GPe has a direct inhibitory output on both excitatory and inhibitory neurons in the cortex ^[37,38]. In addition, deep brain stimulation technology has a significant effect on GPe to alleviate the symptoms of Parkinson's disease ^[39,40,41], which may be regulated through these direct outputs to the cortex. Because many regions of the cortex have been reported to be closely related to the generation of Parkinson oscillations ^[42,43,44], we are interested to know whether the direct inhibitory projection from the basal ganglia to the cortex can explain the origin and disappearance of beta oscillation. Computational models can provide some possible mechanisms for these pathways in inducing and controlling oscillations, which can be confirmed in experiments.

In this paper, we propose a new model that consists of the STN, the GPe, the CEX and the CIN nucleus. In this model, there are direct inhibitory inputs from GPe to CEX and CIN, respectively. We mainly study the critical conditions of oscillation, Hopf bifurcation mechanism, amplitude and frequency characteristics of oscillation in the CEX-CIN circuit.

The new model is described in Section 2, and the stability analysis is given in Sections 3 and 4. The theoretical mechanism and main results are obtained in Sections 5 and 6. Finally, the conclusions are summarized in Section 7.

2.   Model description

Figure 1 is our new model, which consists of cortex and basal ganglia. The cortex includes excitatory nuclei (CEX) and inhibitory nuclei (CIN), and the basal ganglia is composed of subthalamic nucleus (STN) and globus pallidus external (GPe) segment. C represents the constant excitatory input to CEX and STN, and Str represents the inhibitory input from the striatum. The arrow represents the excitatory projection regulated by glutamate, and the dot represents the inhibitory projection regulated by GABAergic. There is a self-feedback connection inside GPe.${T}_{ij}\;\mathrm{a}\mathrm{n}\mathrm{d}\;{w}_{ij}$ represent the time delay of signal transmission and coupling strength on the pathway, respectively, and subscripts i and j indicate the origin and termination positions of the signal, respectively. CEX-CIN and STN-GPe loops are both "excitatory-inhibitory" coupling structures. In this model, there are two inhibitory feedback projections from the GPe to cortex, which is different from many previous studies. For the convenience of the following description, CEX, CIN, STN and GPe are represented by marks E, I, S and G respectively.

We use the following average discharge rate equation to describe the coupling relationship and dynamic behavior between neuron groups in Figure 1 ^{[19,20,32,45,46]}:

where S(t), G (t), E(t) and I(t) represent the average discharge rates of STN, GPe, CEX and CIN, respectively; S(t-T), G(t-T), E(t-T) and I(t-T) represent the corresponding delay discharge rates. ${\tau }_{X}$ is the time constant corresponding to neuron group X. ${F}_{Y}$(Y = S, G, E, I) is the activation function of population Y, showing the effect of synaptic input on the average discharge rate ^[19,20,32].

syt in Eqs (5)–(8) represents all possible inputs of the corresponding neuron group, ${M}_{X}$ is the maximum discharge rate of neuron group X, and ${B}_{X}$ is the basic discharge rate of neuron group X without external input. Similar to previous research ^[19,20,32], we show the curves of the activation function and its derivatives in Figure 2A, B, respectively. The activation functions of different neurons are heterogeneous, and all of them are S-shaped curves. The maximum derivative of the activation functions is 1.

All parameters and values used in this paper are listed in Tables 1 and 2, and they come from previous experiments and model studies ^{[19,20,32,33,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62]}. Figure 3A, B show the simulated discharge rates in Parkinson's state and healthy state, respectively. In Parkinson's state, different neural nuclei show synchronous oscillation (OS). In the healthy state, the system is stable (SS), and the discharge rate of each neuron tends to be constant. Stability and oscillation are two basic states of the model.

3.   Stability analysis of linearized system

First, let ${F}_{j}\left(\mathrm{s}\mathrm{y}\mathrm{t}\right) = \mathrm{s}\mathrm{y}\mathrm{t}\left(\mathrm{j} = \mathrm{S}, \mathrm{G}, \mathrm{E}, \mathrm{I}\right)$ to obtain the linearized system (9)–(12) of Eqs (1)–(4). In this part, we use the Laplace transform and the theory of solutions of ordinary differential equations to obtain the oscillation critical conditions of system (9)–(12). For the feasibility of mathematical derivation, we assume that the time delay and time constant are the same for all projections and neurons, and they are recorded as T and ${\rm{ \mathsf{ τ} }}$ respectively. Then, the following linear system is obtained:

By writing Eqs (9)–(12) in matrix form and shifting the terms, we can get

where

The following operations and properties of the Laplace transform will be used ^[63]:

By performing the corresponding Laplace transform on the matrices other than A and B in Eq (13), we can get

Without losing generality, we set the initial conditions as

Then,

where s is the eigenvalue of the characteristic equation, and I is the identity matrix. Take the determinant on both sides of Eq (21):

i.e.,

Then,

Let

From the knowledge of ordinary differential equations, we know that the stability of system (9)–(12) is determined by the sign of the real part of eigenvalue s in Eq (24). When the real part of s is greater than 0, the system is unstable; when the real part of s is less than 0, the system is stable. Therefore, if we find the case that s is a pure imaginary number, we can get the stable and unstable transition boundary of the system, that is, the stability boundary. When ${\Delta }_{1}$ or ${\Delta }_{2}$ is equal to 0, Eq (24) holds. Next, we will discuss these two situations.

(ⅰ) Assume the characteristic value $\mathrm{s} = \mathrm{i}{\rm{ \mathsf{ λ} }}$ (${\rm{ \mathsf{ λ} }} > 0$). Substitute the eigenvalue into Eq (25) and expand it by the Euler formula:

Let ${\rm{ \mathsf{ τ} }} = 1, \; {\mathrm{T}}_{1} = \frac{\mathrm{T}}{{\rm{ \mathsf{ τ} }}}$, and then Eq (27) becomes

Separating the real part and imaginary part of Eq (28), we can obtain

Add the squares of Eqs (29) and (30), and we can get that

The stability boundary condition can be obtained by substituting ${\rm{ \mathsf{ λ} }}$ into Eq (29):

Then, the conditions for oscillation of the linear system are as follows:

(ⅱ) Assume characteristic value $\mathrm{s} = \mathrm{i}{\rm{ \mathsf{ λ} }}$ (${\rm{ \mathsf{ λ} }} > 0$). We substitute the eigenvalue into Eq (26) and expand it by the Euler formula:

Set ${\rm{ \mathsf{ τ} }} = 1, {\mathrm{T}}_{2} = \frac{\mathrm{T}}{{\rm{ \mathsf{ τ} }}}$, and then Eq (34) can be written as

Next, for the convenience and feasibility of derivation, we assume that the self-feedback connection of GPe is 0. That is, when ${w}_{GG} = 0$, Eq (35) is reduced to

By separating the real part and imaginary part of Eq (36), Eqs (37) and (38) can be obtained:

Add the squares of Eqs (37) and (38) and find the expression of ${\rm{ \mathsf{ λ} }}$,

Replace ${\rm{ \mathsf{ λ} }}$ in Eq (37) with (39), and the stability boundary condition can be obtained:

Then, the conditions for oscillation are as follows:

4.   Mathematical analysis of stability conditions of nonlinear systems

In this part, we use the results of oscillation boundary conditions of linear systems to analyze the stability boundary conditions of nonlinear systems. First, the equilibrium point of the nonlinear system can be obtained by minimizing Eq (42), which is recorded as $\left({\mathrm{S}}^{\mathrm{*}}, {\mathrm{G}}^{\mathrm{*}}, {\mathrm{E}}^{\mathrm{*}}, {\mathrm{I}}^{\mathrm{*}}\right)$.

The characteristic equation of the new system can be obtained by linearizing the activation function at the equilibrium point $\left({\mathrm{S}}^{\mathrm{*}}, {\mathrm{G}}^{\mathrm{*}}, {\mathrm{E}}^{\mathrm{*}}, {\mathrm{I}}^{\mathrm{*}}\right)$:

where

The expressions of ${\mathrm{F}}_{{\mathrm{S}}^{\mathrm{*}}}^{\mathrm{'}}, {\mathrm{F}}_{{\mathrm{G}}^{\mathrm{*}}}^{\mathrm{'}}, {\mathrm{F}}_{{\mathrm{E}}^{\mathrm{*}}}^{\mathrm{'}}, {\mathrm{F}}_{{\mathrm{I}}^{\mathrm{*}}}^{\mathrm{'}}$ are as follows:

Then, using the same steps as dealing with linear systems (Section 3), we can obtain two stability boundary conditions corresponding to nonlinear systems:

The corresponding two oscillation conditions are

5.   Mechanism analysis of nonlinear oscillation boundary conditions

In this part, we mainly analyze and study some possible generation mechanisms of the stability critical condition (50).

5.1.   The effect of the excitatory-inhibitory coupling loop E-I on oscillation

The E and I form an excitatory-inhibitory coupling circuit in this model, and this structure is prone to oscillation. In this section, we study the mechanism of coupling weights in this loop on oscillations.

Figure 4 describes the influence of coupling weight ${w}_{EI}$ on oscillation. Figure 4A shows the theoretical oscillation critical conditions obtained by simulating Eq (50). The system is oscillating above the curve and stable below the curve. Fix parameter $\frac{T}{\tau }$ at an appropriate value (such as 4), and with the increase of ${w}_{EI}$, the state of the system can transition from stable to oscillation (as shown by the arrow in Figure 4A). The influence mechanism of ${w}_{EI}$ on nuclei E and I is different. From the connection structure in Figure 1, we know that the excitability coupling weight ${w}_{EI}$ acts directly on I, and its increase directly promotes the improvement of discharge ability of I. Therefore, with the increase of ${w}_{EI}$, the ADR of I gradually increases, as shown in Figure 4B. When ${w}_{EI}$ increases to a certain value, Hopf bifurcation (HP1) occurs in the state of I, and oscillation occurs in I, as shown in Figure 4C. For the convenience of the following description, we define HP1 to represent the Hopf bifurcation transition between the ordinary stable state (SS) and the oscillating state (OS). ${w}_{EI}$ mainly affects the discharge state of E through pathway $\mathrm{E}\to \mathrm{I}\to \mathrm{E}$, and we can infer that this effect is generally inhibitory from Figure 1. Therefore, with the increase of ${w}_{EI}$, the ADR of E gradually decreases, as shown in Figure 4D. When ${w}_{EI}$ is relatively small, the ADR of E is relatively high, and the E is in a stable state of high discharge rate (HSS). Interestingly, when ${w}_{EI}$ increases to a certain value, Hopf bifurcation phenomenon (HP2) will also appear in E, which represents the transition between the stable state of high discharge rate and oscillation state. Figure 4F is the stable state obtained by taking ${w}_{EI}$ = 8; Figure 4G is the time series diagram of oscillation state obtained by taking ${w}_{EI}$ = 42.

Figure 5 shows the influence of coupling weight ${w}_{IE}$ on oscillation in the system. Figure 5A shows the theoretical oscillation critical condition obtained from Eq (50). The system is oscillating above the curve and stable below the curve. The oscillation critical curve presents a U-shape, and when $\frac{T}{\tau }$ is fixed at a value (such as 2), with the increase of ${w}_{IE}$, the state of the system will be transferred from HSS to OS and finally to SS, as shown by the arrow in Figure 5A. ${w}_{IE}$ is the inhibitory coupling weight directly acting on E, and the increase of ${w}_{IE}$ will have a strong inhibitory effect on the activity of E. ${w}_{IE}$ mainly affects the discharge activity of I through indirect pathway "$\mathrm{I}\to \mathrm{E}\to \mathrm{I}$", and we can infer from Figure 1 that the increase of ${w}_{IE}$ also exerts an inhibitory effect on I. Therefore, with the increase of ${w}_{IE}$, the ADRs of nuclei E and I are decreasing, as shown in Figures 5B and 5C. In Figure 5D, E, we simulate the state bifurcation diagrams of E and I, respectively. When ${w}_{IE}$ is relatively small, the E-I circuit is in HSS; with the increase of ${w}_{IE}$, its inhibitory effect on E and I gradually increased. When ${w}_{IE}$ increases to a certain value, HP2 bifurcation occurs in E-I loop, and HSS transitions to OS. If ${w}_{IE}$ is too large, the oscillation state will disappear, and HP1 bifurcation will appear in E and I. The two-way bifurcation phenomenon shown in Figure 5D, E is consistent with the U-shaped curve in Figure 5A. Figure 5F–H describe three different states of loop E-I. Figure 5F shows the HSS state, with a relatively high average discharge rate. Figure 5G shows the OS state, and Figure 5H shows the SS state.

In addition to the coupling weight, time delay is another key parameter in the model. In Figures 4A and 5A, we can observe that the time delay must reach a certain value before oscillation can occur in the system. Figure 6 describes the effect of time delay on loop E-I oscillation. Figure 6A is a bifurcation diagram of the state of E with an increase in delay. Figure 6B is a state transition diagram of I with the change of delay. They show that when the time delay reaches a certain value, Hopf bifurcation will occur in E and I, and the oscillation amplitude will increase with the increase of time delay. Figure 6C is the stable time series obtained by taking T = 0.003 s. Figure 6D is the oscillation time series obtained by taking T = 0.007 s.

5.2.   Influence of direct and indirect projections on loop E-I oscillation

In this part, we study the influence of direct and indirect feedback inputs from STN and GPe on the oscillation in E and I.

5.2.1.   Influence of the direct inhibitory projections of GPe on oscillation

Figure 7 describes the effect of the direct inhibitory pathway "GPe→E" on oscillation. ${\mathrm{w}}_{\mathrm{G}\mathrm{E}}$ is the coupling weight on the pathway "GPe→E". Figure 7A is the oscillation critical curve simulated by Eq (50). The system is oscillating above the curve and stable below the curve. As shown by the arrow in Figure 7A, when we fix abscissa $\frac{\mathrm{T}}{{\rm{ \mathsf{ τ} }}}$ at an appropriate value (such as 7), with the increase of ${\mathrm{w}}_{\mathrm{G}\mathrm{E}}$, the system can realize bidirectional Hopf bifurcation transition between different stable states and oscillation states. ${\mathrm{w}}_{\mathrm{G}\mathrm{E}}$ is the inhibitory coupling weight directly acting on E; therefore, the increase of ${\mathrm{w}}_{\mathrm{G}\mathrm{E}}$ will lead to the decrease of ADR of E, as shown in Figure 7B. ${\mathrm{w}}_{\mathrm{G}\mathrm{E}}$ mainly affects the activities of I through the pathway "GPe→E $\to \mathrm{I}$". From Figure 1, we can infer that the increase of ${\mathrm{w}}_{\mathrm{G}\mathrm{E}}$ generally exerts an inhibitory effect on I; therefore, the ADR of I also decreases gradually, as shown in Figure 7C. Figure 7D, E show the bifurcation of different oscillation states of E and I, respectively. With the increase of ${\mathrm{w}}_{\mathrm{G}\mathrm{E}}$, loop E-I realizes the synchronous transfer of HSS, OS and SS. When ${\mathrm{w}}_{\mathrm{G}\mathrm{E}}$ is small, its inhibitory effect is relatively small, and E and I are in the HSS state. When ${\mathrm{w}}_{\mathrm{G}\mathrm{E}}$ increases to a certain value, HP2 bifurcation appears in loop E-I, and oscillation begins to appear. When ${\mathrm{w}}_{\mathrm{G}\mathrm{E}}$ is too large, its inhibitory effect on loop E-I is too strong, and the oscillation of E and I disappears with the appearance of HP1. Figure 7F–H show time series diagrams obtained by taking ${\mathrm{w}}_{\mathrm{G}\mathrm{E}}$ as 5, 70 and 100, respectively, and they visually show the different discharge states of loop E-I.

"GPe→I" is a direct inhibitory pathway applied from globus pallidus to I. ${\mathrm{w}}_{\mathrm{G}\mathrm{I}}$ is the coupling weight on the pathway "GPe → I". We simulate the effect of ${\mathrm{w}}_{\mathrm{G}\mathrm{I}}$ on oscillation in Figure 8A. As shown by the arrow in the figure, when the delay is fixed, the increase of ${\mathrm{w}}_{\mathrm{G}\mathrm{I}}$ can make the state of loop E-I transition from oscillation to stability, and the mechanism of oscillation in E and I is different. The GPe mainly acts on E through the pathway "GPe $\to \mathrm{I}\to \mathrm{E}$" and the two coupling weights on this pathway are inhibitory. Therefore, the increase of ${\mathrm{w}}_{\mathrm{G}\mathrm{I}}$ exerts an excitatory effect on E, and the ADR of E increases with the increase of ${\mathrm{w}}_{\mathrm{G}\mathrm{I}}$ (as shown in Figure 8B). When ${\mathrm{w}}_{\mathrm{G}\mathrm{I}}$ increases to a certain value, HP2 bifurcation will occur in E, the oscillation of E disappears, and E enters into HSS state, as shown in Figure 8C. On the contrary, ${\mathrm{w}}_{\mathrm{G}\mathrm{I}}$ exerts a direct inhibitory effect on I; therefore, the ADR of I decreases with the increase of ${\mathrm{w}}_{\mathrm{G}\mathrm{I}}$ (as shown in Figure 8D). When ${\mathrm{w}}_{\mathrm{G}\mathrm{I}}$ increases to a certain value, HP1 bifurcation will occur in I, the oscillation of I disappears, and I enters into SS state, as shown in Figure 8E. Figure 8F is a stable time series diagram of loop E-I, and Figure 8G is an oscillation time series diagram of loop E-I. They are obtained by taking ${\mathrm{w}}_{\mathrm{G}\mathrm{I}}$ = 38.1 and ${\mathrm{w}}_{\mathrm{G}\mathrm{I}}$ = 10 respectively.

5.2.2.   The effect of direct excitatory projection of cortex on oscillation

In this model, we assume that there is a constant excitatory input from cortex to E, and ${\mathrm{w}}_{\mathrm{C}\mathrm{E}}$ represents the coupling weight in the input. Figure 9 describes the influence of coupling weight ${\mathrm{w}}_{\mathrm{C}\mathrm{E}}$ on oscillation. Figure 9A shows the theoretical boundary conditions obtained from Eq (50). The system is oscillating above the curve and stable below the curve. The left side of the curve is SS state, and the right side of the curve is HSS state. When $\frac{\mathrm{T}}{{\rm{ \mathsf{ τ} }}}$ is fixed at a suitable value (such as 3), the increase of ${\mathrm{w}}_{\mathrm{C}\mathrm{E}}$ can realize the bidirectional Hopf bifurcation transition of loop E-I. ${\mathrm{w}}_{\mathrm{C}\mathrm{E}}$ acts directly on E, and its increase will certainly promote the improvement of discharge capacity of E. ${\mathrm{w}}_{\mathrm{C}\mathrm{E}}$ acts on I indirectly through the pathway "E→I" and"E→I" is an exciting pathway. Therefore, the increase of ${\mathrm{w}}_{\mathrm{C}\mathrm{E}}$ also promotes the improvement of the discharge capacity of I. Therefore, as shown in Figure 9B, C, the ADRs of E and I both increase with the increase of ${\mathrm{w}}_{\mathrm{C}\mathrm{E}}$. When ${\mathrm{w}}_{\mathrm{C}\mathrm{E}}$ is small, the discharge capacity of E-I is low, and the circuit is in a stable state. When ${\mathrm{w}}_{\mathrm{C}\mathrm{E}}$ increases to a certain value, HP1 bifurcation occurs in E and I, and the loop begins to oscillate. When ${\mathrm{w}}_{\mathrm{C}\mathrm{E}}$ is too large, the oscillation will disappear, and the loop will enter into HSS state. Figure 9D, E intuitively show the bidirectional Hopf bifurcation phenomenon of E and I.

5.2.3.   Effects of subthalamic nucleus and globus pallidus circuit on oscillations in E and I

The influence of subthalamic nucleus and globus pallidus circuit on E and I is indirect. In this part, we mainly study the influence mechanism of coupling weights related to STN-GPe loop on E-I loop oscillation.

${\mathrm{w}}_{\mathrm{G}\mathrm{S}}$ is an inhibitory coupling weight inside the loop STN-GPe. Figure 10 shows the effect of ${\mathrm{w}}_{\mathrm{G}\mathrm{S}}$ on oscillation. Figure 10A is the oscillation critical curve simulated by Eq (50). The system is oscillating above the curve and stable below the curve. As shown by the arrow in Figure 10A, when we fix $\frac{\mathrm{T}}{{\rm{ \mathsf{ τ} }}}$ at a suitable value (such as 10), the increase of ${\mathrm{w}}_{\mathrm{G}\mathrm{S}}$ can make the state of the system transition from stability to oscillation. From the connection structure in Figure 1, we can analyze that ${\mathrm{w}}_{\mathrm{G}\mathrm{S}}$ affects nuclei E and I through the pathways "GPe→STN→GPe→E" and "GPe→STN→GPe→I", respectively. From the competition of inhibition and excitation coupling weights in these two pathways, we can infer that the increase of ${\mathrm{w}}_{\mathrm{G}\mathrm{S}}$ can promote the discharge of E and I in general. Therefore, the ADR of E and I will increase with the increase of ${\mathrm{w}}_{\mathrm{G}\mathrm{S}}$, as shown in Figure 10B, C. When ${\mathrm{w}}_{\mathrm{G}\mathrm{S}}$ increases to a certain value, the discharge level of loop E-I exceeds a certain threshold, and HP1 bifurcation occurs in E and I, thus starting oscillation. Figure 10D, E are bifurcation diagrams of E and I with the increase of ${\mathrm{w}}_{\mathrm{G}\mathrm{S}}$. In order to observe the stable and oscillating states more intuitively, we take ${\mathrm{w}}_{\mathrm{G}\mathrm{S}}$ = 6 and ${\mathrm{w}}_{\mathrm{G}\mathrm{S}}$ = 85, respectively, to simulate the time series diagram, as shown in Figure 10F, G.

In this model, we assume that the striatum has an inhibitory constant input Str to GPe. Figure 11 shows the effect of Str on oscillation. Figure 11A is the oscillation critical curve simulated by Eq (50). The system is oscillating above the curve and stable below the curve. As shown by the arrow in Figure 11A, when we fix $\frac{T}{\tau }$ at a suitable value (such as 1), the increase of Str can make the state of the system transition from stability to oscillation. From the connection structure in Figure 1, we can analyze that Str affects nuclei E and I through the pathways "Str→GPe→E" and "Str→GPe→I", respectively. The coupling weights in these two pathways are all inhibitory, and we can speculate that the increase of Str can promote the discharge of E and I in general. Therefore, the ADR of E and I will increase with the increase of Str, as shown in Figure 11B, C. When Str increases to a certain value, the discharge level of loop E-I exceeds a certain threshold, and HP1 bifurcation occurs in E and I, thus starting oscillation, as shown in Figure 11D, E. In order to observe the stable and oscillating states more intuitively, we take Str = 3 and Str = 20 to simulate the time series diagrams, as shown in Figure 11F, G.

${\mathrm{W}}_{\mathrm{S}\mathrm{G}}$ is the excitability coupling weight inside the loop STN-Gpe. In Figure 12, we simulate the influence of ${\mathrm{w}}_{\mathrm{S}\mathrm{G}}$ on the oscillation in system E-I. The curve in Figure 12A is the critical condition of Hopf bifurcation, which is simulated by mathematical Eq (50). If the parameters fall above the curve, the system is oscillatory; if the parameters fall below the curve, the system is stable. As shown by the arrow in Figure 12A, when $\frac{\mathrm{T}}{{\rm{ \mathsf{ τ} }}}$ is fixed at a certain value (such as 0.8), the increase of ${\mathrm{w}}_{\mathrm{S}\mathrm{G}}$ can make the state of the system transition from oscillation to stability, and the mechanism of oscillation in E and I is different. It can be seen from Figure 1 that ${\mathrm{w}}_{\mathrm{S}\mathrm{G}}$ mainly affects the active state of I through the inhibitory pathway "Gpe→I" and the increase of ${\mathrm{w}}_{\mathrm{S}\mathrm{G}}$ will strengthen the inhibitory effect of "GPe→I". Therefore, the ADR of I decreases with the increase of ${\mathrm{w}}_{\mathrm{S}\mathrm{G}}$, as shown in Figure 12B. When ${\mathrm{w}}_{\mathrm{S}\mathrm{G}}$ increases to a certain value, the oscillation of I is suppressed, and HP1 bifurcation occurs in I, as shown in Figure 12C. On the other hand, ${\mathrm{w}}_{\mathrm{S}\mathrm{G}}$ affects the active state of E through two pathways: "GPe→E" and "GPe→I→E". Obviously, the increase of ${\mathrm{w}}_{\mathrm{S}\mathrm{G}}$ will strengthen the inhibitory effect of "GPe→E". However, the two coupling weights on the pathway "GPe→I→E" are both inhibitory, and the increase of ${\mathrm{w}}_{\mathrm{S}\mathrm{G}}$ generally strengthens the excitatory effect of "GPe→I→E". Because "I→E" is the internal pathway of loop E-I, we speculate that its impact on E is stronger than the external pathway. Therefore, through the competition between the pathways "GPe→E" and "GPe→I→E", the increase of ${\mathrm{w}}_{\mathrm{S}\mathrm{G}}$ generally exerts an excitatory effect on E, and the ADR of E also increases, as shown in Figure 12D. When ${\mathrm{w}}_{\mathrm{S}\mathrm{G}}$ reaches a certain value, the oscillation state of E shifts to HSS, and HP2 occurs, as shown in Figure 12E. Figure 12F, G show time series diagrams obtained by taking ${\mathrm{w}}_{\mathrm{S}\mathrm{G}}$ = 5 and ${\mathrm{w}}_{\mathrm{S}\mathrm{G}}$ = 20, respectively.

In this model, we assume that there is an inhibitory self-feedback projection in the GPe, and ${\mathrm{w}}_{\mathrm{G}\mathrm{G}}$ is the coupling weight in this projection. Obviously, ${\mathrm{w}}_{\mathrm{G}\mathrm{G}}$ affects the discharge activities of nuclei E and I through inhibitory pathways "Gpe→E" and "Gpe→I", respectively. Because ${\mathrm{w}}_{\mathrm{G}\mathrm{G}}$ itself is inhibitory, the increase of ${\mathrm{w}}_{\mathrm{G}\mathrm{G}}$ exerts an excitatory effect on E and I. Figure 13A shows the oscillation critical conditions simulated on plane (${\mathrm{w}}_{\mathrm{G}\mathrm{G}}$, $\frac{T}{\tau }$). As shown in the figure, with the increase of ${\mathrm{w}}_{\mathrm{G}\mathrm{G}}$, the state of the system will transition from oscillation to HSS. The ADR of E and I gradually increases with the increase of ${\mathrm{w}}_{\mathrm{G}\mathrm{G}}$, as shown in Figure 13B, C. When ${\mathrm{w}}_{\mathrm{G}\mathrm{G}}$ increases to a certain value, HP2 bifurcation occurs in both E and I, as shown in Figure 13D, E. Figure 13F is the oscillation time series obtained by taking ${\mathrm{w}}_{\mathrm{G}\mathrm{G}}$ = 5. Figure 13G is the stable time series diagram of high discharge rate obtained by taking ${\mathrm{w}}_{\mathrm{G}\mathrm{G}}$ = 50.

In this model, we assume that the cortex has a constant input C to both the STN and E. Next, we consider the influence of the change of C on the oscillations of E and I. From the connection structure in Figure 1, we know that the influence of C acting on E on circuit E-I is direct, and the influence of C acting on STN on circuit E-I is indirect. Therefore, we speculate that the direct effect on E is the main one, which is also consistent with the simulation results in Figure 14. C mainly affects I through the excitatory pathway "E→I". Therefore, the increase of C exerts excitatory effects on both E and I. It can be seen from the critical curve in Figure 14A that the increase of C can make the system have bidirectional Hopf bifurcation. When C is small, the system is in a stable state; when C is large, the system is in HSS state. The ADRs of E and I both increase with the increase of C, as shown in Figure 14B, C. Figure 14D, E show two specific state bifurcation diagrams, and the direction of bidirectional Hopf bifurcation is consistent with Figure 14A. Figure 14F–H show state time series diagrams obtained by taking C = 2, C = 9 and C = 14, respectively, which visually show the discharge state of loop E-I under different parameters.

${\mathrm{w}}_{\mathrm{C}\mathrm{S}}$ is the excitatory coupling weight directly acting on STN by the cortex. Obviously, ${\mathrm{w}}_{\mathrm{C}\mathrm{S}}$ mainly acts on nuclei E and I through the pathways "STN→GPe→E" and "STN→GPe→I", respectively, and these two pathways are generally inhibitory. When ${\mathrm{w}}_{\mathrm{C}\mathrm{S}}$ is relatively small, the discharge rate of loop E-I is relatively high, and E and I are in HSS state. When ${\mathrm{w}}_{\mathrm{C}\mathrm{S}}$ increases to a certain value, HSS will be suppressed, and the system will enter the oscillation state, as shown by the arrow in Figure 15A. To sum up, the ADR of E and I will decrease with the increase of ${\mathrm{w}}_{\mathrm{C}\mathrm{S}}$, as shown in Figure 15B, C. Figure 15D, E show two specific state bifurcation diagrams. As shown in the figures, when ${\mathrm{w}}_{\mathrm{C}\mathrm{S}}$ increases to a certain value, HP2 bifurcation occurs in E and I.

Interestingly, when we simulate some parameters, we observe that HP bifurcation may occur at both ends of the steady state, which is a relatively rare phenomenon in previous model studies. As shown in Figure 16A, when ${\mathrm{w}}_{\mathrm{S}\mathrm{G}}$ increases, the oscillation state of E first enters stability after HP2 bifurcation, and then the stable state enters oscillation state with the emergence of HP1. In Figure 16B, we simulate the ADR of E; we observe that the slope of ADR changes greatly at the bifurcation point, as shown by SLP1 and SLP2 in the figure. With the increase of ${\mathrm{w}}_{\mathrm{S}\mathrm{G}}$, HP1 and HP2 appear in the state of I, as shown in Figure 16C. The slope of ADR of I also changes greatly at the bifurcation point, as shown in Figure 16D. Therefore, we speculate that the sudden change of the slope of ADR is also an important reason for bifurcation transition. A similar phenomenon can also be observed in Figure 17.

6.   Influence of different coupling weights on oscillation amplitude and frequency

6.1.   Oscillation amplitude characteristics of E and I

Oscillation amplitude is an important feature of periodic motion. The change of oscillation amplitude is closely related to the evolution of the disease and Parkinson's symptoms. In this part, we mainly discuss the influence of coupling weight on oscillation amplitude.

Figure 18 describes the influence of excitability-inhibition coupling weights (${w}_{EI}$, ${w}_{IE}$) inside loop E-I on oscillation amplitude. Figure 18A, B show the effect of ${w}_{EI}$ and ${w}_{IE}$ on the amplitude of E in three-dimensional and two-dimensional cases, respectively. As shown in Figure 18B, when ${w}_{EI}$ is fixed at a suitable value, the increase or decrease of ${w}_{IE}$ will make the oscillation amplitude disappear, as shown by the two-way arrow in the figure. Figure 18C shows the variation trend of the amplitude of E with the increase of ${w}_{IE}$, which is obtained by taking ${w}_{EI}$ = 18. When ${w}_{IE}$ is relatively small or too large, the oscillation amplitude will be 0. Within the range of oscillation parameters, the amplitude of E first increases and then decreases, which is consistent with the bidirectional bifurcation in Figure 5D. In the oscillation region, the influence of ${w}_{EI}$ on the amplitude of E is relatively small. Figure 18D shows the change trend of the amplitude of E with the increase of ${w}_{EI}$, which is obtained by taking ${w}_{IE}$ = 5. When ${w}_{EI}$ increases to a certain value, the oscillation amplitude begins to appear and E begins to oscillate, which is consistent with the situation described in bifurcation diagram Figure 4E. Figure 18E, F describe the change trend of the oscillation amplitude of I in the plane (${w}_{EI}$, ${w}_{IE}$). The amplitude of I is much higher than that of E. Although the bifurcation mechanism of oscillation in E and I is different, their amplitude change trend under the influence of parameters ${w}_{EI}$ and ${w}_{IE}$ is similar, as shown in Figure 18F. Figure 18G, H show the amplitude changes of I obtained by taking ${w}_{EI}$ = 18 and ${w}_{IE}$ = 5, respectively, and they are consistent with the bifurcation trends shown in Figures 4C and 5E.

In this model, we assume that there are two external constant inputs: striatal inhibitory input Str and cortical excitatory input C. Figure 19 depicts the effect of these two inputs on the oscillation amplitudes of nuclei E and I. Figure 19A is a three-dimensional stereogram of the amplitude of E, and Figure 19B is a corresponding plane graph. As shown by the two-way arrow in Figure 19B, the increase or decrease of C can make the oscillation amplitude disappear. Figure 19C is a specific amplitude change graph, where we fix Str to 25. With the increase of C, the E experiences SS, OS and HSS in turn, which is consistent with the bifurcation analysis in Figure 14D. As shown by the white arrow in Figure 19B, the increase or decrease of Str can also make the oscillation disappear, and this direction depends on the size of C. For example, we take C = 7 and get a change trend of amplitude with Str in Figure 19D. When Str increases to a certain value, the oscillation begins to appear, and the oscillation amplitude basically remains unchanged, which is consistent with the bifurcation result in Figure 11D. Due to the synchronous resonance of E and I, the characteristics of oscillation amplitude in I are similar to those in E, as shown in Figure 19E–H. The amplitude of I is also greater than that of E, but their gap is smaller than that in Figure 18.

${w}_{GI}$ and ${w}_{GE}$ are two direct inhibitory coupling weights feedback from GPe to cortex. Figure 20 simulates the variation trend of the oscillation amplitude of loop E-I with the increase of ${w}_{GI}$ and ${w}_{GE}$. Figure 20A is a three-dimensional rendering. Figure 20B is a graph of the amplitude change of E obtained in the parameter plane (${w}_{GE}$, ${w}_{GI}$). When both ${w}_{GI}$ and ${w}_{GE}$ are small, the oscillation amplitude of E is relatively large. When ${w}_{GI}$ is fixed (such as 12), the increase of ${w}_{GE}$ will increase the amplitude of E from 0. When ${w}_{GE}$ increases to a certain value, the oscillation amplitude will drop to 0, as shown in Figure 20C; this is consistent with the bidirectional bifurcation case described in Figure 7D. When ${w}_{GE}$ = 10 is fixed, the oscillation amplitude will disappear with the increase of ${w}_{GI}$, as shown in Figure 20D; this is consistent with the bifurcation trend in Figure 8C. Figure 20E, F describe the changing trend of the amplitude of I. Under the influence of parameters ${w}_{GI}$ and ${w}_{GE}$, the amplitude difference between E and I is not very obvious, and even in some parameter ranges, the amplitude of I is lower than that of E. This is mainly because the increase of ${w}_{GI}$ and ${w}_{GE}$ both exert inhibitory effects on I.

6.2.   The frequency of system oscillation

As mentioned in the introduction, excessive beta frequency oscillation activity is a typical feature of Parkinson's patients. Figure 21 describes the oscillation frequency of the system under the influence of four different parameters. In the steady state, the oscillation frequency is 0; in the oscillation state, the frequency of the system falls in the range of 13–30 Hz, which is a typical beta frequency band oscillation. Figure 21 shows that the model we establish is effective for studying Parkinson's beta oscillation.

7.   Conclusions

In this paper, we propose a new improved cortex-basal ganglia model to study the origin and control mechanism of Parkinson beta oscillation. In this model, the cortex contains excitatory and inhibitory nuclei. The basal ganglia mainly contains the subthalamic nucleus and globus pallidus external segment, and the striatum is simplified as a constant input. Different from many previous models, in this new model, we assume that there are two direct inhibitory projections from GPe to the two nuclei of the cortex. This model includes two excitatory-inhibitory coupling loops: E-I and STN-GPe. Theoretically, the circuit with this structure is prone to oscillation under certain conditions.

First, we obtain two nonlinear oscillation critical conditions (50) and (51) of systems (1)–(4) by means of quasi linearization, Laplace transform and eigenvalue analysis. These two formulas show the relationship between the parameters in the model and the conditions to be met when the oscillation occurs. Equation (50) mainly describes the relationship between parameters and time delay in loop E-I, which rarely appears in previous model research, and is a new calculation formula. Moreover, we do not need to assume that the self-feedback connection of GPe is 0 to get Eq (50), which is also different from previous studies. Equation (51) is a calculation formula obtained by assuming that the self-feedback connection of GPe is 0, and it mainly describes the relationship between the coupling weight, time delay and time constant of the loop STN-GPe. We mainly use Eq (50) to simulate the critical conditions of oscillation, combined with bifurcation analysis and calculation of average discharge rate to discuss the generation of oscillation and transition dynamic mechanism. Since the oscillation characteristics of E and I nuclei have received little attention in previous studies, our simulation and analysis mainly focus on these two nuclei.

We first discuss the influence of the two coupling weights ${w}_{EI}$ and ${w}_{IE}$ in the loop E-I and the time delay on the oscillation. The critical condition of oscillation is a monotonic or U-shaped curve. The system above the curve is oscillatory, and the system below the curve is stable. The same parameter has different effects on E and I, and the origin mechanism of oscillation in different nuclei may also be different. For example, the increase of ${w}_{EI}$ promotes the discharge of I but suppresses the discharge of E. Average discharge rate (ADR) describes the discharge capacity of nuclei. Generally, when ADR exceeds a certain value, the nucleus will oscillate, and the state of the nucleus will change from the stability of lower ADR to oscillation. We denote this type of Hopf bifurcation as HP1, as shown in Figure 4C. In some cases, if the ADR of the nucleus is too large, the nucleus will be in the discharge saturation state or high discharge rate state; and in order to produce oscillation, it is necessary to reduce the ADR of the nucleus, as shown in Figure 4E. We call the bifurcation in this case HP2. HP1 and HP2 can explain all bifurcation behaviors in this model. The increase of ${w}_{IE}$ exerts both inhibitory effects on E and I, and at this time, the origin mechanism of oscillation in E and I is similar, as shown by the U-shaped critical curve and bidirectional bifurcation diagram in Figure 5. Time delay is an indispensable parameter in this model, as we know from Eqs (52) and (53) that the system will oscillate only when the time delay exceeds a certain value. In Figure 6, we simulate the state transition diagram of E and I with respect to the increase of delay. As shown in the figure, when the time delay is large enough, HP1 bifurcation will occur in E and I, and the oscillation amplitude also increases with the increase of T.

Then, we study the influence of several projections that directly affect the loop E-I on the oscillation. ${w}_{GE}$ and ${w}_{GI}$ are the two inhibitory feedback coupling weights from GPe to E-I loop. ${w}_{GE}$ exerts an inhibitory effect on E and I, and the increase of ${w}_{GE}$ can make the loop state realize the transfer of HSS-OS-SS, as shown in Figure 7. The effect of ${w}_{GI}$ on E and I is just the opposite. In general, the increase of ${w}_{GI}$ has an exciting effect on E and an inhibitory effect on I. Therefore, the oscillation in E can disappear through HP2 bifurcation, and the oscillation in I can disappear through HP1 bifurcation, as shown in Figure 8C, E. In addition, ${w}_{CE}$ is the excitatory projection from cortical tissue to E, which is also the only excitatory projection in this model that has a direct effect on circuit E-I. Obviously, the increase of ${w}_{CE}$ exerts both excitatory effects on E and I, so the loop realizes the transition from SS to OS and then to HSS, as shown in Figure 9D, E.

In addition to the E-I loop, this model also includes the STN-GPe loop, which is also an inhibitory-excitatory coupling structure. The oscillation in E-I may be spontaneous or induced by the oscillation in the loop STN-GPe, at least affected by the loop STN-GPe. We systematically discuss the influence of the coupling weight and input related to the loop STN-GPe on the oscillation. ${w}_{GS}$ and ${w}_{SG}$ are a pair of coupling weights within the loop STN-GPe. After analysis, ${w}_{GS}$ exerts excitatory influence on E and I, and the oscillation of loop E-I is transferred from HP1 bifurcation, as shown in Figure 10. ${w}_{SG}$ has different effects on E and I, as it exerts inhibitory effects on I and excitatory effects on E. Therefore, when ${w}_{SG}$ increases to a certain value, the oscillation in I disappears due to HP1, and the oscillation in E disappears due to HP2, as shown in Figure 12C, E. Str is an inhibitory striatal constant input acting on GPe, and its mechanism is similar to that of ${w}_{GS}$. The increase of Str exerts an excitatory effect on E and I. Therefore, the oscillation mechanism shown in Figure 11 is similar to that of ${w}_{GS}$.

${w}_{GG}$ is the only self-feedback projection in this model. In many previous theoretical studies ^[20], in order to make the derivation feasible, it is always assumed that ${w}_{GG}$ = 0, and the influence of ${w}_{GG}$ on oscillation is ignored. In this model, the Eq (50) we derived can include parameter ${w}_{GG}$, which is an improvement on the previous research. In Figure 13, we simulate the effect of ${w}_{GG}$ on oscillation. In general, the increase of ${w}_{GG}$ exerts an excitatory effect on E and I. When the E-I circuit is in oscillation, the increase of ${w}_{GG}$ will enhance the discharge capacity of the system and enter the HSS state, as shown in Figure 13D, E.

Theoretically, we also observe the phenomenon of secondary bifurcation at both ends of the steady state, as shown in Figures 16 and 17. This is not involved in the previous theoretical model research. Through analysis, we believe that the sudden change of the slope of ADR is also an important reason for Hopf bifurcation.

Finally, we study the characteristics of beta oscillation amplitude. In general, the amplitude is the smallest near the bifurcation point; as the oscillation evolves, the amplitude will gradually increase. This is consistent with the observation in the experiment; the change of oscillation amplitude is closely related to the evolution of Parkinson's disease. Different parameters have different effects on oscillation amplitude. In the parameter pair (${w}_{EI}$, ${w}_{IE}$), ${w}_{IE}$ is relatively sensitive to the influence of amplitude. When both ${w}_{EI}$ and ${w}_{IE}$ are relatively small, the amplitude of E is relatively large. When ${w}_{EI}$ is relatively large, the amplitude of I is also relatively large. In general, under the influence of parameters ${w}_{EI}$ and ${w}_{IE}$, the amplitude of I is much larger than that of E. In the parameter pair (C, Str), Str has little effect on the amplitude within the range of oscillation parameters; the increase of C will increase the amplitude first and then decrease. Under the influence of parameters (C, Str), the amplitude of I is also relatively larger than that of E, but the amplitude difference between E and I is smaller than that in parameter pair (${w}_{EI}$, ${w}_{IE}$). In the parameter pair (${w}_{GE}$, ${w}_{GI}$), the amplitude difference between E and I is small; when both ${w}_{GE}$ and ${w}_{GI}$ are small, the amplitude is large, as shown in Figure 20. We simulate the influence of several key parameters on the oscillation frequency, and we observe that the frequency of the system basically falls between 13 and 30 Hz.

In the study of oscillation mechanism, we also notice that some control methods can make oscillation disappear. As shown by the two-way arrow in Figure 18, the increase or decrease of ${w}_{IE}$ can make the oscillation activity disappear. Obviously, ${w}_{IE}$ directly and indirectly regulates the discharge capacity of E and I, respectively. Similarly, as shown by the two-way black arrow in Figure 19, the role of C is also to directly or indirectly regulate the discharge capacity of circuit E-I. When we adjust C to a proper value, the oscillation will disappear. ${w}_{GE}$ and ${w}_{GI}$ in Figure 20 are two coupling weights that directly adjust the discharge capacity of E and I, and their increase or decrease can also make the oscillation disappear. One of the main principles of deep brain stimulation technology commonly used in clinical medicine is to alleviate or eliminate the oscillatory activity by regulating the electrophysiological activity of neurons at the stimulation site. Therefore, we speculate that the cortex may be an effective stimulation area to control Parkinson's disease ^[64,65,66], which can provide a basis for the study of Parkinson's oscillation regulation strategy.

In this paper, we mainly study the influence of various key physiological parameters on the oscillation activity in the E-I circuit. In this model, in addition to the factors of the E-I circuit itself, the STN-GPe circuit has obvious input to the E and I. Figure 22 depicts some possible relationships of the oscillation modes between four populations. When the loop STN-GPe is in a stable state, the E may be stable or oscillating, as shown in Figure 22A, B. When the loop STN-GPe is in an oscillating state, the E may oscillate or stabilize, as shown in Figure 22C, D. Therefore, we speculate that some oscillations in loop E-I are spontaneous, and some are caused by the oscillation activity in STN-GPe loop. It is meaningful to accurately distinguish these different oscillation origin mechanisms, which will rely on more complex network models and the introduction of heterogeneity parameters.

In conclusion, in the model proposed in this paper, we systematically study the effects of some main parameters on the oscillations of E and I by mathematical derivation, bifurcation analysis and numerical simulation. The advantage of a mathematical formula is that it can directly describe the relationship between various parameters, can quickly simulate the critical conditions of oscillation and is convenient for theoretical mechanism analysis. The theoretical boundary conditions obtained in this paper are consistent with the trend of bifurcation analysis, which can better explain the mechanism of oscillation origin. The symptoms of Parkinson's patients vary from person to person, and the time delay and coupling strength are two key parameters to describe the interaction between neurons. The results obtained in this paper can help further uncover the pathogenesis of Parkinson's disease and provide guidance for the development of experiments.

Acknowledgments

This research was supported by the National Science Foundation of China (Grant No. 11602092), the Natural Science Foundation of Hubei Province (Grant No. 2018CFB628) and the China Postdoctoral Science Foundation (No. 2018M632184).

Conflict of interest

The authors have no relevant financial or nonfinancial interests to disclose.