Emergence of Beta Oscillations of a Resonance Model for Parkinson's Disease

In Parkinson's disease, the excess of beta oscillations in cortical-basal ganglia (BG) circuits has been correlated with normal movement suppression. In this paper, a physiologically based resonance model, generalizing an earlier model of the STN-GPe circuit, is employed to analyze critical dynamics of the occurrence of beta oscillations, which correspond to Hopf bifurcation. With the experimentally measured parameters, conditions for the occurrence of Hopf bifurcation with time delay are deduced by means of linear stability analysis, center manifold theorem, and normal form analysis. It is found that beta oscillations can be induced by increasing synaptic transmission delay. Furthermore, it is revealed that the oscillations originate from interaction among different synaptic connections. Our analytical results are consistent with the previous experimental and simulating findings, thus may provide a more systematic insight into the mechanisms underlying the transient beta bursts.


Introduction
Parkinson's disease ranks as the second most common neurodegenerative disorder after Alzheimer's disease [1][2][3]. According to data from the Parkinson's Foundation [4], more than 10 million people worldwide are suffering from the typical motor symptoms including bradykinesia, muscular rigidity, rest tremor, and postural and gait impairment, as well as nonmotor impairment such as mood and sleep disorders, cognitive decline, urinary symptoms, and incontinence [5,6]. With an incidence ranging from 10-18 per 100000 person years [5], the health of increasingly individuals is debilitated.
Several computational models have been proposed to elucidate the generation of pathologically exaggerated beta oscillations in Parkinson's disease [15][16][17][18][19][20][21]. With focus on the subthalamic nucleus-globus pallidus circuit, Holgado et al. developed a mean-field model and suggested that beta oscillations observed in the basal ganglia arise due to interactions of two nuclei: the subthalamic nucleus (STN) and the globus pallidus pars externa (GPe) [15]. They also found that the Hopf bifurcation can occur when increasing the synaptic weights from healthy to Parkinsonian regimes. Hu et al. [17] used an improved network model consisting of two STN populations and one GP population. Some recent experiments evidenced that in addition to the STN-GP loop, the motor cortex also plays a critical role in the generation of pathological beta oscillations in Parkinson's disease [10,13,22,23]. Pavlides et al. [24] thus proposed a cortex-subthalamic nucleus-globus pallidus model, with the excitatory and inhibitory neuronal populations of the cortex taken into account, and successfully reproduced the beta oscillations in the experimental observation by Tachibana et al. [25]. Since it was proposed under the hypothesis that beta oscillations are generated in the motor cortex, and the basal ganglia resonate to the cortical input, the cortex-subthalamic nucleus-globus pallidus model is also called as resonance model. For brevity, we refer to it as the resonance model from now on.
Note that the resonance model is mainly explored by data fitting and parameter identification [24], thus the fundamental dynamical mechanism by which the beta oscillation is generated has not been disclosed. According to the research with the STN-GP circuit [16], it was found that Hopf bifurcation can induce the pathologically exaggerated beta oscillations. Thus, we wonder whether this is still the dynamical mechanism for the resonance model. We also note that synaptic transmission, the biological process by which a neuron communicates with a target cell across a synapse, always plays a critical part in the neuronal activity [3,26]. For instance, blocking excitatory synaptic transmission could decrease neural firing and reveal spontaneous firing [26]. In particular, experimental evidence suggests that the loss of dopamine in Parkinson's disease could influence dendritic excitability [27]. Recent model investigations [16,17,26] further highlighted that higher synaptic transmission delays and strong synaptic connections between STN and GP populations are beneficial for promoting beta-frequency activity. Therefore, in the present study, we aim to validate the Hopf bifurcation mechanism of the pathological beta oscillation in the resonance model based on the normal form theory, using synaptic transmission delay and synaptic connection weight as bifurcation parameters.

Materials and Methods
2.1. The Model and Existence of Hopf Bifurcation. We begin with a review of the resonance model proposed by Alex Pavlides et al. [24], which investigated the cortico-basalgangalia-thalamic circuit as depicted in Figure 1(a). The model includes two circuits, one composed of interconnected neural populations in STN and GPe, and another composed of excitatory and inhibitory neurons in the cortex. Moreover, the excitatory cortical neurons project excitatory glutamatergic axons to STN. The model is described by a continuum mean-field approach, given by where SðtÞ, GðtÞ, EðtÞ, and IðtÞ are the firing rates of STN, GPe, excitatory, and inhibitory populations, respectively. Str denotes the constant inhibitory input from striatum to GPe, and C denotes the constant component of extrinsic and intrin-sic excitatory input to cortical excitatory neurons. The parameters T ij and w ij represent transmission delay and connection weight, respectively. Here, the subscript i indicates the population from which the signal originates, and the subscript j indicates where the signal is received. τ x denotes the membrane time constants for population x, describing how rapidly the population reacts to its inputs. Notice that the "resonance" is mainly reflected in the hypothesis that the oscillations in basal ganglia resonate to the excitatory cortical input. The terms F X ðX = S, G, E, IÞ are the sigmoid activation functions expressing the relationship between firing rate and synaptic input, shown as follows Here, the constant M X is the maximum firing rate of population X, and B X is the firing rate in the absence of the synaptic input. The curves of these activation functions and their derivatives are shown in Figures 1(b) and 1(c), respectively. The values of all parameters in this paper, except for the transmission delay and connection weights, are summarized in Table 1. It is noteworthy that these values could be estimated directly on the basis of experimental data, more details could be found in Ref. [24].
It is well known that one of the classical mechanism for occurring oscillations is the Hopf bifurcation, in which the attractive limit cycle and asymptotically stable equilibrium point could transform as parameters change. And these two states correspond to the Parkinsonian and healthy state in the STN-GPe model, respectively. Thus, we will investigate the conditions for Hopf bifurcation of the resonance model here. Since the model is difficult to analyze mathematically, we make two simplifications: (i) the membrane time constants τ X ðX = S, G, E, IÞ were taken to have an average value of τ = 10ms; (ii) all transmission delays were taken to be equal, denoted by the single variable T. Then, we could get the following system: denoting the equilibrium point of the system as u 0 = ðS * , G * , E * , I * Þ T , and the model (3) around the equilibrium point can be formally rewritten into with with For the linearized system du/dt = B 1 uðtÞ + B 2 uðt − TÞ, the characteristic equation is given by where represent filtering that occurs when a signal gets through the STN-GPe loop and the cortical excitatory-inhibitory loop, respectively. Thus, the linear instability for Hopf bifurcation can be induced by resonance in filter of the STN-GPe loop or cortical excitatory-inhibitory loop.
Next, we need to verify the transversality condition. Differentiating the two sides of Eq. (7) on time delay, we where Thus, with Therefore, if (H2) P R Q R + P I Q I ≠ 0 holds, Re f ðdλ/dTÞ −1 λ=iω 0 g ≠ 0, the transversality condition for Hopf bifurcation is satisfied.
On the other hand, we need to prove that the remaining roots of Eq. (7) have strictly negative real parts. The following Lemma is used. Lemma 1. [28,29]. Consider the exponential polynomial where τ i ≥ 0ði = 1,⋯,mÞ and p ðiÞ j ðj = 1,⋯,n ; i = 1,⋯,mÞ are constants. As ðτ 1 , τ 2 ,⋯τ m Þ varies, the sum of orders of the zeros of pðλ, e −λτ 1 ,⋯e −λτ m Þ in the open right half plane can change only if a zero appears on or crosses the imaginary axis.
For T = 0, Eq. (7) becomes By the Routh-Hurwitz criterion, we know all roots of Eq. (24) have negative real parts if the following condition holds Through the analysis above, we have if the conditions (H1)-(H3) hold, system (3) undergoes a Hopf bifurcation at the equilibrium u 0 when transmission delay T = T 0 .

Direction and Stability of the Hopf Bifurcation.
In Section 2.1, we have studied the condition for Hopf bifurcation occurring. However, how the system advances towards the parkinsonian state by Hopf bifurcation is not clear. In order to further study the relationship between Hopf bifurcation and the pathological beta oscillation, we turn to the center manifold theorem and normal form method to judge the direction of Hopf bifurcation and the stability of bifurcation periodic solutions at a critical value T 0 . At first, let vðtÞ = uð TtÞ,T = T 0 + μ, μ ∈ R, then the system (4) can be transformed into the following form where with By the Riesz representation theorem, there exists a function ηðθ, μÞ, θ ∈ ½−1, 0, such that In fact, we can choose ηðθ, μÞ = ðT 0 + μÞ½B 1 δðθÞ + B 2 δðθ + 1Þ, where δðθÞ is the Dirac delta function, i.e., For Then, the system (27) can be transformed into the following operator equation: The adjoint operator A * of A is given by According to the discussion in Section 2.1, we know that ±iω 0 T 0 are eigenvalues of Að0Þ and A * ð0Þ, let qðθÞ = ð1, χ, β, γÞ T e iT 0 ω 0 θ ð−1 < θ ≤ 0Þ be the eigenvectors of Að0Þ corresponding to eigenvalue iω 0 T 0 , andq * ðsÞ = ð1/ρÞð1, χ * , β * , γ * Þ T e iT 0 ω 0 s ð0 ≤ s < 1Þ be the eigenvectors of A * ð0Þ corresponding to the eigenvalue −iω 0 T 0 . With a simple computation, we can obtain And from the definition of the bilinear inner product we have such that hq * , qi = 1, hq * , qi = 0.
In the following, we use the center manifold theorem to simplify the system. Note that at the Hopf bifurcation point, the corresponding linear system has a pair of pure imaginary eigenvalues λ = ±iω 0 , and all other eigenvalues have strictly negative real parts. So the whole infinite-dimensional state space C could be decomposed into two complementary subspaces, namely, C = E C + E S [30,31]. Here, E C is the two-dimensional subspace spanned by the eigenvectors corresponding to ±iω 0 , termed the center eigenspace. And E S corresponds to the subspace complementary to E C , in which the real part of all eigenvalues is negative. Then, based on the center manifold theorem, there exist a two-dimensional center manifold C 0, and the dynamical flow of the system (24) on it can be transformed into with where z is the local coordinate for C 0 in the direction of q for the solution of Eq.(35), satisfying and WðzðtÞ, zðtÞÞ is the nonlinear map from E C to E S with Thus, our system in the ðz, WÞ plane reads with Then, we apply the normal form theory to deduce the Poincare normal form for the Hopf bifurcation, i.e., where oðjξj 3 Þ represents all terms of fourth and higher order in jξj, and where g ij ði + j = 2Þ and g 21 can be explicitly determined. Let us rewrite gðz, z, WÞ as follows: with the term oðjzj 3 Þ including all terms of fourth and higher order in jzj; then, the calculation of g ij ði + j = 2Þ is straightforward. By keeping Fð0, v t Þ = f 2 ðv t ð−1ÞÞ and inserting v t ðθÞ = qð θÞzðtÞ + qðθÞ zðtÞ in Eq. (29), we get where On the other hand, considering one can calculate g 21 by two parts. The calculation of ð∂ 3 /∂ z 2 ∂ zÞhq * , f 3 i is as straightforward as the calculation of g 20 . Insertion of v t ðθÞ = qðθÞzðtÞ + qðθÞ zðtÞ into Eq. (30) yields where The calculation of ð∂ 3 /∂z 2 ∂ zÞhq * , f 2 i involves the map Wðz, zÞ from the subspace E C to its complementary subspace E S . With the second equation of Eq. (46) in mind, Taylor expansion of Hðz, zÞ is as follows: and then by comparing the corresponding coefficients of Eq. (56) and Eq. (45), we have where W 20 ð−1Þ and W 11 ð−1Þ can be computed based on the following equations, respectively.
With g ij ði + j = 2Þ and g 21 all available, one can obtain c 1 ð0Þ as Eq. (49), and the model (3) could be transformed into the Poincare form. Based on the results in Ref [30], the coefficient μðεÞ which determine the exists of the periodic solution, the period TðεÞ and the nontrivial Floquet exponent near 0 of the periodic solution are given by with Therefore, we have the following results: (1) The sign of μ 2 determines the direction of the Hopf bifurcation: if μ 2 > 0 (μ 2 < 0), the Hopf bifurcation is supercritical (subcritical) (2) The sign of T 2 determines the period of the bifurcating periodic solutions: if T 2 > 0 (T 2 < 0), the period increases (decreases) (3) The sign of β 2 determines the stability of the bifurcating periodic solutions: if β 2 < 0 (β 2 > 0), the bifurcation periodic solutions are stable (unstable)

Results
Synaptopathy is the earliest step in the Parkinson's disease cascade [32]. To clarify the effects of synaptic transmission delay and connection strength on the onset of beta oscillations, numerical simulations are implemented with model (3) to validate the theoretical results in the previous deduction, with comparison of analytical predictions of Hopf bifurcation types with the numerically calculated bifurcation diagrams that are presented in Figures 2 -4. For convenience, we take the parameters from Ref. [24], that is, w SG = 2:56, w GS = 3:22, w CC = 2:75, w CS = 6:60, w GG = 0:90, C = 277:94 and Str = 40:51: In the calculations, we distinguish between the supercritical and subcritical Hopf bifurcation by ramping up and ramping down the parameter across the critical point with the Euler integration, respectively [33].

The Dependence of Oscillations on Synaptic Transmission
Delay. Experimental data suggested that synaptic transmission delay exists between different neuronal populations [34]. And the time delay could often be a source of oscillation. Therefore, we firstly explore the dependence of oscillation onset on the transmission delay in this resonance model. According to the analytic deduction in Section 2, we know that Hopf bifurcation occurs at T 0 = 3:6807ms, with μ 2 = 5:5252 × 10 −4 , β 2 = −1:4874 × 10 −5 , and T 2 = 8:6223 × 10 −6 . Here, the signs of μ 2 , β 2 , and T 2 imply that the Hopf bifurcation is supercritical, the periodic solution is stable, and the period increases with increasing time delay, respectively. To confirm the theoretical results, let us resort to Figure 2. From the bifurcation diagram (fixed points or local minimums and maximums for each parameter point) against transmission delay in Figure 2(a), we see that the model (3) undergoes a supercritical Hopf bifurcation at T 0 = 3:6807 ms, starting from a stable equilibrium point (Figure 2(b)) and entering into a stable limit loop (Figure 2(c)). The two kinds of attractors are usually referred to as a healthy state and an oscillation state, which falls into the beta band here [16]. Moreover, as the time delay increases, the period gradually becomes large. The simulated results are totally consistent with the theoretical analysis, and this enables us to do cross validation throughout Section 3, although we skip the discussion. Besides, from Figure 2(c), one clearly sees that beta oscillations occur in the excitatory neuronal population of cortex first, and then the basal ganglia resonates at the same frequency. This explains why the model (3) was called as the "resonance" model [24]. Hence, synaptic transmission delay between neuronal populations in the basal ganglia and cortex needs to be sufficiently long to allow them to "charge" and increase their firing rate, which is just coincident with the results in Refs [12,13,17]. Figure 2(d) exhibits the corresponding phase portraits, in which the red curve converges to a point, representing the health state, while the blue curve converges to a limit cycle, corresponding to the Parkinsonian state.

The Dependence of Oscillations on Connection Weights.
Although most of the parameters in the resonance model could be fixed based on the experimental data in Ref. [25], the synaptic connection weights w ij cannot be estimated from experimental studies directly, which are usually estimated by fitting the model to experimental recordings. Moreover, it is well known that the depletion of dopamine in basal ganglia maybe related to the synaptic connection strength [35]. As a result, it should be interesting to explore the effect of the synaptic connection strength on the generation of pathological beta oscillation. For this purpose, let us fix the transmission delay T = 3:6708ms and show in Figure 3 the bifurcation diagrams of the firing rate against w SG , which controls the strength of the excitatory input from STN to GPe. Figure 3, the model (3) undergoes a supercritical Hopf bifurcation when the connection weight w SG is changed from the healthy to the Parkinsonian parameter. A beta oscillation of small amplitude appears after the destabilization of the steady state, and the oscillation amplitude increases as the weight w SG enlarges. Thus, pathological beta oscillation could be significantly attenuated or restrained in the resonance model by blocking the STN-GPe connection, which in agreement with the experimental observation reported in Ref. [25]. In addition, as seen from Figures 3(a1)-(a3), the Hopf bifurcation value becomes smaller with increasing w CS , but the mean firing rates of STN and GPe populations get large as the weight w CS enlarges, as shown in Figures 3(b1)-(b3). It reveals that the excitatory input from the cortex to the subthalamic nucleus also affects the neuronal activity in basal ganglia.

As shown in
Next, we examine the impact of the connection weight w CS , which links the STN-GPe circuit and cortical circuit, on the beta oscillation onset. Figures 4(a1)-(a3) exhibit the range of the firing rate of the STN population after the initial transient response as a function of the connection weight w CS . It reveals that there are two bifurcation points appear as w CS increases linearly, consisting with the theoretical prediction. When the excitatory input from the cortex to STN is too strong or too weak, the firing rates of STN and GPe converge to stable equilibrium points. While the beta oscillations are generated when the connection w CS lies in an moderate region, the amplitude firstly increases and then declines until zero. That may be the reason why blocking the excitatory input to STN could abolish beta oscillation in STN [25]. At the same time, it further demonstrates that the parameter range for beta oscillations gets bigger when the connection w SG increases. Figures 4(b1)-(b3)) depict the corresponding time series of firing rates of STN and GPe at the left Hopf bifurcation points, respectively.
3.3. Codimension Two Bifurcation Analysis. The above numerical results show that the critical points for Hopf bifurcation in the resonance model can be influenced not only by STN-GPe connection but also the excitatory input from the cortex to STN. To illustrate how the beta oscillation onset depends on both the STN-GPe and cortex-STN connection strength, we depict the codimension-two bifurcation diagram in Figure 5(a). It is clear that when the two-dimensional bifurcation parameter ðw SG , w CS Þ is located in domain I, the system converges to a stable equilibrium point, but excessive oscillations at beta frequencies occur when the parameter is within domain II. Hence, the system always stays at a healthy state for small w CS or w SG . This is equally to say that the beta oscillations cannot be generated in the single STN-GPe circuit but can originate from interaction among different neuronal population circuits, which is in agreement with Ref [24]. Therefore, blocking the connection between the STN-GPe circuit and cortical circuit may restrain the appearance of beta oscillation.
The codimension-two bifurcation diagram with the connection weight w CS and the transmission delay T is shown in Figure 5(b), where the parameter region is separated by the Hopf bifurcation curve into parts: one is related with steady firing rate and the other is about the oscillatory firing rate. What is more, the oscillatory region can be divided into three parts: alpha oscillation (marked as II with 8-13 Hz oscillation frequencies), beta oscillation (III, 13-30 Hz), and gamma oscillation (IV, larger than 30 Hz). With w CS = 8:0 fixed, the variation of the nonzero oscillation frequency with transmission delay in Figure 5(c) exhibits the transition among oscillations at gamma, beta, and alpha band frequency in turn, and the firing rate directly evolves from a steady value into the gamma oscillation, but the beta band frequency occupies a more wide region. For intuitiveness, the firing rates of STN and GPe with different band frequencies are exemplified in Figure 5(d). From Figures 5(b) and 5(c), it is clear that a moderate range of synaptic delay is responsible for the emergence of beta oscillation, which contains the model parameters in Refs. [15,24].

Discussion
We have investigated critical conditions for pathological beta oscillation onset in the resonance model based on the center manifold theorem and normal form analysis. It is confirmed that the model undergoes a supercritical Hopf bifurcation as the synaptic transmission delay increases, which governs the transitions from the healthy state to the Parkinsonian state. It is found that a strong excitatory connection from STN to GPe is favorable for the generation of beta oscillations, while excessive excitatory input from cortex to STN would suppress beta oscillations. Particularly, the codimension-two bifurcation diagram suggests that the beta oscillation onset depends on the interaction of the STN-GPe circuit and Cortex-STN synaptic connection. Our investigation has demonstrated that a suitable transmission delay is responsible for the emergence of the beta oscillation. The investigation could be inspiring for clinical physician in treating Parkinsonian patients.
In the near future, this study can be extended to generalized models with more biological conditions such as the feedback connection from the STN-GPe circuit to the cortex. Note that the model of this study considers only four populations, namely, the excitatory population and the inhibitory population of cortex, the subthalamic nucleus and globus pallidus external segment, thus, for more insight in this regard, one may consider the more complicated model [21] which can take striatum and globus pallidus internal segment into account as well. In addition, the effect of the synaptic plasticity and the environmental fluctuations on the onset of beta oscillation should also be worthy to be explored.

Data Availability
The data of parameters used to support the findings of this study are included within the article.

Conflicts of Interest
The authors declare no competing financial interests.