# Emergence of Beta Oscillations of a Resonance Model for Parkinson’s Disease

Authors: Yaqian Chen; Junsong Wang; Yanmei Kang; Muhammad Bilal Ghori
Journal: Neural Plasticity
 (2020)

---

## Abstract

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.


---

## Body

## 1. Introduction

Parkinson’s disease ranks as the second most common neurodegenerative disorder after Alzheimer’s disease [1–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.A commonly acknowledged opinion is that these symptoms are associated with the degeneration of dopaminergic neurons in the substantia nigra par compacta, which releases the neurotransmitter dopamine to the basal ganglia [7]. Particularly, excessive beta oscillations (13-30 Hz) have been observed in the local field potential recording from the basal ganglia of Parkinson patients with symptom severity [6–9]. The magnitude of beta oscillations is even related to the severity and degree of bradykinetic or akinetic motor symptoms and rigidity [10, 11]. Some motor symptoms can be ameliorated by suppressing these oscillations by deep brain stimulation (DBS) or dopaminergic medications [12–14]. The recent clinical findings are persuading more attention to deeply explore the mechanisms of beta oscillations associated with Parkinson’s disease.Several computational models have been proposed to elucidate the generation of pathologically exaggerated beta oscillations in Parkinson’s disease [15–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.


## 2. 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-basal–gangalia-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
(1)τSS′=FSwCSEt−TCS−wGSGt−TGS−St,τGG′=FGwSGSt−TSG−wGGGt−TGG−Str−Gt,τEE′=FE−wCCIt−TCC+C−Et,τII′=FIwCCEt−TCC−It,where St, Gt, Et, and It 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 intrinsic excitatory input to cortical excitatory neurons. The parameters Tij and wij 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.Figure 1
(a) Connectivity in the cortical-basal-ganglia-thalamic circuit in the resonance model. This model includes sigmoid functions for STN, GPe, excitatory, and inhibitory populations, which describe the input-output relationship of neurons in the populations as shown in Figure1(b). Here, wij and Tij denote the connection strength and synaptic delay between neural populations i and j, respectively. The arrows represent the excitatory input, and the solid points represent the inhibitory input. (b) Output from the activation function FXinX=S,G,E,I, given by Eq. (2). (c) Derivatives of the activation function from Figure 1(b).
(a)(b)(c)The termsFXX=S,G,E,I are the sigmoid activation functions expressing the relationship between firing rate and synaptic input, shown as follows
(2)FXin=MX1+MX−BX/BXexp−4in/MX,X=S,G,E,I.Here, the constant MX is the maximum firing rate of population X, and BX 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 Table1. It is noteworthy that these values could be estimated directly on the basis of experimental data, more details could be found in Ref. [24].Table 1
The parameter values used in this paper.
ParameterValueParameterValueτS12.8 msτG20 msτE10—20 msτI10—20 msBS10 spk/sBG20 spk/sBE0—20 spk/sBI0—20 spk/sMS300 spk/sMG400 spk/sME50—80 spk/sMI20—330 spk/sStr40.51 spk/sC277.94 spk/sIt 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τXX=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:
(3)τS′=FSwCSEt−T−wGSGt−T−St,τG′=FGwSGSt−T−wGGGt−T−Str−Gt,τE′=FE−wCCIt−T+C−Et,τI′=FIwCCEt−T−It.denoting the equilibrium point of the system asu0=S∗,G∗,E∗,I∗T, and the model (3) around the equilibrium point can be formally rewritten into
(4)dudt=B1ut+B2ut−T+fut−T,with ut=St−S∗,Gt−G∗,Et−E∗,It−I∗T. Here, f=fS,fG,fE,fIT is deduced from the higher-order terms of the vector field and B1=−1/τI4×4(5)B2=0a12a130a21a2200000a3400a430,with
(6)a12=4S∗2wGSτMS2−4S∗wGSτMS,a13=−4S∗2wCSτMS2+4S∗wCSτMS,a21=−4G∗2wSGτMG2+4G∗wSGτMG,a22=4G∗2wGGτMG2−4G∗wGGτMG,a34=4E∗2wCCτME2−4E∗wCCτME,a43=−4I∗2wCCτMI2+4I∗wCCτMI.For the linearized systemdu/dt=B1ut+B2ut−T, the characteristic equation is given by
(7)Δλ=λ+1/τ−a12exp−λT−a13exp−λT0−a21exp−λTλ+1/τ−a22exp−λT0000λ+1/τ−a34exp−λT00−a43exp−λTλ+1/τ=Δ1λΔ2λ=0,where
(8)Δ1λ=λ+1/τ2−λ+1/τa22exp−λT−a12a21exp−2λT,Δ2λ=λ+1τ2−a34a43exp−2λTrepresent 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.As well known, if Eq. (7) has a pair of purely imaginary roots ±iωω>0, the stability should be changed, and a Hopf bifurcation may occurs. So, we consider Δ1iω=0 and Δ2iω=0 next.
(i)
Suppose thatλ=iωω>0 is a root of Δ1λ, i.e.,(9)Δ1λexpλT=−iωa22−a22τ+−ω2+1τ2+i2ωτcosωT+isinωT−a12a21cosωT−isinωT=0.Separating real parts and imaginary parts of Eq. (9), one can obtain
(10)−ω2+1τ2sinωT+2ωτcosωT=ωa22−a12a21sinωT,−ω2+1τ2cosωT−2ωτsinωT=a22τ+a12a21cosωT.Then, one can further obtain(11)sinωT=−a22τ4ω3−τ2a22+τ4a12a21a22ωτ4ω4+2τ2ω2+1−τ4a122a212,(12)cosωT=a22τ3ω2+τa22+τ3a12a21a22τ4ω4+2τ2ω2+1−τ4a122a212.According to Eqs. (11) and (12) and sin2ωT+cos2ωT=1, we have
(13)ω8+k1ω6+k2ω4+k3ω2+k4=0,where
(14)k1=4/τ2−a222,k2=6τ4−3a222τ2−2a12a21a12a21+a222,k3=4τ6−3a222τ4−4a12a21a12a21+a222τ2−a12a21a222,k4=1τ8−a222τ6−2a12a21a12a21+a222τ4−a12a21a222τ2+a12a214.Letz=ω2, then Eq. (13) becomes
(15)z4+k1z3+k2z2+k3z+k4=0.Here, we give the following assumption (H1). Eq. (15) has at least one positive root. Without loss of generality, we assume that it has four positive roots, denoting by zk,k=1,2,⋯,4. Then, Eq. (13) has four positive roots as well, namely, ωk=zk,k=1,2,⋯,4. According to Eq. (12), we can get the corresponding critical value of time delay
(16)Tkj=1ωkarccosa22τ3ωk2+τa22+τ3a12a21a22τ4ωk4+2τ2ωk2+1−τ4a122a212+2jπωk,k=1,2,3,4;j=0,1,2,⋯Then±iωk is a pair of purely imaginary roots of Eq. (7) with Tkj.
(ii)
Suppose thatλ=iωω>0 is a root of Δ2λ, i.e.,(17)ω4+2ω2τ2+1τ4−a34a432=0.So, we know thatω=−a34a43−1/τ2 and the corresponding time delay are
(18)T5j=−12ωarcsin2ωτa34a43+jπω,j=0,1,2,⋯Finally, let the critical time delay asT0=mink=1,2,3,4,5Tkj and the corresponding purely imaginary roots as ±iω0.Next, we need to verify the transversality condition. Differentiating the two sides of Eq. (7) on time delay, we have
(19)dλdT−1=−d1Δ1λ+d2Δ2λd3Δ1λ+d4Δ2λ,where
(20)d1=2λ+1/τ+Ta34a43exp−2λT,d2=2λ+1/τ+a22exp−λTTλ+T/τ−1+2Ta12a21exp−2λT,d3=2λa34a43exp−2λT,d4=λλ+1/τa22exp−λT+2λa12a21exp−2λT.Thus,(21)RedλdTλ=iω0−1=−PRQR+PIQIQR2+QI2,with
(22)PR=−d1IΔ1I+d1RΔ1R−d2IΔ2I+d2RΔ2R,PI=d1IΔ1R+d1RΔ1I+d2IΔ2R+d2RΔ2I,QR=−d3IΔ1I+d3RΔ1R−d4IΔ2I+d4RΔ2R,QI=d3IΔ1R+d3RΔ1I+d4IΔ2R+d4RΔ2I,d1I=2ω0−2Ta34a43sin2ω0T,d1R=2/τ+2Ta34a43cos2ω0T,d2I=2ω0+Ta22ω0cosω0T−a22T/τ−1sinω0T−2Ta12a21sin2ω0T,d2R=2/τ+Ta22ω0sinω0T+a22T/τ−1cosω0T+2Ta12a21cos2ω0T,d3I=2ω0a34a43cos2ω0T,d3R=2ω0a34a43sin2ω0T,d4I=ω0a22τcosω0T+ω02a22sinω0T+2ω0a12a21cos2ω0T,d4R=−ω02a22cosω0T+ω0a22τsinω0T+2ω0a12a21sin2ω0T,Δ1I=2ω0τ−ω0a22cosω0T+a22τsinω0T+a12a21sin2ω0T,Δ1R=−ω02+1τ2−ω0a22sinω0T−a22τcosω0T−a12a21cos2ω0T,Δ2I=2ω0τ+a34a43sin2ω0T,Δ2R=−ω02+1τ2−a34a43cos2ω0T.Therefore, if (H2)PRQR+PIQI≠0 holds, Redλ/dTλ=iω0−1≠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
(23)pλ,e−λτ1,⋯e−λτm=λn+p10λn−1+⋯pn−10λ+pn0+p11λn−1+⋯pn−11λ+pn1e−λτ1+⋯+p1mλn−1+⋯pn−1mλ+pnme−λτm,where τi≥0i=1,⋯,m and pjij=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.ForT=0, Eq. (7) becomes
(24)λ4+d3λ3+d2λ2+d1λ+d0=0,with
(25)d0=1/τ4−a22/τ3−a12a21+a34a43/τ2−a22a34a43/τ+a12a21a34a43,d1=4/τ3−3a22/τ2−2a12a21+a34a43/τ+a22a34a43,d2=4/τ2−3a22/τ−a12a21+a34a43,d3=4/τ−a22.By the Routh–Hurwitz criterion, we know all roots of Eq. (24) have negative real parts if the following condition holds
(26)H3di>0i=0,1,⋯3,d3d2>d1,d3d2d1>d12+d32d0.Through the analysis above, we have if the conditions (H1)-(H3) hold, system (3) undergoes a Hopf bifurcation at the equilibriumu0 when transmission delay T=T0.


### 2.2. Direction and Stability of the Hopf Bifurcation

In Section2.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 T0. At first, let vt=uTt,T=T0+μ,μ∈R, then the system (4) can be transformed into the following form
(27)v̇t=Lμvt+Fμ,vt,where Lμ:C⟶R4,F:R×C⟶R4are given by
(28)Lμvt=T0+μB1vt0+B2vt−1,Fμ,vt=T0+μfvt−1=T0+μf2vt−1+f3vt−1,with
(29)f2vt−1=f2S,f2G,f2E,f2IT=cSwGS2v22t−1+wCS2v32t−1−2wGSwCSv2t−1v3t−1cGwSG2v12t−1+wGG2v22t−1−2wSGwGGv1t−1v2t−1cEwCC2v42t−1cIwCC2v32t−1,(30)f3vt−1=f3S,f3G,f3E,f3IT=eSwCS3v33t−1−wGS3v23t−1+3wCS2wGSv32t−1v2t−1−3wGS2wCSv22t−1v3t−1eGwSG3v13t−1−wGG3v23t−1+3wSG2wGGv12t−1v2t−1−3wGG2wSGv22t−1v1t−1−eEwCC3v43t−1eIwCC3v33t−1,(31)cX=−8X∗MX−X∗2X∗−MXτMX4,(32)eX=32X∗MX−X∗MX2+6X∗2−6MXX∗3τMX6X=S,G,E,I.By the Riesz representation theorem, there exists a functionηθ,μ,θ∈−1,0, such that
(33)Lμφ=∫−10dηθ,μφθ,φ∈C.In fact, we can chooseηθ,μ=T0+μB1δθ+B2δθ+1, where δθ is the Dirac delta function, i.e.,
(34)δθ=1,θ=00,θ≠0.Forφ∈C1−1,0,R4,, we define
(35)Aμφ=dφdθ,−1≤θ<0∫−10dηθ,μφθ,θ=0,(36)Rμφ=0,−1≤θ<0Fμ,φ,θ=0.Then, the system (27) can be transformed into the following operator equation:(37)v̇t=Aμvt+Rμvt.The adjoint operatorA∗ of A is given by
(38)A∗μψ=−dψds,0≤s<1∫−10dηTs,μψ−s,s=0.According to the discussion in Section2.1, we know that ±iω0T0 are eigenvalues of A0 and A∗0, let qθ=1,χ,β,γTeiT0ω0θ−1<θ≤0 be the eigenvectors of A0 corresponding to eigenvalue iω0T0, andq∗s=1/ρ1,χ∗,β∗,γ∗TeiT0ω0s0≤s<1 be the eigenvectors of A∗0 corresponding to the eigenvalue −iω0T0. With a simple computation, we can obtain
(39)χ=a21e−iω0T0iω0+1/τ−a22e−iω0T0,β=1/τ+iω0a13e−iω0T0−a12a21e−iω0T0a131/τ+iω0−a22e−iω0T0,γ=a43a13−a12a21a43e−2iω0T0a131/τ+iω01/τ+iω0−a22e−iω0T0,χ∗=−iω0+1/τa21eiω0T0,β∗=−iω0+1/τa13eiω0T0−iω0+1/τ2−a34a43e2iω0T0,γ∗=a34a13e2iω0T0−iω0+1/τ2−a34a43e2iω0T0.And from the definition of the bilinear inner product(40)ψs,φθ=ψ¯0φ0−∫θ=−10∫ξ=0θψ¯ξ−θdηθ,0φξdξ,we have
(41)ρ¯=1+χχ¯∗+ββ¯∗+γγ¯∗+T0e−iω0T0a21χ¯∗+a12+a22χ¯∗χ+a13+a43γ¯∗β+a34β¯∗γ,such that q∗,q=1, q∗,q¯=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=EC+ES [30, 31]. Here, EC is the two-dimensional subspace spanned by the eigenvectors corresponding to ±iω0, termed the center eigenspace. And ES corresponds to the subspace complementary to EC, in which the real part of all eigenvalues is negative. Then, based on the center manifold theorem, there exist a two-dimensional center manifold C0, and the dynamical flow of the system (24) on it can be transformed into
(42)vtθ=qθzt+q¯θz¯t+Wzt,z¯t,with
(43)zt=q∗θ,vt,Wzt,z¯t=vtθ−ztqθ−z¯tq¯θ=vtθ−2Reztqθ,where z is the local coordinate for C0 in the direction of q for the solution of Eq.(35), satisfying
(44)żt=q∗θ,v̇t=q∗θ,Aμvt+Rμvt=iT0ω0z+q¯∗0F0,Wt,0+2Reztq0,and Wzt,z¯t is the nonlinear map from EC to ES with
(45)Wzt,z¯t=W20θz22+W11θzz¯+W02θz¯22+⋯.Thus, our system in thez,W plane reads
(46)żt=iT0ω0z+gz,z¯Ẇt=AW+Hz,z¯,with
(47)gz,z¯=q∗θ,Rμvt=q¯∗0F0,Wt,0+2Reztq0,Hz,z¯=F0,qz+q¯z¯+W−q∗,F0,qz+q¯z¯+Wq−q¯∗,F0,qz+q¯z¯+Wq¯.Then, we apply the normal form theory to deduce the Poincare normal form for the Hopf bifurcation, i.e.,(48)ξ̇=iω0ξ+c10ξξ2+oξ3,where oξ3 represents all terms of fourth and higher order in ξ, and
(49)c10=i2ω0g20g11−2g112−g0223+g212,where giji+j=2 and g21 can be explicitly determined. Let us rewrite gz,z¯,W as follows:
(50)gz,z¯,W=g20z22+g11zz¯+g02z¯22+g21z2z¯2+oz3,with the term oz3 including all terms of fourth and higher order in z; then, the calculation of giji+j=2 is straightforward.By keepingF0,vt=f2vt−1 and inserting vtθ=qθzt+q¯θz¯t in Eq. (29), we get
(51)g20=∂2g0,0∂z2=T0ρ¯∂2f2S∂z2+χ¯∗∂2f2G∂z2+β¯∗∂2f2E∂z2+γ¯∗∂2f2I∂z2,g02=∂2g0,0∂z¯2=T0ρ¯∂2f2S∂z¯2+χ¯∗∂2f2G∂z¯2+β¯∗∂2f2E∂z¯2+γ¯∗∂2f2I∂z¯2,g11=∂2g0,0∂z∂z¯=T0ρ¯∂2f2S∂z∂z¯+χ¯∗∂2f2G∂z∂z¯+β¯∗∂2f2E∂z∂z¯+γ¯∗∂2f2I∂z∂z¯,where
(52)∂2f2S∂z2=2cSe−2iT0ω0χ2wGS2+β2wCS2−2χβwGSwCS,∂2f2G∂z2=2cGe−2iT0ω0wSG2+χ2wGG2−2χwSGwGG,∂2f2E∂z2=2cEwCC2γ2e−2iT0ω0,∂2f2I∂z2=2cIwCC2β2e−2iT0ω0,∂2f2X∂z¯2=conj∂2f2X∂z2X=S,G,E,I,∂2f2S∂z∂z¯=cS2χχ¯wGS2+2ββ¯wCS2−χβ¯+χ¯βwGSwCS,∂2f2G∂z∂z¯=cG2wSG2+2χχ¯wGG2−χ+χ¯wGSwCS,∂2f2E∂z∂z¯=2cEwCC2γγ¯,∂2f2I∂z∂z¯=2cIwCC2ββ¯.On the other hand, considering(53)g21=∂3g0,0∂z2∂z¯=∂3∂z2∂z¯q∗,f3+∂3∂z2∂z¯q∗,f2,one can calculate g21 by two parts. The calculation of ∂3/∂z2∂z¯q∗,f3 is as straightforward as the calculation of g20. Insertion of vtθ=qθzt+q¯θz¯t into Eq. (30) yields
(54)∂3∂z2∂z¯q∗,f3=T0ρ¯∂3f3S∂z2∂z¯+χ¯∗∂3f3G∂z2∂z¯+β¯∗∂3f3E∂z2∂z¯+γ¯∗∂3f3I∂z2∂z¯,where
(55)∂2f3S∂z2∂z¯=6eSe−iT0ω0β2β¯wCS3−χ2χ¯wGS3+wCS2wGS2ββ¯χ+β2χ¯−wGS2wCS2χχ¯β+χ2β¯,∂2f3G∂z2∂z¯=6eGe−iT0ω0wSG3−χ2χ¯wGG3+wSG2wGG2χ+χ¯−wGG2wSG2χχ¯+χ2,∂2f3E∂z2∂z¯=−6eEe−iT0ω0γ2γ¯wCC3,∂2f3I∂z2∂z¯=6eIe−iT0ω0β2β¯wCC3.The calculation of∂3/∂z2∂z¯q∗,f2 involves the map Wz,z¯ from the subspace EC to its complementary subspace ES. With the second equation of Eq. (46) in mind, Taylor expansion of Hz,z¯ is as follows:
(56)Hz,z¯=H20z22+H11zz¯+H02z¯22+oz2,and then by comparing the corresponding coefficients of Eq. (56) and Eq. (45), we have
(57)A−2iT0ω0W20=−H20,AW11=−H11,A+2iT0ω0W02=−H02,where W20−1 and W11−1 can be computed based on the following equations, respectively.W20θ=ig20/T0ω0q0eiT0ω0θ+ig¯20/3T0ω0q¯0e−iT0ω0θ+E1e2iT0ω0θ,(58)W11θ=−ig11T0ω0q0eiT0ω0θ+ig¯11T0ω0q¯0e−iT0ω0θ+E2,with
(59)E1=2e−2iT0ω02iω0+1/τ−a12e−2iω0T0−a13e−2iω0T00−a21e−2iω0T02iω0+1/τ−a22e−2iω0T000002iω0+1/τ−a34e−2iω0T000−a43e−2iω0T02iω0+1/τ−1cSwGS2χ2+wCS2β2−2wGSwCSχβcGwSG2+wGG2χ2−2wSGwGGχcEwCC2γ2cIwCC2β2,E2=−1/τa12a130a21−1/τ+a220000−1/τa3400a43−1/τ−1cSwGS2χ2+wCS2β2−2wGSwCSχβ¯+χ¯βcGwSG2+wGG2χ2−4wSGwGGReχcEwCC2γ2cIwCC2β2.Then, inserting(60)vtθ=qθzt+q¯θz¯t+W20θ2zt2+W11θztz¯t+W02θ2z¯t2into Eq. (29) to get
(61)∂3∂z2∂z¯q∗,f2=T0ρ¯∂3f2S∂z2∂z¯+χ¯∗∂3f2G∂z2∂z¯+β¯∗∂3f2E∂z2∂z¯+γ¯∗∂3f2I∂z2∂z¯,where
(62)∂2f2S∂z2∂z¯=2cSwGS22χe−iT0ω0W112−1+χ¯eiT0ω0W202−1+wCS22βe−iT0ω0W113−1+β¯eiT0ω0W203−1−2wGSwCSχe−iT0ω0W113−1+βe−iT0ω0W112−1−wGSwCSχ¯eiT0ω0W203−1+β¯eiT0ω0W202−1,∂2f2G∂z2∂z¯=2cGwSG22e−iT0ω0W111−1+eiT0ω0W201−1+wGG22χe−iT0ω0W112−1+χ¯eiT0ω0W202−1−2wSGwGGe−iT0ω0W112−1+χe−iT0ω0W111−1−wSGwGGeiT0ω0W202−1+χ¯eiT0ω0W201−1,∂2f2E∂z2∂z¯=2cEwCC22γe−iT0ω0W114−1+γ¯eiT0ω0W204−1,∂2f2I∂z2∂z¯=2cIwCC22βe−iT0ω0W113−1+β¯eiT0ω0W203−1.Withgiji+j=2 and g21 all available, one can obtain c10 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
(63)με=μ2ε2+Oε3,Tε=2πω01+T2ε2+Oε3,βε=β2ε2+Oε3,with
(64)μ2=−Rec10Reλ′,T2=−Imc10+μ2Imλ′ω0,β2=2Rec10.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 ofT2 determines the period of the bifurcating periodic solutions: if T2>0 (T2<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)


## 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-basal–gangalia-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
(1)τSS′=FSwCSEt−TCS−wGSGt−TGS−St,τGG′=FGwSGSt−TSG−wGGGt−TGG−Str−Gt,τEE′=FE−wCCIt−TCC+C−Et,τII′=FIwCCEt−TCC−It,where St, Gt, Et, and It 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 intrinsic excitatory input to cortical excitatory neurons. The parameters Tij and wij 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.Figure 1
(a) Connectivity in the cortical-basal-ganglia-thalamic circuit in the resonance model. This model includes sigmoid functions for STN, GPe, excitatory, and inhibitory populations, which describe the input-output relationship of neurons in the populations as shown in Figure1(b). Here, wij and Tij denote the connection strength and synaptic delay between neural populations i and j, respectively. The arrows represent the excitatory input, and the solid points represent the inhibitory input. (b) Output from the activation function FXinX=S,G,E,I, given by Eq. (2). (c) Derivatives of the activation function from Figure 1(b).
(a)(b)(c)The termsFXX=S,G,E,I are the sigmoid activation functions expressing the relationship between firing rate and synaptic input, shown as follows
(2)FXin=MX1+MX−BX/BXexp−4in/MX,X=S,G,E,I.Here, the constant MX is the maximum firing rate of population X, and BX 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 Table1. It is noteworthy that these values could be estimated directly on the basis of experimental data, more details could be found in Ref. [24].Table 1
The parameter values used in this paper.
ParameterValueParameterValueτS12.8 msτG20 msτE10—20 msτI10—20 msBS10 spk/sBG20 spk/sBE0—20 spk/sBI0—20 spk/sMS300 spk/sMG400 spk/sME50—80 spk/sMI20—330 spk/sStr40.51 spk/sC277.94 spk/sIt 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τXX=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:
(3)τS′=FSwCSEt−T−wGSGt−T−St,τG′=FGwSGSt−T−wGGGt−T−Str−Gt,τE′=FE−wCCIt−T+C−Et,τI′=FIwCCEt−T−It.denoting the equilibrium point of the system asu0=S∗,G∗,E∗,I∗T, and the model (3) around the equilibrium point can be formally rewritten into
(4)dudt=B1ut+B2ut−T+fut−T,with ut=St−S∗,Gt−G∗,Et−E∗,It−I∗T. Here, f=fS,fG,fE,fIT is deduced from the higher-order terms of the vector field and B1=−1/τI4×4(5)B2=0a12a130a21a2200000a3400a430,with
(6)a12=4S∗2wGSτMS2−4S∗wGSτMS,a13=−4S∗2wCSτMS2+4S∗wCSτMS,a21=−4G∗2wSGτMG2+4G∗wSGτMG,a22=4G∗2wGGτMG2−4G∗wGGτMG,a34=4E∗2wCCτME2−4E∗wCCτME,a43=−4I∗2wCCτMI2+4I∗wCCτMI.For the linearized systemdu/dt=B1ut+B2ut−T, the characteristic equation is given by
(7)Δλ=λ+1/τ−a12exp−λT−a13exp−λT0−a21exp−λTλ+1/τ−a22exp−λT0000λ+1/τ−a34exp−λT00−a43exp−λTλ+1/τ=Δ1λΔ2λ=0,where
(8)Δ1λ=λ+1/τ2−λ+1/τa22exp−λT−a12a21exp−2λT,Δ2λ=λ+1τ2−a34a43exp−2λTrepresent 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.As well known, if Eq. (7) has a pair of purely imaginary roots ±iωω>0, the stability should be changed, and a Hopf bifurcation may occurs. So, we consider Δ1iω=0 and Δ2iω=0 next.
(i)
Suppose thatλ=iωω>0 is a root of Δ1λ, i.e.,(9)Δ1λexpλT=−iωa22−a22τ+−ω2+1τ2+i2ωτcosωT+isinωT−a12a21cosωT−isinωT=0.Separating real parts and imaginary parts of Eq. (9), one can obtain
(10)−ω2+1τ2sinωT+2ωτcosωT=ωa22−a12a21sinωT,−ω2+1τ2cosωT−2ωτsinωT=a22τ+a12a21cosωT.Then, one can further obtain(11)sinωT=−a22τ4ω3−τ2a22+τ4a12a21a22ωτ4ω4+2τ2ω2+1−τ4a122a212,(12)cosωT=a22τ3ω2+τa22+τ3a12a21a22τ4ω4+2τ2ω2+1−τ4a122a212.According to Eqs. (11) and (12) and sin2ωT+cos2ωT=1, we have
(13)ω8+k1ω6+k2ω4+k3ω2+k4=0,where
(14)k1=4/τ2−a222,k2=6τ4−3a222τ2−2a12a21a12a21+a222,k3=4τ6−3a222τ4−4a12a21a12a21+a222τ2−a12a21a222,k4=1τ8−a222τ6−2a12a21a12a21+a222τ4−a12a21a222τ2+a12a214.Letz=ω2, then Eq. (13) becomes
(15)z4+k1z3+k2z2+k3z+k4=0.Here, we give the following assumption (H1). Eq. (15) has at least one positive root. Without loss of generality, we assume that it has four positive roots, denoting by zk,k=1,2,⋯,4. Then, Eq. (13) has four positive roots as well, namely, ωk=zk,k=1,2,⋯,4. According to Eq. (12), we can get the corresponding critical value of time delay
(16)Tkj=1ωkarccosa22τ3ωk2+τa22+τ3a12a21a22τ4ωk4+2τ2ωk2+1−τ4a122a212+2jπωk,k=1,2,3,4;j=0,1,2,⋯Then±iωk is a pair of purely imaginary roots of Eq. (7) with Tkj.
(ii)
Suppose thatλ=iωω>0 is a root of Δ2λ, i.e.,(17)ω4+2ω2τ2+1τ4−a34a432=0.So, we know thatω=−a34a43−1/τ2 and the corresponding time delay are
(18)T5j=−12ωarcsin2ωτa34a43+jπω,j=0,1,2,⋯Finally, let the critical time delay asT0=mink=1,2,3,4,5Tkj and the corresponding purely imaginary roots as ±iω0.Next, we need to verify the transversality condition. Differentiating the two sides of Eq. (7) on time delay, we have
(19)dλdT−1=−d1Δ1λ+d2Δ2λd3Δ1λ+d4Δ2λ,where
(20)d1=2λ+1/τ+Ta34a43exp−2λT,d2=2λ+1/τ+a22exp−λTTλ+T/τ−1+2Ta12a21exp−2λT,d3=2λa34a43exp−2λT,d4=λλ+1/τa22exp−λT+2λa12a21exp−2λT.Thus,(21)RedλdTλ=iω0−1=−PRQR+PIQIQR2+QI2,with
(22)PR=−d1IΔ1I+d1RΔ1R−d2IΔ2I+d2RΔ2R,PI=d1IΔ1R+d1RΔ1I+d2IΔ2R+d2RΔ2I,QR=−d3IΔ1I+d3RΔ1R−d4IΔ2I+d4RΔ2R,QI=d3IΔ1R+d3RΔ1I+d4IΔ2R+d4RΔ2I,d1I=2ω0−2Ta34a43sin2ω0T,d1R=2/τ+2Ta34a43cos2ω0T,d2I=2ω0+Ta22ω0cosω0T−a22T/τ−1sinω0T−2Ta12a21sin2ω0T,d2R=2/τ+Ta22ω0sinω0T+a22T/τ−1cosω0T+2Ta12a21cos2ω0T,d3I=2ω0a34a43cos2ω0T,d3R=2ω0a34a43sin2ω0T,d4I=ω0a22τcosω0T+ω02a22sinω0T+2ω0a12a21cos2ω0T,d4R=−ω02a22cosω0T+ω0a22τsinω0T+2ω0a12a21sin2ω0T,Δ1I=2ω0τ−ω0a22cosω0T+a22τsinω0T+a12a21sin2ω0T,Δ1R=−ω02+1τ2−ω0a22sinω0T−a22τcosω0T−a12a21cos2ω0T,Δ2I=2ω0τ+a34a43sin2ω0T,Δ2R=−ω02+1τ2−a34a43cos2ω0T.Therefore, if (H2)PRQR+PIQI≠0 holds, Redλ/dTλ=iω0−1≠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
(23)pλ,e−λτ1,⋯e−λτm=λn+p10λn−1+⋯pn−10λ+pn0+p11λn−1+⋯pn−11λ+pn1e−λτ1+⋯+p1mλn−1+⋯pn−1mλ+pnme−λτm,where τi≥0i=1,⋯,m and pjij=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.ForT=0, Eq. (7) becomes
(24)λ4+d3λ3+d2λ2+d1λ+d0=0,with
(25)d0=1/τ4−a22/τ3−a12a21+a34a43/τ2−a22a34a43/τ+a12a21a34a43,d1=4/τ3−3a22/τ2−2a12a21+a34a43/τ+a22a34a43,d2=4/τ2−3a22/τ−a12a21+a34a43,d3=4/τ−a22.By the Routh–Hurwitz criterion, we know all roots of Eq. (24) have negative real parts if the following condition holds
(26)H3di>0i=0,1,⋯3,d3d2>d1,d3d2d1>d12+d32d0.Through the analysis above, we have if the conditions (H1)-(H3) hold, system (3) undergoes a Hopf bifurcation at the equilibriumu0 when transmission delay T=T0.


## 2.2. Direction and Stability of the Hopf Bifurcation

In Section2.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 T0. At first, let vt=uTt,T=T0+μ,μ∈R, then the system (4) can be transformed into the following form
(27)v̇t=Lμvt+Fμ,vt,where Lμ:C⟶R4,F:R×C⟶R4are given by
(28)Lμvt=T0+μB1vt0+B2vt−1,Fμ,vt=T0+μfvt−1=T0+μf2vt−1+f3vt−1,with
(29)f2vt−1=f2S,f2G,f2E,f2IT=cSwGS2v22t−1+wCS2v32t−1−2wGSwCSv2t−1v3t−1cGwSG2v12t−1+wGG2v22t−1−2wSGwGGv1t−1v2t−1cEwCC2v42t−1cIwCC2v32t−1,(30)f3vt−1=f3S,f3G,f3E,f3IT=eSwCS3v33t−1−wGS3v23t−1+3wCS2wGSv32t−1v2t−1−3wGS2wCSv22t−1v3t−1eGwSG3v13t−1−wGG3v23t−1+3wSG2wGGv12t−1v2t−1−3wGG2wSGv22t−1v1t−1−eEwCC3v43t−1eIwCC3v33t−1,(31)cX=−8X∗MX−X∗2X∗−MXτMX4,(32)eX=32X∗MX−X∗MX2+6X∗2−6MXX∗3τMX6X=S,G,E,I.By the Riesz representation theorem, there exists a functionηθ,μ,θ∈−1,0, such that
(33)Lμφ=∫−10dηθ,μφθ,φ∈C.In fact, we can chooseηθ,μ=T0+μB1δθ+B2δθ+1, where δθ is the Dirac delta function, i.e.,
(34)δθ=1,θ=00,θ≠0.Forφ∈C1−1,0,R4,, we define
(35)Aμφ=dφdθ,−1≤θ<0∫−10dηθ,μφθ,θ=0,(36)Rμφ=0,−1≤θ<0Fμ,φ,θ=0.Then, the system (27) can be transformed into the following operator equation:(37)v̇t=Aμvt+Rμvt.The adjoint operatorA∗ of A is given by
(38)A∗μψ=−dψds,0≤s<1∫−10dηTs,μψ−s,s=0.According to the discussion in Section2.1, we know that ±iω0T0 are eigenvalues of A0 and A∗0, let qθ=1,χ,β,γTeiT0ω0θ−1<θ≤0 be the eigenvectors of A0 corresponding to eigenvalue iω0T0, andq∗s=1/ρ1,χ∗,β∗,γ∗TeiT0ω0s0≤s<1 be the eigenvectors of A∗0 corresponding to the eigenvalue −iω0T0. With a simple computation, we can obtain
(39)χ=a21e−iω0T0iω0+1/τ−a22e−iω0T0,β=1/τ+iω0a13e−iω0T0−a12a21e−iω0T0a131/τ+iω0−a22e−iω0T0,γ=a43a13−a12a21a43e−2iω0T0a131/τ+iω01/τ+iω0−a22e−iω0T0,χ∗=−iω0+1/τa21eiω0T0,β∗=−iω0+1/τa13eiω0T0−iω0+1/τ2−a34a43e2iω0T0,γ∗=a34a13e2iω0T0−iω0+1/τ2−a34a43e2iω0T0.And from the definition of the bilinear inner product(40)ψs,φθ=ψ¯0φ0−∫θ=−10∫ξ=0θψ¯ξ−θdηθ,0φξdξ,we have
(41)ρ¯=1+χχ¯∗+ββ¯∗+γγ¯∗+T0e−iω0T0a21χ¯∗+a12+a22χ¯∗χ+a13+a43γ¯∗β+a34β¯∗γ,such that q∗,q=1, q∗,q¯=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=EC+ES [30, 31]. Here, EC is the two-dimensional subspace spanned by the eigenvectors corresponding to ±iω0, termed the center eigenspace. And ES corresponds to the subspace complementary to EC, in which the real part of all eigenvalues is negative. Then, based on the center manifold theorem, there exist a two-dimensional center manifold C0, and the dynamical flow of the system (24) on it can be transformed into
(42)vtθ=qθzt+q¯θz¯t+Wzt,z¯t,with
(43)zt=q∗θ,vt,Wzt,z¯t=vtθ−ztqθ−z¯tq¯θ=vtθ−2Reztqθ,where z is the local coordinate for C0 in the direction of q for the solution of Eq.(35), satisfying
(44)żt=q∗θ,v̇t=q∗θ,Aμvt+Rμvt=iT0ω0z+q¯∗0F0,Wt,0+2Reztq0,and Wzt,z¯t is the nonlinear map from EC to ES with
(45)Wzt,z¯t=W20θz22+W11θzz¯+W02θz¯22+⋯.Thus, our system in thez,W plane reads
(46)żt=iT0ω0z+gz,z¯Ẇt=AW+Hz,z¯,with
(47)gz,z¯=q∗θ,Rμvt=q¯∗0F0,Wt,0+2Reztq0,Hz,z¯=F0,qz+q¯z¯+W−q∗,F0,qz+q¯z¯+Wq−q¯∗,F0,qz+q¯z¯+Wq¯.Then, we apply the normal form theory to deduce the Poincare normal form for the Hopf bifurcation, i.e.,(48)ξ̇=iω0ξ+c10ξξ2+oξ3,where oξ3 represents all terms of fourth and higher order in ξ, and
(49)c10=i2ω0g20g11−2g112−g0223+g212,where giji+j=2 and g21 can be explicitly determined. Let us rewrite gz,z¯,W as follows:
(50)gz,z¯,W=g20z22+g11zz¯+g02z¯22+g21z2z¯2+oz3,with the term oz3 including all terms of fourth and higher order in z; then, the calculation of giji+j=2 is straightforward.By keepingF0,vt=f2vt−1 and inserting vtθ=qθzt+q¯θz¯t in Eq. (29), we get
(51)g20=∂2g0,0∂z2=T0ρ¯∂2f2S∂z2+χ¯∗∂2f2G∂z2+β¯∗∂2f2E∂z2+γ¯∗∂2f2I∂z2,g02=∂2g0,0∂z¯2=T0ρ¯∂2f2S∂z¯2+χ¯∗∂2f2G∂z¯2+β¯∗∂2f2E∂z¯2+γ¯∗∂2f2I∂z¯2,g11=∂2g0,0∂z∂z¯=T0ρ¯∂2f2S∂z∂z¯+χ¯∗∂2f2G∂z∂z¯+β¯∗∂2f2E∂z∂z¯+γ¯∗∂2f2I∂z∂z¯,where
(52)∂2f2S∂z2=2cSe−2iT0ω0χ2wGS2+β2wCS2−2χβwGSwCS,∂2f2G∂z2=2cGe−2iT0ω0wSG2+χ2wGG2−2χwSGwGG,∂2f2E∂z2=2cEwCC2γ2e−2iT0ω0,∂2f2I∂z2=2cIwCC2β2e−2iT0ω0,∂2f2X∂z¯2=conj∂2f2X∂z2X=S,G,E,I,∂2f2S∂z∂z¯=cS2χχ¯wGS2+2ββ¯wCS2−χβ¯+χ¯βwGSwCS,∂2f2G∂z∂z¯=cG2wSG2+2χχ¯wGG2−χ+χ¯wGSwCS,∂2f2E∂z∂z¯=2cEwCC2γγ¯,∂2f2I∂z∂z¯=2cIwCC2ββ¯.On the other hand, considering(53)g21=∂3g0,0∂z2∂z¯=∂3∂z2∂z¯q∗,f3+∂3∂z2∂z¯q∗,f2,one can calculate g21 by two parts. The calculation of ∂3/∂z2∂z¯q∗,f3 is as straightforward as the calculation of g20. Insertion of vtθ=qθzt+q¯θz¯t into Eq. (30) yields
(54)∂3∂z2∂z¯q∗,f3=T0ρ¯∂3f3S∂z2∂z¯+χ¯∗∂3f3G∂z2∂z¯+β¯∗∂3f3E∂z2∂z¯+γ¯∗∂3f3I∂z2∂z¯,where
(55)∂2f3S∂z2∂z¯=6eSe−iT0ω0β2β¯wCS3−χ2χ¯wGS3+wCS2wGS2ββ¯χ+β2χ¯−wGS2wCS2χχ¯β+χ2β¯,∂2f3G∂z2∂z¯=6eGe−iT0ω0wSG3−χ2χ¯wGG3+wSG2wGG2χ+χ¯−wGG2wSG2χχ¯+χ2,∂2f3E∂z2∂z¯=−6eEe−iT0ω0γ2γ¯wCC3,∂2f3I∂z2∂z¯=6eIe−iT0ω0β2β¯wCC3.The calculation of∂3/∂z2∂z¯q∗,f2 involves the map Wz,z¯ from the subspace EC to its complementary subspace ES. With the second equation of Eq. (46) in mind, Taylor expansion of Hz,z¯ is as follows:
(56)Hz,z¯=H20z22+H11zz¯+H02z¯22+oz2,and then by comparing the corresponding coefficients of Eq. (56) and Eq. (45), we have
(57)A−2iT0ω0W20=−H20,AW11=−H11,A+2iT0ω0W02=−H02,where W20−1 and W11−1 can be computed based on the following equations, respectively.W20θ=ig20/T0ω0q0eiT0ω0θ+ig¯20/3T0ω0q¯0e−iT0ω0θ+E1e2iT0ω0θ,(58)W11θ=−ig11T0ω0q0eiT0ω0θ+ig¯11T0ω0q¯0e−iT0ω0θ+E2,with
(59)E1=2e−2iT0ω02iω0+1/τ−a12e−2iω0T0−a13e−2iω0T00−a21e−2iω0T02iω0+1/τ−a22e−2iω0T000002iω0+1/τ−a34e−2iω0T000−a43e−2iω0T02iω0+1/τ−1cSwGS2χ2+wCS2β2−2wGSwCSχβcGwSG2+wGG2χ2−2wSGwGGχcEwCC2γ2cIwCC2β2,E2=−1/τa12a130a21−1/τ+a220000−1/τa3400a43−1/τ−1cSwGS2χ2+wCS2β2−2wGSwCSχβ¯+χ¯βcGwSG2+wGG2χ2−4wSGwGGReχcEwCC2γ2cIwCC2β2.Then, inserting(60)vtθ=qθzt+q¯θz¯t+W20θ2zt2+W11θztz¯t+W02θ2z¯t2into Eq. (29) to get
(61)∂3∂z2∂z¯q∗,f2=T0ρ¯∂3f2S∂z2∂z¯+χ¯∗∂3f2G∂z2∂z¯+β¯∗∂3f2E∂z2∂z¯+γ¯∗∂3f2I∂z2∂z¯,where
(62)∂2f2S∂z2∂z¯=2cSwGS22χe−iT0ω0W112−1+χ¯eiT0ω0W202−1+wCS22βe−iT0ω0W113−1+β¯eiT0ω0W203−1−2wGSwCSχe−iT0ω0W113−1+βe−iT0ω0W112−1−wGSwCSχ¯eiT0ω0W203−1+β¯eiT0ω0W202−1,∂2f2G∂z2∂z¯=2cGwSG22e−iT0ω0W111−1+eiT0ω0W201−1+wGG22χe−iT0ω0W112−1+χ¯eiT0ω0W202−1−2wSGwGGe−iT0ω0W112−1+χe−iT0ω0W111−1−wSGwGGeiT0ω0W202−1+χ¯eiT0ω0W201−1,∂2f2E∂z2∂z¯=2cEwCC22γe−iT0ω0W114−1+γ¯eiT0ω0W204−1,∂2f2I∂z2∂z¯=2cIwCC22βe−iT0ω0W113−1+β¯eiT0ω0W203−1.Withgiji+j=2 and g21 all available, one can obtain c10 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
(63)με=μ2ε2+Oε3,Tε=2πω01+T2ε2+Oε3,βε=β2ε2+Oε3,with
(64)μ2=−Rec10Reλ′,T2=−Imc10+μ2Imλ′ω0,β2=2Rec10.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 ofT2 determines the period of the bifurcating periodic solutions: if T2>0 (T2<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)


## 3. 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, wSG=2.56,wGS=3.22,wCC=2.75,wCS=6.60,wGG=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].Figure 2
(a) Bifurcation diagram of the firing rate of the STN populationSt against the transmission delay T. The critical point is consistent with the theoretical result, as indicated by the black dashed line, where the steady state is stable at left and unstable at right. (b) Time series of the firing rate for T=2.5ms. (c) Time series of firing rate for T=8.0ms. (d) The corresponding phase diagram in which the red curve corresponds to Figure 2(b), and the blue curve corresponds to Figure 2(c).
(a)(b)(c)(d)Figure 3
(a) Bifurcation diagrams of the firing rate of the STN populationSt against the connection weight wSG and (b) time series of S and G at critical points for three selected values of wCS ((a1)wCS=4.0, (a2)wCS=6.6, (a3)wCS=10.0). The critical points are consistent with the theoretical results, as indicated by the black dashed lines, where the steady states are stable at left and unstable at right.
(a)(b)Figure 4
(a) Bifurcation diagrams of the firing rate of the STN populationSt against the connection weight wCSand (b) time series of S and G at the left critical points for three selected values of wSG. ((a1)wSG=1.85, (a2)wSG=2.56, (a3)wSG=5.00). The critical points are consistent with the theoretical results, as indicated by the black dashed lines.
(a)(b)


### 3.1. 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 T0=3.6807ms, with μ2=5.5252×10−4, β2=−1.4874×10−5, and T2=8.6223×10−6. Here, the signs of μ2, β2, and T2 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 T0=3.6807ms, 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 Figure2(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.


### 3.2. 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 wij 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 wSG, which controls the strength of the excitatory input from STN to GPe.As shown in Figure3, the model (3) undergoes a supercritical Hopf bifurcation when the connection weight wSG 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 wSG 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 wCS, but the mean firing rates of STN and GPe populations get large as the weight wCS 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.Next, we examine the impact of the connection weightwCS, 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 wCS. It reveals that there are two bifurcation points appear as wCS 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 wCS 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 wSG 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 Figure5(a). It is clear that when the two-dimensional bifurcation parameter wSG,wCS 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 wCS or wSG. 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.Figure 5
(a) Boundary line of Hopf bifurcation in the parameter spacewSG,wCS, where the red points correspond to Hopf bifurcation points in Figure 3, and the blue points correspond to Hopf bifurcation points in Figure 4. (b) Boundary line of Hopf bifurcation in the parameter space wCS,T, where the black curve is the Hopf bifurcation curve, and the equilibrium point is stable in domain I, corresponding to the healthy-type behavior. The unstable region includes three parts, oscillations with an α‐band onset frequency in domain II, a β‐band frequency in domain III, and a γ‐band frequency in domain IV. (c) The evolution of the oscillation frequency via the synaptic delay T with wCS=8.0. (d) Time series of S and G with parameters wCS=8.0 and T=1.0ms in the first panel, T=13.0ms in the second panel, T=4.0ms in the third panel, and T=3.0ms in the last panel, which corresponds to the red points in Figure 5(c), respectively.
(a)(b)(c)(d)The codimension-two bifurcation diagram with the connection weightwCS 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 wCS=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].


## 3.1. 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 T0=3.6807ms, with μ2=5.5252×10−4, β2=−1.4874×10−5, and T2=8.6223×10−6. Here, the signs of μ2, β2, and T2 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 T0=3.6807ms, 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 Figure2(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.


## 3.2. 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 wij 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 wSG, which controls the strength of the excitatory input from STN to GPe.As shown in Figure3, the model (3) undergoes a supercritical Hopf bifurcation when the connection weight wSG 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 wSG 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 wCS, but the mean firing rates of STN and GPe populations get large as the weight wCS 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.Next, we examine the impact of the connection weightwCS, 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 wCS. It reveals that there are two bifurcation points appear as wCS 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 wCS 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 wSG 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 Figure5(a). It is clear that when the two-dimensional bifurcation parameter wSG,wCS 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 wCS or wSG. 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.Figure 5
(a) Boundary line of Hopf bifurcation in the parameter spacewSG,wCS, where the red points correspond to Hopf bifurcation points in Figure 3, and the blue points correspond to Hopf bifurcation points in Figure 4. (b) Boundary line of Hopf bifurcation in the parameter space wCS,T, where the black curve is the Hopf bifurcation curve, and the equilibrium point is stable in domain I, corresponding to the healthy-type behavior. The unstable region includes three parts, oscillations with an α‐band onset frequency in domain II, a β‐band frequency in domain III, and a γ‐band frequency in domain IV. (c) The evolution of the oscillation frequency via the synaptic delay T with wCS=8.0. (d) Time series of S and G with parameters wCS=8.0 and T=1.0ms in the first panel, T=13.0ms in the second panel, T=4.0ms in the third panel, and T=3.0ms in the last panel, which corresponds to the red points in Figure 5(c), respectively.
(a)(b)(c)(d)The codimension-two bifurcation diagram with the connection weightwCS 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 wCS=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].


## 4. 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.



---

*Source: 8824760-2020-11-30.xml*