Suppression of beta oscillations by delayed feedback in a cortex-basal ganglia-thalamus-pedunculopontine nucleus neural loop model

Excessive neural synchronization of neural populations in the beta (β) frequency range (12–35 Hz) is intimately related to the symptoms of hypokinesia in Parkinson’s disease (PD). Studies have shown that delayed feedback stimulation strategies can interrupt excessive neural synchronization and effectively alleviate symptoms associated with PD dyskinesia. Work on optimizing delayed feedback algorithms continues to progress, yet it remains challenging to further improve the inhibitory effect with reduced energy expenditure. Therefore, we first established a neural mass model of the cortex-basal ganglia-thalamus-pedunculopontine nucleus (CBGTh-PPN) closed-loop system, which can reflect the internal properties of cortical and basal ganglia neurons and their intrinsic connections with thalamic and pedunculopontine nucleus neurons. Second, the inhibitory effects of three delayed feedback schemes based on the external globus pallidum (GPe) on β oscillations were investigated separately and compared with those based on the subthalamic nucleus (STN) only. Our results show that all four delayed feedback schemes achieve effective suppression of pathological β oscillations when using the linear delayed feedback algorithm. The comparison revealed that the three GPe-based delayed feedback stimulation strategies were able to have a greater range of oscillation suppression with reduced energy consumption, thus improving control performance effectively, suggesting that they may be more effective for the relief of Parkinson’s motor symptoms in practical applications.


Introduction
There are several hypotheses for the development of Parkinson's disease (PD) symptoms, and it is widely believed that the development of the symptoms is closely related to abnormalities in dopamine neurons.For example, the neuronal death hypothesis, the oxidative stress hypothesis, and the inflammation hypothesis [1] suggest that oxidative stress or inflammation causes neuronal damage or death, respectively, resulting in decreased dopaminergic levels and PD symptoms.The neuronal synaptic dysregulation hypothesis suggests that neuronal synaptic dysregulation leads to decreased dopaminergic levels and the development of PD symptoms.Excessive rhythms in the pathological beta (β) band  in the brain disrupt the transmission function of thalamic (Th) neurons, leading to PD symptoms such as body tremor and dyskinesia [2][3][4][5][6], thus suppressing β oscillations that play a crucial role in alleviating the local firing state of neurons and thus effectively treating PD.Deep brain stimulation (DBS), a neuromodulation technique that regulates neuronal activity by applying electrical pulses to neuronal populations in the basal ganglia (BG) [2,[7][8][9][10], is a standard method for treating PD.However, the fact that permanent strong high-frequency stimulation is required for the implementation of DBS leads to its poor performance in terms of efficient use of energy [11][12][13][14], which may be risky in clinical applications.Therefore, the further development of DBS is constrained.
Delayed feedback control strategy is a closed-loop control method to achieve suppression of β oscillations through a desynchronization mechanism, which can effectively alleviate PD symptoms by modulating the stimulation signal according to neuronal activity [15].Since Pyragas [16] proposed that it has a desynchronization mechanism [17][18][19], the delayed feedback strategy has received wide attention.Rosenblum and Pikovsky [20] developed the view of Pyragas by proposing a delayed feedback control method that removes collective synchronization, and they found that time-delayed feedback can be used to suppress pathological brain rhythms.Hövel and Schöll [21] explored the effect of timedelayed feedback control on steady-state stability and demonstrated that time-delayed feedback control is equivalent to a tool used to stabilize unstable states.Popovych et al. [17] proposed a nonlinear delayed feedback stimulation scheme, resulting in better adaptation to parameter changes for desynchronization.Later, Popovych et al. [22,23] applied the delayed feedback control scheme in a computational model for PD including the subthalamic nucleus (STN)-external globus pallidum (GPe) circuit and have shown that different types of delayed feedback algorithms are effective in eliminating synchronization of neuronal populations.Furthermore, they found that nonlinear differential delayed feedback has a larger parameter space, while linear differential delayed feedback strategies are more effective in suppressing more intense synchronization.Daneshzand et al. [24] found solutions that could improve desynchronization efficiency while reducing energy consumption by adjusting the frequency of the stimulus signal.Popovych and Tass [25] in their recent work compared the desynchronization effects of continuous and adaptive high-frequency DBS in a STN-GPe neuronal network model, and they found that adaptive DBS was more effective in suppressing abnormal synchronization.
Much progress has been made in the work on optimizing the delayed feedback strategy; however, the previous work has mostly focused on algorithm improvements where the stimulus signal is only obtained from the local field potential (LFP) of the neurons in the STN [22], with little involvement of the GPe.Considering that it is still challenging to further improve the inhibition effect with lower energy expenditure, it is hoped that more effective parameter regimes can be found at reduced energy consumption to address the difficulties in parameter setting.Therefore, GPe was also considered as the stimulus target and detection target, combined with STN to form four types of target combinations [26], and the advantages of GPe-based delayed feedback schemes were demonstrated by comparison.
Most of the previous studies were based on STN-GPe networks [22,23,25].In order to simulate a more realistic brain and to take into account the important roles of the cerebral cortex, BG, and PPN, we developed a cortex-basal ganglia-thalamus-pedunculopontine nucleus (CBGTh-PPN) neuronal network model.We understand that biophysical models [27] can represent neuronal features in detail, including the electrical and chemical properties of neurons, and in addition, can represent neuronal interactions well.However, using biophysical models to simulate large-scale neural networks is computationally intensive and difficult to explain the complex dynamics of individual neurons.Population-level models [28] differ in many ways from biophysical models in that instead of modeling individual neurons, they represent the population of neurons as a statistical distribution.The advantage is that it is simpler, but has the potential to oversimplify complex neural interactions.The modeling philosophy of the neural mass model [29,30] differs from both biophysical and population-level models.It is a representation of a large population of neurons as a single, simplified entity.Although the details of individual neurons may be individual, the model is simpler, computationally efficient, and commonly used to simulate large-scale brain networks.To better simulate neuronal LFPs, we chose to use the neural mass model.
In this study, we first developed CBGTh-PPN neuronal network model and described four control strategies.The model developed can reflect the interaction between some neuronal populations in the brain and can acquire LFP signals from STN and GPe.Then, we explored the effects of each of the four schemes on alleviating the pathological activity.In addition, we analyzed the effects of feedback gain and time delay in the suppression of pathological oscillations by the delayed feedback scheme.

Structure of the cortex-basal ganglia-thalamus-pedunculopontine nucleus closed-loop system
The CBGTh-PPN circuit is shown in Fig. 1a and contains three parts: cortex, basal gangliathalamus (BGTh), and PPN, which are constitutively connected through excitatory and inhibitory projections between them.The cortex contains three populations of neurons, namely excitatory interneurons (EI), pyramidal neurons (PY), and inhibitory interneurons (II).Yan et al. [31] investigated the function of EI by bifurcation analysis and found that EI induces transition behaviors under both Hopf type and fold limit cycle bifurcations, suggesting that EI has an impact on neurological diseases such as PD and epilepsy.Therefore, EI was added to the cortex.The BGTh contains seven populations of neurons, namely the fast-spiking interneuron (FSI), medium spiny-projection neuron (MSN) (includes D1 and D2 types), STN, GPe, internal globus pallidum (GPi), and Th.In addition, it was clinically shown that despite stimulation of the STN/GPi, early-onset gait disturbances and treatmentresistant gait freezing may benefit from additional low-frequency stimulation of the PPN [32][33][34][35].This suggests that PPN has the potential to improve gait freezing and reduce falls in some PD patients.Given the complexity of the control of PD by PPN, it is considered to be a separate component.Furthermore, the inclusion of PPN and EI in the model is also expected to simulate a comprehensive model of PD.
In the cortex, the PY nucleus forms an excitatory-inhibitory network with II and an excitatory-excitatory network with EI [36].The PY nucleus is also the output nucleus of the cortex, which regulates the direct and indirect pathways in the BG by projecting synaptic connections to the BGTh and forms the hyper-direct pathway with the STN [6,36].In BGTh, the structure of the striatum model proposed by Moran et al. [37] is referenced; it contains FSI and MSN [38][39][40].MSNs express two types of dopamine receptors (D1 and In the striatum, all three neuronal populations have selfinhibitory connections and the MSNs receive inhibitory projections from the FSI.In the direct pathway, MSN D1 has a direct effect on GPi.In the indirect pathway, MSN D2 regulates the output of BG through the STN-GPe network.According to biophysical properties, GPe has inhibitory projections to FSI [38], forming local circuits between GPe, FSI, and MSN D2 .In addition to the populations of neurons in the striatum that form a local circuit with GPe, STN forms an essential local neuronal circuit with GPe, which are the two nuclei used as stimulation targets in the delayed feedback strategy.Among them, the STN has excitatory projections to the GPe, while the GPe has inhibitory projections to itself and the STN [41].GPi is the output nucleus of BG, which receives excitatory projections from the STN and inhibitory projections from GPe and MSN D1 [36,38].Finally, the Th receives inhibitory projections from the GPi and sends excitatory projections to the EI nucleus in the cortex, thus forming a cortex-basal ganglia-thalamus (CBGTh) neural circuit.
There is also a connection between the PPN nucleus and the BGTh, which forms an excitatory-excitatory neural circuit with the STN and an excitatory-inhibitory neural circuit with the GPi [34].At the same time, the PPN also sends excitatory projections to the Th [32,33,35], thus eventually forming a CBGTh-PPN neural circuit.Figure 1b shows a simplified depiction of this system.

Neural mass model of cortex-basal ganglia-thalamus-pedunculopontine nucleus neural circuit
A neural mass model is commonly used to describe the behavior of neuronal populations in the brain.It is based on the idea that individual neurons can be grouped into populations, and the activity of each population can be described by a set of differential equations [42,43].In the CBGTh-PPN system, each neuronal population is represented by a linear block and a nonlinear block, as shown in Fig. 2. They are represented by T i and S i functions [44,  45], respectively, where i ∈ {EI, PY, II, FSI, MSN D1 , MSN D2 , STN, GPe, GPi, Th, PPN}.T i converts the average pulse concentration x i ( i ∈ {EI, PY, II, FSI, MSN D1 , MSN D2 , STN, GPe, GPi, Th, PPN}) picked up from other nodes to the average membrane potential y i ( i ∈ {EI, PY, II, FSI, MSN D1 , MSN D2 , STN, GPe, GPi, Th, PPN}).S i converts the average membrane potential y i to the average pulse density x i and delivers it to other nodes via excitatory or inhibitory projections.
As shown in Figs.1a and 2, there are 31 connections between different populations of neurons in the CBGTh-PPN model.The parameter K i−j indicates the strength of coupling between presynaptic nucleus i to postsynaptic nucleus j .Because of the presence of many neuronal populations in the model, all connections are divided into internal and external connections for convenience, as shown in Fig. 1b.Specifically, internal connection K 1 contains the coupling strengths within the cortex and within BGTh, denoted by K Cor and K BGTh , respectively.The external connections contain the coupling strength between the cortex and BGTh ( K 2 ) and the coupling strength between BGTh and PPN ( K 3 ).Moreo- ver, K 2 contains K Cor-BGTh (i.e., the coupling strength from cortex to BGTh) and K BGTh-Cor (i.e., the coupling strength from BGTh to cortex), and K 3 contains K BGTh-PPN (i.e., the cou- pling strength from BGTh to PPN) and K PPN-BGTh (i.e., the coupling strength from PPN to BGTh).Inspired by Jansen and Rit [43], the 31 connections between different neuronal populations in our model were all determined by some ratio of these six basic parameters, as shown in Table 1.When the neural state changes, the six coupling parameters are able to reflect the trend of the coupling strengths within or between the systems.
To characterize the dynamics of individual nucleus in the CBGTh-PPN neural network, the linear and nonlinear elements of each nucleus were modeled using second-order transfer functions and sigmoid functions, respectively [36,46].In detail, T i takes the form of where s is the Laplace variable, X i (s) and Y i (s) are the Laplace trans- formations of x i and y i , respectively, and R i and r i are the parameters of the linear block.Based on control theory, the transfer function is a mathematical representation describing the relationship between the input and output of a linear element T and can therefore be trans- formed into a differential equation to gain the kinetics in the time domain.Thus, the transfer function can be rewritten as , using the differential operator d∕dt and replacing the Laplace variable s to obtain Then, ̇y is used in place of dy∕dt , and ÿ is used in place of d 2 y∕dt 2 , and we can obtain . Combining the expressions of x i and the coupling strength K i−j , the differential equation is used to establish the connection between y i and x i (as shown in Eq. ( 1)), and the dynamics of each nucleus can be obtained by numerical solution.The x i can be expressed in terms of the nonlinear operator S i , which is given by the sig- moid function where b i and c i are the parameters of S i .
For the PPN: where y i j in Eq. ( 1) represents the output from presynaptic nucleus i to postsynaptic nucleus j .= ±1 , = 1 represents the excitatory connection between i and j , = −1 rep- resents the inhibitory connection between the two, and N represents the total number of (1) Coupling strengths between BGTh and PPN presynaptic nuclei owned by the prominent postsynaptic nucleus j .It should be noted that m i in Eq. ( 2) represents the stimulus signal obtained from the delayed feedback control policy.Table 2 shows the parameters in the model, referring to the model parameters of Liu et al. [47,48].The equations were solved using the forward Euler method with a time step of dt = 0.1 ms and a duration of 5 s.The simulation experiments were performed in MATLAB 2018b.

Description of four delayed feedback stimulation schemes
In the study of Holt et al. [49,50], the feedback signal was only obtained from the LFP of the STN and then the stimulus was only delivered to the STN itself, without considering the GPe as the feedback signal or the stimulus target.Therefore, we refer to the strategy that does not involve GPe as the conventional strategy.Based on the established model, we will discuss four delayed feedback stimulation schemes, including the traditional STNbased scheme and three GPe-based schemes, as shown in Fig. 3.We would refer to panels a-d in Fig. 3 as schemes A-D in the rest of the section.The linear delayed feedback algorithm [22,51] is as follows: where y j ( j ∈ (STN, GPe) ) is the LFPs generated by nucleus j , m i ( i ∈ (STN, GPe) ) is the stimulus signal received by population i , G represents the feedback gain, and represents the time delay.The feedback signal is generated by the LFPs of the STN in scheme A and scheme B, and is generated by the LFPs of the GPe in scheme C and scheme D, while the feedback signal is received by the STN in scheme A and scheme D and by the GPe in scheme B and scheme C. The strategy was considered effective when the stimulation signal was able to suppress the pathological oscillations of PD.

Evaluation indicators
In order to verify the suppression effect of the four delayed feedback control methods, the following three metrics were calculated and compared.Power spectrum: The power spectrum can show the oscillation characteristics in the PD state.The Fast Fourier transform based on the LFP of each neuron population can generate the corresponding power spectrum, and then, the direct current component is removed by subtracting the average LFP [52].Here, a time period of 0.5 to 1.5 s is chosen.Oscillation index ( I j OS ): In order to more accurately detect the level of oscillation of the nucleus after receiving a stimulus, we define an index of oscillation.We will perform the calculation at the fourth to fifth second after the implementation of the scheme, which allows us to effectively exclude interference.The formula is as follows: Energy consumption index ( I EN ): Energy consumption is an index to evaluate the per- formance of the delayed feedback stimulation strategy, and it is calculated as follows:

Oscillations of the CBGTh-PPN system in normal and Parkinsonian states
In order to understand the control effect of linear delayed feedback, we expected the model to reproduce the oscillations in the pathological state of PD.Here, Since there are many neu- ronal populations in the model, only the results for PY, STN, and PPN are shown here.From the connections between nucleus, we can see that the coupling strengths between The feedback signal is generated from the LFPs of the STN and the simulation is delivered to the GPe.c Scheme C. Similar to scheme A, the delayed feedback signal is obtained from the LFPs of the GPe and received by the GPe.d Scheme D. Contrary to scheme B, the feedback signal is generated from the LFPs of the GPe and the simulation is delivered to the STN the cerebral cortex and BGTh are related to K 1 and K 2 , but not to K 3 .For the PPN, K 3 is connected to both K 1 and K 2 .So we first explored the relationship between K 1 and K 2 , and then explored the relationship between K 3 and K 1 and K 2 , respectively.As shown in Fig. 4, in the normal state, the values of K 1 corresponding to both PY and STN are small, and neither nucleus showed enhanced oscillations.Meanwhile, to better simulate the Parkinsonian oscillations in the cortex, we chose larger values of K 1 .Since K 2 has less influence on the state of the nucleus and considering the combined effect of K 2 on PY and PPN, the selected values of K 2 are larger in the normal state than in the pathological state.The smaller K 3 was selected to simulate the normal state according to the variation pattern of K 3 in biology.Therefore, the values corresponding to ( K 1 , K 2 , K 3 ) in the normal and pathological states are (3, 9.75, 3) and (60, 5, 10), respectively.In Fig. 4, the values of K 1 , K 2 , and K 3 are marked with different colors of A, B, and C for the normal and Parkin- sonian states, respectively.
By changing the values of the coupling strengths K 1 , K 2 , and K 3 , the normal and patho- logical states of each nucleus can be demonstrated, as shown in Fig. 5. From Fig. 5a, it can be seen that, unlike the normal state, all nucleus in the pathological state exhibit significantly enhanced oscillations.Spectral analysis shows that in the pathological state, as the spectral power of the neuronal population increases, more peaks appear in the band, as shown in Fig. 5b.This is consistent with the experimental results of Asadi et al. [53] and McCarthy et al. [54,55].Here, we only show the peak power and peak frequency of the neuronal populations in BGTh and PPN.Since the oscillations are closely related to the dynamic characteristics of the limit cycle, the original stable equilibrium point changes to a stable limit loop when the coupling strengths is enhanced.Figure 5c gives the trajectories in the phase diagram corresponding to the four neuronal populations of STN, Th, GPe, and PPN, from which the activities of the populations in normal and pathological states can be more clearly seen.It can be seen that the model we built can reproduce the neural characteristics of all nucleus in normal and Parkinsonian states.Therefore, we will further investigate the suppression of oscillations in Parkinsonian state by delayed feedback stimulation schemes based on STN and GPe nucleus.

Implementation of four delayed feedback stimulation schemes
We next investigated the effects of four different delayed feedback stimulation schemes of A, B, C, and D on the β oscillations in the CBGTh-PPN model.The stimulation signals of the four different schemes are shown in Fig. 6, corresponding to the four schemes in Fig. 3, respectively.As can be seen from Fig. 6, all four stimulus signals fluctuate above and below zero and show different waveforms.Considering the strong excitatory and inhibitory connections between neuronal populations, the feedback gain G corresponded to the four schemes were set to 0.4, -0.3, -0.8, and 0.4, respectively.From Fig. 5, it can be seen that in the Parkinsonian state, the neuron population will have multiple peaks with 22 Hz as the oscillation period, and it is more appropriate to choose a value close to half of the oscillation period as the time delay in the delayed feedback control algorithm.Therefore, the delay η = 12 ms was firstly set, however, when investigating scheme A and scheme C, it was found that the oscillation control effect was not satisfactory when η = 12 ms, so the delays of scheme A and scheme C were adjusted to 17 ms and 6 ms, respectively.So, the four delay feedback strategies were finally set to 17 ms, 12 ms, 6 ms, and 12 ms, respectively.The comprehensive effect of time delay and feedback gain G with respect to the control signal will be investigated in Sect.3.4.
The control strategy was switched at t = 1 s.It can be seen from Fig. 7 that when a delayed feedback control strategy based on STN or GPe was implemented in the system, the pathological oscillation suppression of neuronal populations can be successfully achieved.As many populations of neurons were present in the neural mass model, we only showed the suppression effect of Th and PPN in addition to STN and GPe.It can also be found from the spectrum analysis that the neuronal populations change from multiple peaks to normal power after the control signal was applied.

The role of feedback gain G in four linear delayed feedback stimulation schemes
To investigate the effect of the feedback gain G on pathological oscillations for the four different schemes, we calculated the minimum and maximum output values of the neuronal populations in the last second and gave the bifurcation diagrams of the STN, GPe, Th, and PPN under different strategies as shown in Fig. 8.The part of the graph where the minimum and maximum values overlap represents the effective suppression of oscillations.By testing, we choose the value of feedback gain G between −10 and 10 to explore.In this range, since STN has excitatory projections to other neuronal populations, positive G is valid for m STN of scheme A and scheme D, and since GPe has inhibi- tory projections to other neuronal populations, negative G is valid for m GPe of scheme B and scheme C [26].It means that the valid range exists for scheme A and scheme D only when G is at a positive value, while scheme B and scheme C exist only when G is at a negative value.From Fig. 8, it can be found that the oscillation can be effec- tively suppressed only when G is in a reasonable range, while it cannot be suppressed in  other ranges, and even enhanced oscillation can occur.And there were different effective suppression ranges for the four control schemes.The results showed that the effective G ranges corresponding to the delayed feedback control schemes A and D are 0.303 to 0.404 and 0.303 to 0.960, respectively, and the effective G ranges corresponding to schemes B and C are -0.253 to -0.909 and -0.606 to -3.030, respectively.Therefore, compared with the first conventional scheme A, the feedback gain G of the latter three schemes has a larger effective range, and the scheme C has the widest range of feedback gain G.
The average values of the oscillation index I STN os , I GPe os , I Th os , and I PPN os for STN, GPe, Th, and PPN in the four delayed feedback strategies are given in Fig. 9.It can be seen that the effective region of the feedback gain G shown in Fig. 8 is also shown in Fig. 9. Therefore, it showed that the three strategies B, C, and D are more robust to the range of variation of G compared to the conventional strategy A. Meanwhile, Fig. 9 also gives the energy con- sumption corresponding to the four schemes after their implementation.As can be seen, the energy consumption of strategies B, C, and D was slightly reduced compared to that of scheme A. Thus, it was demonstrated that schemes B, C, and D have a wider range of feedback gain and can be performed without increasing the energy consumption.

The exploration of the parameter space for time delay η and feedback gain G
In the linear delay feedback strategy, time delay also plays a role that should not be underestimated [52].The combined effect of time delay and feedback gain G on the four schemes is given in Fig. 10.Only the effects of and G on the STN and GPe nucleus were shown here.The first two columns of Fig. 10 represent the combined effect of and G on the oscillation index I STN OS and I GPe OS of the two nuclei.It can be seen that there are several regions in each plot where oscillations are suppressed, and we called them oscillation suppression regions.The color of the regions represented the change of I STN OS and I GPe OS under the combined influence of and G .And the white color outside the region indicated the area where the oscillations expand in the pathological state.
To explore the role of time delay, we considered a range of several oscillation periods.It can be seen that three oscillation suppression regions appear with increasing time delay, but the center of the time delay differs for each region.It can be found that scheme A has a center point corresponding to 17 ms, while schemes B, C, and D have centers of 12 ms, 6 ms, and 12 ms, respectively, probably because the GPe was self-repressed compared to the STN, thus reducing the time of oscillation suppression.And it can be found that scheme C has the shortest time delay, which means that when both the stimulus nucleus and the feedback nucleus are GPe, the propagation speed is faster than the other three schemes, which can directly and quickly affect the oscillation suppression.From the first two columns of Fig. 10, it can be observed that the area of the first oscillation suppression region of scheme B, C, and D was much larger than that of the oscillation suppression region of scheme A. In other words, the last three schemes had more robust oscillation suppression regions and better suppression, which indicates that schemes B, C, and D seem to have better suppression than scheme A and perhaps have better robustness in applications.
The third column in Fig. 10 represents the combined effect of and G on the energy consumption I EN of the four schemes.The oscillation suppression regions in Fig. 10 rep- resent the regions where the energy consumption values used for oscillation suppression It is evident that the oscillation suppression regions of schemes B, C, and D also have a much larger area than scheme A. In other words, the latter three schemes do not increase the energy consumption compared with scheme A, and even reduce it while expanding the effective suppression space.Therefore, it was again verified that schemes B, C, and D have better robustness.
In the recent work of Popovych et al. [22,23,25], by obtaining feedback signals from the LFP of STN neurons and applying stimulation to STN neurons in STN-GPe networks found that both smooth and pulsatile multisite linear delay feedback (MLDF) stimulation could effectively remove abnormal synchronization of STN neurons, and found that adaptive DBS (aDBS) and adaptive pulsatile linear delayed feedback stimulation (apLDF) were more effective than conventional high-frequency DBS (cDBS) in suppressing abnormal synchronization.Based on a comprehensive exploration of four protocols containing either the STN or GPe, we found that the strategy of using the GPe as a source of feedback signal or as a stimulus target was more robust to changes of delay and gain parameters.Crompe et al. [56] experimentally found that GPe activity has an important influence on the expression and propagation of β oscillations in PD.They demonstrated that GPe is well able to coordinate the firing frequency and synchronization changes of cortex-basal ganglia (CBG) circuits in both normal and Parkinsonian states, emphasizing the key role played by GPe.Xie et al. [57] found that GPe-DBS could significantly improve the electrical activity of GPi through rat experiments, indicating that GPe may have an important impact on the treatment of PD.Our computational results validate the important influence of GPe on abnormal neural dynamics in PD and also suggest that delayed feedback scheme that includes only STN is more sensitive to parameter changes.Therefore, our computational results provide evidence that the subsequent application of the DBS strategy should consider not only the role of STN but also the role of GPe appropriately, which may be able to reduce energy consumption and thus effectively achieve improved control performance.

Conclusions
To enhance the control performance of linear delayed feedback strategy, a larger range of parameters needs to be found.Given the special role of PPN, the neural mass model of the closed-loop system of cortex-basal ganglia-thalamus-pedunculopontine nucleus was firstly proposed in this paper to simulate the real brain dynamic characteristics.Based on the exploration of the coupling strengths between and within the neuronal populations in the model, we reproduced the oscillations of the model in normal and pathological states.Second, the special role of the GPe nucleus was taken into account by adding it to the linear delayed feedback schemes.In other words, three stimulation strategies were investigated using GPe as stimulus nucleus and feedback nucleus.And it is compared with the conventional strategy of using only the STN as the stimulus nucleus and the feedback nucleus.The simulation results show that all three GPe-based schemes can effectively suppress the oscillation of frequency in PD state.And compared with the conventional scheme, the GPe-based scheme can not only expand the parameter space with reduced time delay, but also effectively reduce the energy consumption.This shows that the GPe-based delayed feedback control scheme has high robustness.GPe may play an important role in future clinical applications.However, our conclusions need to be further validated by clinical trials.We hope our findings can bring some insight into the treatment of PD.

Fig. 1 a
Fig. 1 a Structure of the CBGTh-PPN loop model.It consists of three parts: the cortex, the BGTh, and the PPN.The cortex contains three populations of neurons, namely EI, PY, and II.The BGTh contains seven populations of neurons, namely FSI, MSN D1 , MSN D2 , STN, GPe, GPi, and Th.The red lines with triangular arrows indicate glutamate-mediated excitatory projections, and the green lines with circles indicate inhibitory projections mediated by GABA receptors.b A sketch of the CBGTh-PPN model

Fig. 2
Fig.2The neural mass model of CBGTh-PPN.Each of these neurons consists of a linear part T i and a nonlinear part block S i

4 dtFig. 3
Fig. 3 Four different delayed feedback schemes.a-d Scheme A, B, C, and D, respectively.a Scheme A. Delayed feedback signals are obtained from the LFPs of the STN and received by the STN.b Scheme B.The feedback signal is generated from the LFPs of the STN and the simulation is delivered to the GPe.c Scheme C. Similar to scheme A, the delayed feedback signal is obtained from the LFPs of the GPe and received by the GPe.d Scheme D. Contrary to scheme B, the feedback signal is generated from the LFPs of the GPe and the simulation is delivered to the STN

Fig. 4
Fig.4 Combined effect of coupling strengths K 1 , K 2 , and K 3 on the peak power and peak frequency of PY, STN, and PPN populations.a, b Effect of coupling strengths K 1 and K 2 on the peak power and peak frequency of PY and STN neuron populations.c The effect of coupling strengths K 2 and K 3 on the peak power and peak frequency of PPN.d The effect of K 1 and K 3 on the peak power and peak frequency of PPN

Fig. 5
Fig. 5 Corresponding output and power spectral density of neuronal populations in neural circuits in normal and pathological states.a The output of all populations in the CBGTh-PPN loop in normal and pathological states.b The power spectral densities of neuronal populations in BGTh and PPN in two states.The phase trajectories of the four populations in STN, Th, GPe, and PPN are shown in c.The orange solid line in a-c represents the normal state, and the blue solid line represents the Parkinsonian state

Fig. 6
Fig. 6 Stimulus signals for the four delayed feedback control strategies.They correspond to the four schemes in Fig. 3.All four stimulus signals fluctuate up and down around the zero point and have different waveforms

Fig. 7
Fig. 7 The suppression effects after applying the delayed feedback control schemes.Each scheme demonstrates the suppression effect on STN, GPe, Th, and PPN nucleus.And the feedback gain G of the four strategies A, B, C, and D are 0.4, -0.3, -0.8, 0.4, respectively, and the corresponding time delay is 17 ms, 12 ms, 6 ms, and 12 ms, respectively

Fig. 8
Fig. 8 The effect of feedback gain G on the four delayed feedback schemes.The time delays corresponding to the four schemes A, B, C, and D are 17 ms, 12 ms, 6 ms, 12 ms, respectively.All schemes are implemented on the nucleus at t = 1 s.The bifurcation diagrams of the four populations under the four strategies are obtained by calculating the minimum values of the maximum values of the STN, GPe, Th, and PPN in the last second, and thus, the effective ranges of the feedback gains corresponding to the four schemes are 0.303 to 0.404, -0.253 to -0.909, -0.606 to -3.030, and 0.303 to 0.960, respectively

Fig. 9
Fig. 9 The average values of the oscillation indexes I STN os , I GPe os , I Th os , and I PPN os for the STN, GPe, Th, and PPN are given.And the energy consumption I EN of each of the four schemes is given in the bottom panels

Fig. 10
Fig. 10 Oscillation suppression regions.It reflects the combined influences of time delay and feedback gain on the oscillation index and energy consumption of the four schemes.The first two columns indicate the effect of and G on the oscillation index I STN OS and I GPe OS of STN and GPe, and the last column indicates the effect of and G on the energy consumption I EN of STN and GPe