Dynamical mechanism of parkinsonian beta oscillation in a heterogenous subthalamopallidal network

Dysfunction of basal ganglia is associated with the pathogenesis of Parkinson's disease including alteration of firing rate and excessive beta-band (13-30 Hz) synchronization activity. Neuronal heterogeneity enriches dynamics of external globus pallidus (GPe), especially showing significant differences in firing alterations under pathological state. The precise mechanism underlying these neural signatures remains elusive. To address this, we propose a subthalamopallidal network containing two classes of GPe neurons, calcium-binding protein parvalbumin (PV) and Lim homeobox (Lhx6) GPe. Our results show that Lhx6 GPe tends to rein in synchronous behavior and abnormal activity of PV GPe. Under the pathological condition, the alteration of synaptic coupling in a heterogenous pallidal network manifests itself as a direct increase of inhibitory input to PV GPe or an indirect elevation of Lhx6 GPe firing rate. These essentially enhance the inhibition of PV GPe, which results in beta-band synchronous bursting. The subthalamic nucleus (STN) is instrumental in stabilizing the spiking sequence of GPe neurons, inhibiting abnormal synchronous oscillations both in control and pathological conditions. In a dopamine-depleted state, the PV GPe-PV GPe pathway notably impacts the enhancement of beta rhythmic oscillations in STN-GPe circuit. Besides, the synaptic coupling in heterogenous pallidal and STN-GPe affect the propagation of abnormal rhythms in pallidal and subthalamopallidal networks, respectively. Our study highlights the pivotal role played by PV GPe in producing and amplifying pathological oscillatory behavior. Our results suggest that STN prevents abnormal oscillatory rhythm in GPe, providing a novel insight into the precise mechanism by which STN affects pallidal activity.


Introduction
Degeneration of dopaminergic neurons in the substantia nigra pars compacta (SNc) is a major trigger for Parkinson's disease (PD) [1], which leads to persistent alterations in local firing rates [2] and global oscillatory behavior among basal ganglia (BG) nuclei [3,4]. The imbalance of direct and indirect pathway firing activity has been confirmed in rodents [5,6] and subsequent optogenetic studies [7]. The classical computational model based on a direct/indirect pathway theory is established to reveal the intrinsic mechanism of pathological changes in firing rates [2]. Exaggerated synchronized oscillatory activities in the beta (b) frequency range (13)(14)(15)(16)(17)(18)(19)(20)(21)(22)(23)(24)(25)(26)(27)(28)(29)(30) across BG nuclei are a well-established phenomenon in experimental parkinsonism [8] and PD patients [9][10][11]. Deep brain electrode recordings of STN and GPi from PD patients demonstrate high coherence with electroencephalogram (EEG) in the beta band [9]. And the depth recordings of STN from PD patients have demonstrated exaggerated local field potential (LFP) bursts of activity at frequencies in the low beta band [10,11]. This abnormal dynamical behavior is correlated with motor impairment [12], but the potential mechanism remains unknown. From the perspective of the computational model, there are multiple microcircuits in the BG network that could manufacture the beta-band oscillation which incorporates inhibitory feedback. It is generally thought that STN-GPe circuit orchestrates a key role in b oscillatory activity [13][14][15]. The physiological experiments reveal richer synaptic connections in BG [16,17], hence the BG circuit is much more complex than the previous model describes. Therefore, the generators have also been formed by cortex [5], striatal inhibitory microcircuit [18], and pallidostriatal loop [19]. It is a critical approach for better understanding the pathogenesis to explore how the firing rate and oscillatory synchronization interact and coordinate.
Recent identification of cell-type diversity in BG enables the emergence of intricate neural circuits [20,21], which challenges the classical model and the proposed hypothesis of the origin of beta-oscillation. The heterogeneity of GPe neurons is identified by firing patterns in the electrophysiological experiment [21][22][23][24], and it has revealed at least two distinct subtypes of pallidal neurons in immunocytochemical studies [16,25]. Genetically, distinct cell populations are defined as calcium-binding protein parvalbumin (PV) and Lim homeobox (Lhx6) in GPe [26]. PV GPe neuron which exhibits the same as the classic GPe neuron accumulates dorsolaterally, whereas Lhx6 GPe neuron which has lower firing rates gather on the ventral medial side [26][27][28][29]. At any moment, there will be significant differences in firing rate among heterogeneous GPe neurons [21], indicating the different roles in aberrant beta-band activity. According to the classical direct/indirect theory, GPe is only considered a relay station in the indirect pathway, which receives projections from striatum and STN and inhibits STN, internal globus pallidus (GPi), and substantia nigra pars reticulata (SNr) [1,2]. While the synaptic connections between STN and the heterogeneous GPe are different [24]: PV GPe inhibits all STN, whereas Lhx6 GPe inhibits the peripheral STN [26,30]; STN gives rise to stronger inputs to PV GPe than Lhx6 GPe, and STN preferentially activates PV GPe neurons [31,32]. Previous studies have found that Lhx6 GPe also receives afferent projections from cortex [29,33] and has extensive projections on SNc [26][27][28], striatum [16,30,34], and others [28,29]. Accounting for these observations, GPe coordinates and integrates information from more brain regions than the indirect pathway itself [35,36]. According to recent physiological experiments, cell-specific pallidal interventions can restore the divergences in firing rates of distinct GPe, which attenuates the pathological activity of BG [37,38]. It indicates that heterogeneous GPe may be of major significance in producing and modulating anomalous oscillations, which cannot be captured by the classical model of GPe simply modeled as a homogeneous nucleus.
STN has been a proven effective therapeutic target for deep brain stimulation (DBS) for PD [39][40][41], indicating the significance of STN in electrophysiological abnormalities. Since the origin of the pathological oscillations cannot be determined yet, the mechanism by which STN-DBS treatment is effective and ineffective is also elusive to dissect. It is generally thought that STN-GPe circuit orchestrates a key role in abnormal BG dynamics, which is identified as the generator and amplifier [14,15,42]. STN receives enhanced inhibition from GPe to generate beta oscillations, and are repeatedly enhanced within the STN-GPe network [14]. Based on the bifurcation analysis of the neuronal model Ca dynamics, GPe exhibit intermittent phase-locking with the input from STN in the beta band [15]. Recent experiments have found that inhibition of the STN has slight effects on cortical b oscillations [43] and exacerbates abnormal firing in GPe [44], suggesting that the STN may not be the origin of the abnormal rhythm. And STN is found to have a stabilizing effect on pallidal spike sequences and achieves stable movements [45]. These experiments indicate that the STN may function to halt abnormal rhythms, challenging the hypothesis that the STN-GPe is the origin of abnormal oscillations.
The microcircuit of heterogeneous GPe is shown to independently generate b oscillations [46,47]. However, studies of computational models do not integrate heterogeneous GPe well, and in particular, the role of heterogeneous GPe in STN-GPe loops is not clear. To investigate the conditions under which beta oscillation is modulated in STN-GPe circuit, we develop a biophysical model of the subthalamopallidal network by adapting previous GPe cellular models [48]. Our model incorporates Lhx6 GPe units and their sparse and structural connectivity into STN-GPe circuit ( Fig. 1). Through studying, we find that heterogeneous GPe neurons play different roles in aberrant beta-band activity. Increased inhibitory input to PV GPe leads to synchronous bursts in the beta band, while Lhx6 GPe tends to suppress synchronous behavior and abnormal activity. Besides, STN seems to play a different role in stabilizing the spiking sequence of GPe neurons and inhibiting abnormal synchronous oscillations than previous models. These results highlight the roles of the heterogeneous GPe in modulating beta oscillations and STN in affecting pallidal activity.
The present paper is organized as follows. First, we present the details of the heterogenous subthalamopallidal network and individual neuron dynamics.
Then, we characterize the effects of heterogenous GPe on b synchronous activities in the constrained network to reveal the critical role of PV GPe. We analyze the oscillations in the STN-GPe circuit and explore the effects of neuronal heterogeneity. Based on the pathological pallidal networks, we investigate how STN affects aberrant oscillatory behavior. We simulate the propagation of b oscillations in terms of leading in external oscillatory input to subthalamopallidal network. Finally, we discuss our results and limitations.

Network model
The subthalamopallidal model involves STN and two subtypes of GPe, i.e., PV GPe and Lhx6 GPe neurons in Fig. 1A. For simplicity, we assume that each nucleus is modeled as the same size, which contains only ten neurons. And the network of each nucleus is a periodic structure with the first neuron adjacent to the tenth. Neurons in the subthalamopallidal network are connected in a structured sparse mode (Fig. 1B), similar to previous studies [46,48].
Modifications are made such that the numbers of projections between STN and heterogeneous GPe are equal. Each STN neuron sends excitatory projection to Fig. 1 A Connectivity schematic diagram of the subthalamopallidal network model. The model is mainly composed of two nuclei: STN and GPe in the blue dotted box. GPe is constituted by calcium-binding protein parvalbumin globus pallidus pars external (PV GPe) and Lim homeobox 6 globus pallidus pars external (Lhx6 GPe). Every black line with circle represents an inhibitory current, and each red arrow indicates an excitatory current. B Sparse connections within the subthalamopallidal model. The network is composed of STN, PV GPe and Lhx6 GPe with 10 cells in each neuronal population. Each STN neuron sends excitatory projection to two adjacent PV GPe and two closest Lhx6 GPe neurons. Each PV GPe cell projects inhibition on two adjacent STN, two Lhx6 GPe and two PV GPe neurons. And each Lhx6 GPe cell inhibits two adjacent STN, two PV GPe, and two Lhx6 GPe cells, respectively. (Color figure online) two adjacent PV GPe and two Lhx6 GPe neurons; and each PV GPe and Lhx6 GPe cell projects inhibition on two neighboring STN cells (Fig. 1B). Although the synaptic connections between STN and the heterogeneous GPe are different [26,30]. The model simplifies the connection between STN and the two types of GPe, and the synaptic differences are mainly reflected in the synaptic strength.
The heterogeneous GPe network is simplified to have an equal number of projections between PV GPe and Lhx6 GPe neurons, which is discussed below. For connections between heterogeneous GPe neurons, each PV GPe neuron projects inhibition on two Lhx6 GPe neurons, the closest neuron and a neuron skipping two located nearest to it. And each Lhx6 GPe cell inhibits two PV GPe cells, skipping the two located nearest to it (Fig. 1B). For the connections between homogeneous GPe neurons, each GPe cell also sends an inhibited signal to its one immediate GPe neighbor and a GPe neuron which skips the located nearest to it shown in Fig. 1B.

Individual neuron dynamics
All neurons in the computational model are described with the single-compartment conductance-based Hodgkin-Huxley (HH) models developed by Terman et al. [14] and So et al. [48]. The membrane potential V i with i 2 STN; PV GPe; Lhx6 GPe f g of each cell is specified as following: C m dV i dt ¼ ÀI Na À I K À I L À I T À I Ca À I AHP À I syn þ I app þ I ex : ð1Þ The membrane capacitance C m is scaling to 1 Á lF=cm 2 for every neuron. The membrane currents of STN and heterogeneous GPe neurons contain the fast sodium current I Na , the fast potassium current I K , the leak current I L are described by equations: The high-threshold Ca 2þ current I Ca and the lowthreshold T-type Ca 2þ current I T for STN are: whereas the currents in PV GPe and Lhx6 GPe take the simpler forms: where g X and E X with X 2 L; Na; K; Ca; T f gare the maximal conductance and reversal potential of each current, respectively. The dynamics of each slowly operating gating variable evolves according to the kinetic equation: where / x is the scaling time constant of the variable x.
The voltage-dependent time constant of the variable x can be written as follows: x 2 h; n; c; r f g ð6Þ where h s x denotes the voltage at which the time constant is midway between its maximum and minimum values, and r s x represents the slope factor for the voltage dependence of the time constant. For all gating variables x ¼ m, h, n, c, a, s, r, the steady-state function x 1 is determined as: Þ ; x 2 m; h; n; c; a; s; r f g ð7Þ where h x is the half activation/inactivation voltage and r x is the slope factor. The T-type current for STN inactivation variable b is described as follows: The Ca 2þ -activated after-hyperpolarization K þ current I AHP is defined as follows: where g AHP is the maximal conductance and k 1 is the dissociation constant of the I AHP . Ca ½ is the intracellular concentration of Ca 2þ , which is governed by the following equation: where e is a constant that describes the effects of buffers, cell volume and the molar charge of calcium, and k Ca denotes the calcium pump rate constant. The applied current I app in the model is a constant background excitation term that represents the sum of all exogenous synaptic inputs from other brain regions. It is applied to the STN and heterogeneous GPe neurons to maintain baseline firing rates of 60 Hz for PV GPe, 18 Hz for Lhx6 GPe, and 13 Hz for STN [22,23].
The numerical values of parameters of neuron models are given in Tables 1, 2. All parameters are obtained from previous studies, which measured directly in physiological experiments or by fitting physiological experimental data [14,48]. Except for the parameters for Lhx6 GPe are discussed in the following.

Synaptic dynamics
Synaptic currents I syn in the model are given by the following equation: where g syn is the corresponding synaptic strength, V j is the membrane potential of the postsynaptic neuron and E syn is the corresponding synaptic reversal potential. The gating variable S for the synaptic transmission is the sum: where N is the number of presynaptic neurons and S k!j describes the kinetics of the gating variable,for each pair of presynaptic neuron k and postsynaptic neuron j. S k!j is a function of the presynaptic voltage V k , and its dynamics depend on the dynamics of the presynaptic neuron k.
For STN efferent, it contains the excitatory projection to PV GPe and Lhx6 GPe, i.e., I snpv and I snlh . And it used the second order alpha synapse: where u t ð Þ ¼ 1, if the presynaptic cell k crosses a threshold of -10 mV at time t; otherwise u t ð Þ ¼ 0. For PV GPe efferent, the synaptic currents include I pvpv , I pvlh , and I pvsn .Similarly, for Lhx6 GPe, it includes I lhpv , I lhlh , and I lhsn . The synaptic variable S k!j in heterogeneous GPe units obeys the following first-order kinetics: where V k is the membrane potential of the presynaptic neuron, a and b are the opening and closing rates of channels, respectively, and H 1 is given by: The specific parameter for synapses dynamics involved is shown in Tables 3 and 4. All parameters are obtained from previous studies, which measured directly in physiological experiments or by fitting physiological experimental data [14,46,48].
External excitatory oscillatory input (e.g., from cortex) and inhibitory input (e.g., from striatum) are modeled as an ungated channel with current governed by where V is the postsynaptic voltage and V cat is a cation reversal potential of 0 mV and À85 mV corresponding to excitatory and inhibitory inputs, respectively. g ex is governed by following: Table 1 Parameters used for maximal conductance (g X in mSÁcm À2 ), neuronal and synaptic reversal potentials (V X in mV), calcium dynamics and applied current (I app in lA Á cm À2 ) for STN, PV GPe and Lhx6 GPe neuron models where f is the frequency of the oscillatory input, and r is a random variable subject to the normal distribution which simulates random channel fluctuations. The oscillation index is defined as the proportion of the area under the spectral power analysis curve of the membrane potential of PV GPe neuron in the b frequency range (i.e., between 13 to 30 Hz) to the area for the whole frequency range (i.e., from 1 to 300 Hz). All oscillation indexes are the average oscillation indexes of ten neurons. The frequency corresponding to the maximum value in the spectral power analysis is considered as the dominant frequency of firing activity.
The LFP signal is a proven index to quantify clinical symptoms [49], and it is also chosen as the indicator for feedback control [50,51]. LFP is widely modeled by transmembrane currents, which are the sum of membrane capacitance current and ionic currents [52][53][54][55]. It is proved that the fixed linear combination of the synaptic currents provided an accurate LFP proxy in the point-like leaky integrate-and-fire (LIF) network models [56], and also confirmed in the estimation of EEG with LFP proxies [57]. Here, we ignore the location of the electrodes and the spatial structure of the network and use a fixed linear combination with coefficients all 1 to approximate the LFP. Therefore, the LFP proxy is modeled as the linear combination of membrane capacitance current and ionic currents as follows:

Parameters for heterogeneous GPe
Significant heterogeneity in the intrinsic electrophysiological properties of GPe neurons is found [21], but the distribution and composition of ion channels in heterogeneous GPe neurons are almost indistinguishable [58]. Whereas most current models of heterogeneous GPe are similar to previous models of GPi and GPe, distinguished only by different external inputs, Lhx6 GPe receives smaller external inputs to get a lower firing rate [46]. However, recent experiments have found that Lhx6 GPe has a richer external excitatory input like the connections of cortical-GPe  compared to PV GPe [29]. Therefore, the difference in heterogeneous GPe firing rates reproduced by differences in external inputs may not be physiologically consistent.
Recent experiments revealed that heterogeneous GPe neurons differed significantly in both firing rate and action potential (AP) width [22,26]. While previous studies considered GPe as a homogeneous nucleus of a prototype GPe (PV GPe), and the model parameters used previously do not capture these differences [48]. This heterogeneity may be reflected not only in differences in synaptic connectivity between GPe neurons but also in differences in intrinsic neuronal properties. It is found that the activation kinetics of different GPe neurons did not differ [59]. The physiological heterogeneity of GPe neurons can be replicated in simulations with different conductance densities in models, and the conductance of Na and K ions strongly influence the action potential (AP) widths [60]. A new fit for both types of GPe models is performed in a recent model by varying the conductance density parameters of Na and K in the GPe model to match the statistical properties of the firing spikes [61]. Its spontaneous firing rates are consistent with experiments, but the parameters which are related to exaggerated AP widths are removed. In contrast, the GPe electrophysiological data finds that Lhx6 GPe has a wider spike width than PV GPe [22,26]. Here, we adjust the conductance density parameters of Na and K (corresponding to g Na and g K in Eq. 2) of the Lhx6 GPe neuron model to make the simulation more compatible with physiology. The spontaneous firing is more consistent with physiological experiments ( Fig. 2A (a1), B) than the classic GPe model [48] ( Fig. 2A (a3), B). And the AP width of Lhx6 GPe is larger than that of PV GPe (Fig. 2C (c1,  c4)), while the classic GPe model (Fig. 2C (c3, c6)) and previous STN-GPe model ( Fig. 2C (c2, c5)) cannot find significant differences in AP widths.

Local circuitry in GPe network
Almost all GPe neurons have local axonal collateral branches [16], which form a complex network with local GABAergic synaptic connections. PV GPe neurons have longer local lateral branches and a greater number of axons than other GPe neurons [16], and PV GPe acts as a switch to control the neuronal activity of Lhx6 GPe [31,62]. However, the inhibition synapses between heterogeneous GPe have been modeled as simplified structures with the same number of connections, without considering the different numbers of synapses in physiology. Experimental data suggests that the average number ratio of pallidofugal axon buttons for each prototypic neuron to other GPe neurons is close to 5:3 [63]. To explore the effect of heterogeneous GPe connection structures on mutual inhibition, only the firing activities of PV GPe and Lhx6 GPe are considered. The excitatory input of STN is simplified to an excitatory white noise that brings the firing behavior of GPe closer to physiological states Here, two connection patterns between heterogeneous GPe, an approximate physiological synaptic scaling (Fig. 3B), and a simplified connectivity pattern in the model are simulated (Fig. 3A). The difference between the two connectivity modalities is that the number of physiological PV GPe-Lhx6 GPe channels (PV-Lh) is larger than that of Lhx6 GPe-PV GPe channels (Lh-PV) (Fig. 3B). The change in connection pattern is essentially consistent with the model in which PV-Lh has higher synaptic conductance than Lh-PV. In the biparametric diagram of PV-Lh synaptic conductance and applied current for Lhx6 Table 4 Parameters used for synaptic conductance (g X in mS Ácm À2 ), synaptic reversal potentials (E X in mV) for STN, PV GPe (PV) and Lhx6 GPe (Lh) neuron models

Synapses
Parameters Synapses Parameters GPe, the increase of connections for the same applied current leads to smaller firing rates of Lhx6 GPe neurons (Fig. 3D) than that of the model (Fig. 3C).
When the applied current of Lhx6 GPe is 20 lA.cm -2 , the discharge mode of PV GPe gradually enters betaband oscillation and then tends to lower frequency oscillation, and finally the discharging activity of Lhx6 GPe is completely inhibited by PV GPe (Fig. 3E, F). The change in connection mode does not affect the global trend of the firing pattern of PV GPe (Fig. 3F). Due to the enhanced inhibition of Lhx6 GPe, the parameter range of PV GPe appears to b oscillation expands, which is reflected in the region of higher oscillation index in the biparametric figure becomes larger (Fig. 3F). Hence, simplifying the synapse number in the model does not change the qualitative analysis and validates that the new parameters of Lhx6 GPe can reproduce the modulation of pathological beta oscillations by changes of mutual inhibition in previous studies [46]. Simulations are implemented in Matlab R2019a with equations solved using the forward Euler method with a time step of 0.01 ms. The initial value of the membrane potential is chosen to be close to the resting potential and contains random values, allowing the neurons within the nucleus to maintain an asynchronous firing state initially. Outcome measures are computed discarding the first 300 ms of the simulations to ensure that the activity stabilized from the initial conditions.
The control state represents the physiological state with the same firing rate and discharge pattern as in the physiological experiments [22,23]. This represents a healthy state compared to the Parkinsonian state.

Results
It is generally thought that STN-GPe circuit is a pacemaker in BG since the excitatory STN and inhibitory GPe are found to form a feedback system to engage synchronized bursting in vitro model [13]. The heterogeneous GPe circuit also generates abnormal oscillation in previous studies [46,47], which  [61], g Na ¼ 120; g K ¼ 30 (a3) [48]. C Effects of different parameters on the AP widths. Spontaneous firing activity (c1-c3) and the firing activity in GPe network (c4-c6) corresponding to different parameters (a1-a3). The black line represents the baseline AP width. The dark blue line and orange dashed line respectively represent the half-width of PV GPe and Lhx6 GPe, and the red Ã indicates a significant difference. (Color figure online) challenges popular belief. There is persistent GABAergic tension in the GPe following dopamine depletion [64], and aberrant b synchronization within BG has recently been found to be enhanced by afferent inhibition of PV GPe [43]. Both suggest that PV GPe may have an indispensable role in generating b oscillations. Recent experiments have found that inhibition of STN has fewer effects on cortical b oscillations [43], and STN has a stabilizing effect on the firing spike sequence of GPe [45]. It indicates that STN exhibits a different function contrary to classic STN-GPe model. In this section, we explore the roles of heterogeneous GPe neurons in producing beta-band oscillatory bursts, and the distinct effect of STN.

Mutual inhibitions between GPe are sufficient to produce b oscillation
There are significant differences in mutual inhibition between heterogeneous GPe [31,35]. Previous studies have found that altered mutual inhibition within GPe generates and regulates its b oscillations [46,47] (Fig. 3E). However, the association between afferent inhibition of PV GPe and b oscillations has not been explored in detail in the current study. As mentioned in Sect. 2.5, only the firing activities of PV GPe and Lhx6 GPe are considered and STN is simplified to an excitatory white noise. First, we explore the effect of mutual inhibition g pvlh and g lhpv changes within heterogeneous GPe networks. Figure 4A(a1) shows the firing properties of the PV GPe and Lhx6 GPe neurons under the control state, and the corresponding values of g pvlh = 1.2 and g lhpv ¼ 0: 2 are given in Table 4. Minimizing the PV-Lh synaptic strength individually (g pvlh ¼ 0:5), the Lhx6 GPe is less inhibited with its firing rate increasing (Fig. 4B(a2)), and the inhibitory projection to the PV GPe is enhanced. Bursting occurs in both PV GPe and Lhx6 GPe and spreads through Lh-PV and PV-Lh channels (Fig. 4A(a2)), thus extensive beta oscillations appear in heterogeneous GPe (Fig. 4C(a2)). Enhancing Lh-PV synaptic strength (g lhpv ¼ 0:8) alone also increases inhibitory projections to the PV GPe, decreasing the firing rate of the PV GPe, which in turn promotes a higher firing rate of the Lhx6 GPe (Fig. 4B(a3)). Only significant b oscillations arise within PV GPe (Fig. 4C(a3)), and Lhx6 GPe is subjected to strong inhibitory projections and b oscillations from the PV GPe, but its firing rate rises minimally (Fig. 4B(a3)) and there is no beta oscillation within Lhx6 GPe (Fig. 4A(a3)). Combining the above two patterns simultaneously increasing Lh-PV and decreasing PV-Lh synaptic strength (g pvlh ¼ 0:5, g lhpv ¼ 0:8), PV GPe is greatly inhibited, Lhx6 GPe is less inhibited in Fig. 4A(a4). Lhx6 GPe exhibited a higher firing rate than PV GPe ( Fig. 4B(a4)), generating strong b oscillations throughout the heterogeneous GPe loop (Fig. 4 A(a4), C(a4)). According to Fig. 4C, the dominant frequency of b oscillations shift toward higher frequencies after the increase in Lh-PV synapses. As a result, the elevated Lhx6 GPe firing activity and the decreased PV GPe firing activity suggest potential abnormalities in the dynamic behavior within GPe. It suggests that the increase of afferent inhibition to PV GPe leads to abnormal beta oscillations, is consistent with an experiment in which increased inhibitory input to PV GPe is a necessary condition for the emergence of anomalous b synchronization in BG [43]. In previously established models, the synaptic conductance of GPe-GPe may be elevated at lower dopamine levels [65], which generates pathological synchronous oscillations in STN-GPe circuit [66]. Previous studies of heterogeneous GPe have found that increased self-inhibitory synaptic conductance of GPe generates beta oscillations in the GPe microcircuit [46], but the role of increased self-inhibition of GPe has not been determined separately. Similar to Sect. 3.1, only PV GPe and Lhx6 GPe are considered in the loop, and the STN is modeled as white noise. Self-inhibition is an important component of inhibitory input of PV GPe, we determine the relation between inhibition on PV GPe and beta oscillation by separately elevating self-inhibitory synaptic conductance in PV GPe and Lhx6 GPe. In Fig. 5A and B, orange lines represent spike raster and LFP of Lhx6 GPe neurons while dark blue lines represent that of PV GPe. Figure 5A(a1) shows the firing activities of all PV GPe and Lhx6 GPe neurons with 0.25 s simulation under the control state, and the corresponding values of g pvpv = 1 and g lhlh ¼ 1 are given in Table 4. Comparing Fig. 5 A(a1) and (a2), the enhanced Lhx6 GPe-Lhx6 GPe (Lh-Lh) synaptic conductance (g lhlh ¼ 2) causes a decrease in the firing rate of Lhx6 GPe. And PV GPe still maintains high-frequency discharging with more regular firing activity ( Fig. 5 A(a2)). Due to Lhx6 GPe being strongly inhibited by PV GPe, changing the self-inhibitory synapse of Lhx6 GPe alone cannot disrupt the firing pattern of PV GPe (Fig. 5 A(a2), C(a2)). Lhx6 GPe neurons inhibit the synchronous tendency of GPe itself and promote irregular activities in PV GPe neurons, in agreement with previous studies [58] that neuronal heterogeneity can reduce the synchronization of the nucleus itself.
After blocking the PV-Lh channel (g pvlh ¼ 0) in mutual inhibition to remove the effect of PV GPe on Lhx6 GPe and enhancing the Lh-Lh synaptic conductance (g lhlh ¼ 2), the firing rate of Lhx6 GPe increases dramatically as well as significant bursting behavior in both Lhx6 GPe and PV GPe (Fig. 5B(b4), C(b4)). Enhanced inhibition of PV GPe by Lhx6 GPe, and propagation of bursts of activity at frequencies in the beta band into the PV GPe (Fig. 5D(b4)). Therefore, joint self-inhibition of Lhx6 GPe and mutual inhibition may greatly affect the dynamics of GPe.
In the biparametric plots Fig. 5E, the increase of g pvlh and g lhlh leads to a decrease in the firing rate of Lhx6 GPe. Then, PV GPe is subject to a weaker inhibitory projection from Lhx6 GPe, causing a decrease in the number of spikes within a burst of PV GPe and an increase in spike train variability. When both PV-Lh and Lh-Lh synaptic conductance are at the small value, Lhx6 GPe is less inhibited and exhibits a high firing rate, providing a strong inhibitory projection on PV GPe, resulting in bursting activity in PV GPe with the main frequency located in the b band ( Fig. 5E(II)). When both PV-Lh and Lh-Lh synaptic conductance are at large values, the discharging behavior of Lhx6 GPe is completely inhibited by PV GPe, and PV GPe maintains c band firing (Fig. 5E(I)). When the combination of g pvlh and g lhlh is located in the white region in Fig. 5E, neurons in PV GPe maintain synchronous bursting. Along the white region in Fig. 5E from left to right, g lhlh decreases and g pvlh increases, and the firing rate of Lhx6 GPe decreases due to receive enhanced inhibition from PV GPe. The firing pattern of Lhx6 GPe breaks away from bursting activity (II) to single spike firing (I), and then the synchronous bursting behavior transits to irregular discharging (Fig. 5E). When Lhx6 GPe transitions to single spike firing in the white region ( Fig. 5E (I)), PV GPe also behaves as Lhx6 GPe and pallidal network is in a strongly synchronized state, the main frequency moves from b (II) to c band (I) (Fig. 5E).
According to the biparametric plots Fig. 5F, increasing Lh-Lh synaptic conductance (g lhlh ) cannot disrupt the high-frequency discharge of the PV GPe, and the maintenance of a lower discharge frequency of Lhx6 GPe when g lhpv is small (I) (Fig. 5A(a2)). As the Lh-PV synaptic conductance (g lhpv ) increases, Lhx6 GPe strongly inhibits the PV GPe and significant bursting activities occur in PV GPe. When g lhlh [ 1, the firing mode of PV GPe breaks away from c band firing to b band bursting (II) and then moves to below b band bursting (III) (Fig. 5F). When g lhlh is at a lower level, PV GPe changes from high-frequency firing into a bursting activity at the frequency of b band(II) and then changes to single spike firing with the primary frequency of c band (I). The presence of single spike firing in Lhx6 GPe when g lhlh and g lhpv are approximated, meanwhile bursting activity with two spikes occurs in PV GPe (II) (Fig. 5F). It is confirmed that the firing rate of Lhx6 GPe increases, and enhances afferent inhibition on PV GPe leading to abnormal b bursting.
Since the altered PV-Lh channel can generate extensive b oscillations among GPe, it is possible to transfer the b oscillations to the pallidal network. Considering g pvpv ¼ 1:5 and g pvlh ¼ 0:5, it appears that the firing rate of Lhx6 GPe is greater than that of PV GPe (Fig. 5B(b3)). PV GPe is subjected to stronger inhibitory projections from Lhx6 GPe, causing the dominant frequency of bursting behavior to shift toward lower frequencies (Fig. 5D(b1), (b3)), but the oscillation index decreases compared to increasing PV-PV synapses (g pvpv ¼ 1:5) alone (Fig. 5C(b1),  (b3)). The effects of joint self-inhibition and mutual inhibition differ between PV GPe and Lhx6 GPe because the afferent inhibition of PV GPe changes differently. Indeed, heterogeneous GPe neurons play different roles in aberrant beta-band activity, only the enhancement of PV-PV synaptic promotes abnormal beta behavior. While the enhancement of Lh-Lh essentially reduces inhibitory input to PV GPe, reflecting the tendency of neuronal heterogeneity to reduce synchronization of the nuclei itself [58] STN) and Lh-PV channels (g lhsn ¼ g lhpv ¼ 0) in STN-GPe to exclude the effect of neuronal heterogeneity. From the STN-PV bidirectional synaptic biparametric diagram (Fig. 6A), we find that enhancing the PV GPe-STN (PV-STN) (g pvsn ) and attenuating the STN-PV channel (g snpv ) could promote b oscillations within GPe.
Enhancing the STN-PV (g snpv ¼ 1) synapse promotes the firing of PV GPe neurons (Fig. 6C(c1)), and indirectly increases the self-inhibitory input of it to produce abnormal b oscillations (Fig. 6B, D(c1)). However, the STN is still strongly inhibited by the PV GPe ( Fig. 6C(c1)), so this faciliatory effect is weaker. On top of this, weakening the PV-STN synapse (g pvsn ¼ 0:5) promotes higher firing activities of STN (Fig. 6C(c3)) to further enhance the excitatory input. This indirect facilitation leads to higher firing rates of PV GPe (Fig. 6C(c3)) and generates b oscillations (Fig. 6B, D(c3)). The addition of Lhx6 GPe neurons suppress the abnormal oscillations in PV GPe (Fig. 6B, D(c2, c4)). However, it only has a slight effect on the firing rate of STN (Fig. 6C(c1, c3)), indicating that Lhx6 GPe may not affect b oscillations in PV GPe by directly changing the firing activity of STN.
To explore how Lhx6 GPe affects abnormal rhythms inside PV GPe, Lh-STN and Lh-PV synapses are blocked separately. The Lh-STN channel directly b Fig. 5 A The spike raster of PV GPe (dark blue) and Lhx6 GPe (orange) neurons with 0.25 s simulation under control state of pallidal network (a1), g lhlh ¼ 2 and g pvlh ¼ 1:2 (a2). B The LFP of PV GPe (dark blue) and Lhx6 GPe (orange) neurons with 1 s simulation under g pvpv ¼ 1:5 and g lhpv ¼ 0:2 (b1), g pvpv ¼ 1:5 and g lhpv ¼ 0 (b2), g pvpv ¼ 1:5 and g pvlh ¼ 0:2 (b3) and g lhlh ¼ 2 and g pvlh ¼ 0 (b4). The oscillation index (C) and spectral power analysis (D) of PV GPe correspond to states (a1-a2, b1-b4). The area between the dotted lines represents affects the firing activity of STN, while the Lh-PV channel increases the afferent inhibition of PV GPe. When g snpv ¼ 1 (Fig. 6C(c1)), the beta oscillations in PV GPe and firing activity of STN are slightly attenuated by further inhibition of Lh-STN synapse (g lhsn ¼ 0:2), while the addition of the Lh-PV synapse (g lhpv ¼ 0:2) directly increases the inhibitory input of PV GPe to promote b activities (Fig. 6B). For the simulated state in Fig. 6C(c3) (g pvsn ¼ 0:5; g snpv ¼ 1), the STN-Lh synapse cannot significantly alter the discharge activity and affect b oscillations in STN-GPe network. The addition of Lh-PV inhibits PV GPe activity and thus affects b oscillations (Fig. 6B). Lhx6 GPe suppresses abnormal oscillations generated by STN and PV GPe, and thus neuronal heterogeneity has the tendency to inhibit correlations outside the nucleus.

STN suppresses the anomalous oscillatory behavior in GPe
Although STN-GPe has been generally assumed as a generator of b oscillations, optogenetic suppression of STN has been found to rarely affect broad b oscillations and exacerbate abnormal activity in GPe recently [43,44]. STN could stabilize the pallidal spike sequences and prolong movement time [45]. This indicates that STN has a different role in anomalous beta oscillations than that in previous studies. Based on the previous simulations, we choose several representative states in self and mutual inhibition variations to explore the effect of adding STN on the anomalous oscillatory dynamics in GPe. In the control state of the pallidal network, the spectrogram of the PV GPe has no significant power in the b band and the dominant frequency is in the c band ( Fig. 7A(a1)).
When g pvpv ¼ 1:5, g pvsn ¼ 1, g lhsn ¼ 0:2, the increase of PV-PV leads to a significant bursting behavior in PV GPe (Fig. 7B(b1)). Meanwhile, STN neuron exhibits regular spiking and Lhx6 GPe is strongly inhibited by PV GPe showing sparse firing activity as mentioned before (Fig. 5B(b1)). After adding STN-PV to the heteronomous GPe network (g snpv ¼ 0:3), the firing activity of STN and Lhx6 GPe become irregular, the number of spikes in the burst decreases, and the spike sequences become more variable in PV GPe (Fig. 7B(b3). Introducing STN under pathological conditions (g snpv ¼ 1, g snlh ¼ 0:5, g pvsn ¼ 0:5, g lhsn ¼ 0:1) to the pallidal network simulated in Fig. 7B(b1), the occurrence of burst in PV GPe significantly weakens (Fig. 7B(b2). And efferent inhibitions of PV GPe decrease, leading to higher firing rates in STN and Lhx6 GPe. To combine these two conditions, STN-PV and PV-STN under pathological conditions (g snpv ¼ 1, g pvsn ¼ 0:5) are considered in Fig. 7B(b4). The beta band power of PV GPe and the firing rate of Lhx6 GPe decrease through blocking STN-Lh and Lh-STN ( Fig. 7B(b2, b4)). The pathological connections between STN and GPe relieve STN from strong inhibition from PV GPe, further inhibiting b bursting behavior in PV GPe. Hence, STN strongly affects the firing activities of PV GPe with enhanced selfinhibition, especially through the STN-PV channel.
Mutual inhibition of PV-Lh attenuates (g pvlh ¼ 0:5) leading to a significantly higher power in the low b band (Fig. 7A(a2)). While adding the connection of STN to GPe, the b band firing activity is significantly weaker (Fig. 7A(a3)). In Fig. 7A(a4), PV GPe maintains b band behavior and the dominant frequency shifts towards high beta bands with blocking STN-PV synapse, suggesting STN-PV is a main pathway for STN suppressing the b oscillatory activity in GPe. When these two changes in inhibition are considered jointly (g pvpv ¼ 1:5, g pvlh ¼ 0:2), PV GPe has higher energy in the low b band (Fig. 7A(a5)). Introducing excitatory projections of the STN to heterogeneous GPe network, the b band power of the PV GPe and strong synchronous discharge activity decrease (Fig. 7A(a6)). Compared to Fig. 7A(a6), blocking the STN-Lh channel leads to a more regular bursting behavior and significant b band power in Fig. 7A(a7). However, the connection of STN to GPe is introduced under pathological conditions (g snpv ¼ 1, g snlh ¼ 0:5, g pvsn ¼ 0:5, g lhsn ¼ 0:1), and STN strongly affected the firing activity of PV GPe, changing from bursting in the low b band to irregular bursting in the high b band (Fig. 7A(a8)).
When the b frequency excitatory and inhibitory oscillation are delivered to GPe (Fig. 8A(a2, a3)), b propagates throughout the loop in both control and pathological conditions, while the oscillations are highly amplified in pathological networks ( Fig. 8B(a2, a3)). The oscillatory inputs applied to GPe change its firing activity little and the spectrograms indicate a slight elevation of b band power in control networks ( Fig. 8 B(a2, a3)). However, these inputs modulate the b band firing activity of the STN, and these activities may be entrained by the stimulus frequency. After applying excitatory and inhibitory oscillatory input to GPe, STN still maintains a high level of power near the external oscillation input frequency. In pathological states, changes in mutual inhibition and self-inhibitory synaptic connections (g pvpv ¼ 1:5., g pvlh ¼ 0:2) within the GPe promote b. rhythm activity and propagate within the STN-GPe loop, which in turn affects BG dynamics. The generated b oscillation modulates the external b inputs in the pallidal network and the firing activity of GPe neurons changing from weak firing at about 20 Hz (frequency of external input) to strong activities in the broad beta band. Meanwhile, the synaptic coupling of the GPe to the STN decreased (g pvsn ¼ 0:5) directly affects the firing activity of the STN, further promoting the propagation of b. rhythm activity.
In summary, PV-PV is responsible for the enhancement of beta rhythms in the loop. PV-Lh promotes b oscillations in the pallidal network, especially in Lhx6 GPe neurons, and this propagation feeds back to regulate oscillations within PV GPe. Altered inhibitory PV-STN increases the firing rate of STN and promotes the propagation of abnormal oscillations in the STN-GPe loop.
All major results are shown with different parameters in Table 5, and all parameters changes in results refer to Eq. 11. Table 5 All parameters change in Results and the corresponding descriptions or major results

Discussion
This paper proposes an STN-GPe model incorporating Lhx6 GPe units to reveal a novel mechanism through which the heterogeneous neuron influences the functional coupling in the subthalamopallidal network (Fig. 1). Results reveal that enhanced afferent inhibition of the PV GPe generates b oscillations (Fig. 4), which are necessary for the emergence of abnormal oscillations in GPe. Neuronal heterogeneity promotes irregular firing behavior of GPe to inhibit synchronous activity ( Fig. 5A(a2)) and reduces the correlation between STN and GPe (Fig. 6B, C). Besides, STN neurons in control stabilize the firing sequence of GPe (Fig. 7B), and it suppresses abnormal synchronous activity in GPe (Fig. 7A, C). STN inhibits rhythmic oscillations in a heterogenous pallidal network with STN-PV dominating the effect (Fig. 7C). External oscillations from the STN can be transmitted to GPe, but the effect is weak and the appearance of b oscillations is almost absent in GPe (Fig. 8C). In contrast, oscillations transmitted from the GPe propagate in subthalamopallidal network and are significantly enhanced as well as extensively spread in the pathological state (Fig. 8D). The self-inhibitory enhancement of PV GPe significantly amplifies the incoming b oscillations, while synaptic connections PV-Lhx6 and PV-STN mainly influence the extent of b oscillation propagation. As a whole, the simulation results show that PV GPe promotes pathological synchronous oscillations in the PD state, and STN operates as a nucleus capable of preventing abnormal rhythmic. Several computational models have been used to study infectious diseases, which build differential and difference equations based on susceptible, infected, recovered and cross-immune individuals to explore the dynamics of different populations [71][72][73]. In contrast, computational model of PD, which are based on different study subjects, address neuronal modeling of deep brain nuclei within individuals to study the dynamic changes of neurons and nuclei [14,15,42]. Our study probes anew the role of the clascal STN-GPe loop in b oscillations, similar to previous results as generators and amplifiers, but the roles of the STN and GPe are different from previous perceptions [14,48]. In our study, we find that the role of the STN is weaken the anomalous synchronization activity within GPe (Fig. 7). Consistent with the results in the experiment [45], STN has a stabilizing effect on pallidal spike sequences, providing a new explanation for the selection of STN as a DBS target. And the propagation effect of STN is blocked after DBS and strongly dampens the anomalous oscillations in GPe.
The firing rates of heterogeneous GPe neurons differ significantly in both physiological and pathological states [21], suggesting different roles of heterogeneous GPe in aberrant beta-band activity. The change in the firing rate of heterogeneous GPe neurons could reflect synaptic changes, and the greater firing rate of Lhx6 GPe than PV GPe [37] is a case that could indicate the generation of b oscillations in GPe (Fig. 4A(a4)), but not a necessary condition. Abnormal alterations in synaptic coupling within the heterogeneous pallidal network can affect beta oscillations in PD (Figs. 4,5). For self-inhibition in homogeneous GPe, only the enhancement of PV-PV synaptic promotes abnormal beta behavior ( Fig. 5B(b1)), while the enhancement of Lh-Lh facilitates irregular firing of PV GPe (Fig. 5A(a2)), reflecting the tendency of neuronal heterogeneity to reduce synchronization of the nuclei. For mutual inhibition between heterogeneous GPe, the increase of afferent inhibition to PV GPe leads to abnormal beta oscillations (Fig. 4). Indeed, the origin of b oscillations is probably the inhibitory BG central GPe. The PV GPe exhibits high-frequency firing and is subject to enhanced afferent inhibition leading to its persistent b oscillations, and attenuating inhibitory inputs can control b oscillations well.
However, some limitations of the model remain. First, the pathological oscillations considered here are established by changing the weight of synaptic coupling within the GPe network or between STN-GPe. And the STN in the model only shows sparse firing with no bursting behavior due to the enhanced inhibition of GPe, which may limit the synchronous oscillation phenomenon of the STN-GPe. In further studies, the model could be extended. We only consider the local STN-GPe loop which incorporates heterogeneous GPe, while STN-GPe is influenced by more inputs in the basal ganglia-thalamus-cortical (BGTC) circuit. Heterogeneous STN-GPe loops could be integrated into the BGTC model, and the dynamical behavior may be more abundant tending to the real situation. The simplified modeling of anomalous oscillatory inputs ignores the interaction of external nuclei with STN-GPe and the facilitating and regulating effects on oscillatory inputs. Actually, expiments [8] and computational studies [42,70] have revealed that hyperdirect pathway has an important role in propagating and promoting abnormal oscillations. Through reproducing the physiological phenomenon of exaggerated phase-amplitude coupling (PAC), it is found that the hyperdirect path from cortex to STN drives the beta oscillations of STN [42]. The multiple effects of cortical inputs to the STN on synchronized beta-band oscillations were explored through resonance between external inputs and internally generated beta oscillations accordingly [74]. However, GPe also receives direct excitatory cortical input, and sends inhibitory projection to cortex [29,33], suggesting that the cortical-GPe pathway may also be an hyperdirect channel. In the future, the hyperdirect pathway may extend to encompass the direct connection between the cortex and GPe leading to a complex coupling effect in the cortex-STN-GPe circuit. The individual anisotropy of Parkinson's patients has a strong influence on the origin of beta oscillations, making it difficult to test the hypothesis of the origin. This paper may be important for future hypotheses on the pathogenesis of PD based on neuronal heterogeneity and provide a new perspective on the precise mechanism by which STN affects pallidal activity.